epsilon input file support (NanoComp#240)

* Add epsilon_input_file support and test * Make all_freqs an ndarray
bencbartlett · Mar 22, 2018 · 095bdfd · 095bdfd
1 parent aff0d33
commit 095bdfd
Show file tree

Hide file tree

Showing 4 changed files with 56 additions and 17 deletions.
diff --git a/libpympb/pympb.cpp b/libpympb/pympb.cpp
@@ -332,14 +332,12 @@ int mode_solver::mean_epsilon(symmetric_matrix* meps, symmetric_matrix *meps_inv
     q.y = p.y + neighbors[dimensions - 1][i][1] * d2;
     q.z = p.z + neighbors[dimensions - 1][i][2] * d3;
 
-    // TODO: shift_to_unit_cell runs forever if size has a 0 component
     geometry_lattice.size.x = geometry_lattice.size.x == 0 ? 1e-20 : geometry_lattice.size.x;
     geometry_lattice.size.y = geometry_lattice.size.y == 0 ? 1e-20 : geometry_lattice.size.y;
     geometry_lattice.size.z = geometry_lattice.size.z == 0 ? 1e-20 : geometry_lattice.size.z;
 
     z = shift_to_unit_cell(q);
 
-    // TODO
     geometry_lattice.size.x = geometry_lattice.size.x == 1e-20 ? 0 : geometry_lattice.size.x;
     geometry_lattice.size.y = geometry_lattice.size.y == 1e-20 ? 0 : geometry_lattice.size.y;
     geometry_lattice.size.z = geometry_lattice.size.z == 1e-20 ? 0 : geometry_lattice.size.z;
@@ -847,9 +845,13 @@ void mode_solver::reset_epsilon() {
     (dimensions > 2) ? mesh_size : 1,
   };
 
-  // TODO: Support epsilon_input_file
-  // get_epsilon_file_func(epsilon_input_file, &d.epsilon_file_func, &d.epsilon_file_func_data);
-  // get_epsilon_file_func(mu_input_file, &d.mu_file_func, &d.mu_file_func_data);
+  if (!epsilon_input_file.empty()) {
+      default_material = meep_geom::make_file_material(epsilon_input_file.c_str());
+  }
+
+  // TODO: support mu_input_file
+  // if (!mu_input_file.empty()) {
+  // }
 
   meep::master_printf("Initializing epsilon function...\n");
   set_maxwell_dielectric(mdata, mesh, R, G, dielectric_function, mean_epsilon_func,
@@ -862,8 +864,6 @@ void mode_solver::reset_epsilon() {
                    static_cast<void *>(this));
     eps = true;
   }
-  // destroy_epsilon_file_func_data(d.epsilon_file_func_data);
-  // destroy_epsilon_file_func_data(d.mu_file_func_data);
 }
 
 bool mode_solver::has_mu() {
@@ -1132,12 +1132,6 @@ void mode_solver::solve_kpoint(vector3 kvector) {
     free(deflation.S);
   }
 
-  // TODO
-  // if (num_write_output_vars > 0) {
-  //   //  clean up from prev. call
-  //   destroy_output_vars();
-  // }
-
   iterations = total_iters; /* iterations output variable */
 
   set_kpoint_index(kpoint_index + 1);

diff --git a/python/solver.py b/python/solver.py
@@ -83,7 +83,7 @@ def __init__(self,
         self.verbose = verbose
         self.parity = ''
         self.iterations = 0
-        self.all_freqs = []
+        self.all_freqs = None
         self.freqs = []
         self.band_range_data = []
         self.eigensolver_flops = 0
@@ -596,7 +596,7 @@ def run_parity(self, p, reset_fields, *band_functions):
 
         start = time.time()
 
-        self.all_freqs = []
+        self.all_freqs = np.zeros((len(self.k_points), self.num_bands))
         self.band_range_data = []
 
         init_time = time.time()
@@ -651,14 +651,14 @@ def run_parity(self, p, reset_fields, *band_functions):
                 self.output_mu()  # and mu too, if we have it
 
         if self.num_bands > 0:
-            for k in k_split[1]:
+            for i, k in enumerate(k_split[1]):
                 self.current_k = k
                 solve_kpoint_time = time.time()
                 self.mode_solver.solve_kpoint(k)
                 self.iterations = self.mode_solver.get_iterations()
                 print("elapsed time for k point: {}".format(time.time() - solve_kpoint_time))
                 self.freqs = self.get_freqs()
-                self.all_freqs.append(self.freqs)
+                self.all_freqs[i, :] = np.array(self.freqs)
                 self.band_range_data = self.update_band_range_data(self.band_range_data,
                                                                    self.freqs, k)
                 self.eigensolver_iters += [self.iterations / self.num_bands]

diff --git a/python/tests/data/eps_input_file_test.h5 b/python/tests/data/eps_input_file_test.h5
diff --git a/python/tests/mpb.py b/python/tests/mpb.py
@@ -1130,6 +1130,51 @@ def get_efields(ms, band):
 
         self.check_fields_against_h5(ref_path, result.ravel(), suffix='-new')
 
+    def test_epsilon_input_file(self):
+        ms = self.init_solver(geom=False)
+        eps_fn = 'eps_input_file_test.h5'
+        ms.epsilon_input_file = os.path.join(self.data_dir, eps_fn)
+
+        ms.run_te()
+
+        expected_freqs = np.array([
+            0.0, 0.5543986136451342, 0.7613327775255415, 0.7613339178956054,
+            0.8940893915924257, 0.998342969572652, 0.9983441882455961, 1.0747466061007138
+        ])
+
+        expected_gap_list = [
+            (3.848610367089048e-5, 0.5781812856814899, 0.5781815082009817),
+            (1.4651880980150234, 0.8051999839699242, 0.8170847453549156),
+            (0.75255857475812, 1.0580309832489785, 1.0660233597945266),
+        ]
+
+        expected_brd = [
+            ((0.0, mp.Vector3(0.0, 0.0, 0.0)),
+             (0.4970977843772992, mp.Vector3(0.5, 0.5, 0.0))),
+            ((0.4402896410505961, mp.Vector3(0.5, 0.0, 0.0)),
+             (0.5781812856814899, mp.Vector3(0.5, 0.5, 0.0))),
+            ((0.5781815082009817, mp.Vector3(0.5, 0.5, 0.0)),
+             (0.761332777525562, mp.Vector3(0.0, 0.0, 0.0))),
+            ((0.6689126424359774, mp.Vector3(0.5, 0.5, 0.0)),
+             (0.8051999839699242, mp.Vector3(0.30000000000000004, 0.30000000000000004, 0.0))),
+            ((0.8170847453549156, mp.Vector3(0.30000000000000004, 0.30000000000000004, 0.0)),
+             (0.8940893915924548, mp.Vector3(0.0, 0.0, 0.0))),
+            ((0.8826671164993868, mp.Vector3(0.30000000000000004, 0.30000000000000004, 0.0)),
+             (1.0014926328155058, mp.Vector3(0.5, 0.0, 0.0))),
+            ((0.8832199143682116, mp.Vector3(0.5, 0.5, 0.0)),
+             (1.0580309832489785, mp.Vector3(0.5, 0.0, 0.0))),
+            ((1.0660233597945266, mp.Vector3(0.2, 0.2, 0.0)),
+             (1.087345829555555, mp.Vector3(0.5, 0.5, 0.0))),
+        ]
+
+        self.check_band_range_data(expected_brd, ms.band_range_data)
+        self.compare_arrays(expected_freqs, ms.all_freqs[-1])
+
+        self.check_gap_list(expected_gap_list, ms.gap_list)
+
+        pt = ms.get_epsilon_point(mp.Vector3(0.5, 0.5))
+        self.assertEqual(pt, 1.0)
+
 
 if __name__ == '__main__':
     unittest.main()