Enable Eigenmode Features with Dispersive Materials

python/Makefile.am

@@ -18,12 +18,14 @@ endif # WITH_MPI
 
 if WITH_MPB
   BINARY_GRATING_TEST = $(TEST_DIR)/binary_grating.py
+  DISPERSIVE_EIGENMODE_TEST = $(TEST_DIR)/dispersive_eigenmode.py
   KDOM_TEST = $(TEST_DIR)/kdom.py
   MODE_COEFFS_TEST = $(TEST_DIR)/mode_coeffs.py
   MODE_DECOMPOSITION_TEST = $(TEST_DIR)/mode_decomposition.py
   WVG_SRC_TEST = $(TEST_DIR)/wvg_src.py
 else
   BINARY_GRATING_TEST =
+  DISPERSIVE_EIGENMODE_TEST =
   KDOM_TEST =
   MODE_COEFFS_TEST =
   MODE_DECOMPOSITION_TEST =
@@ -41,6 +43,7 @@ TESTS =                                   \
     $(TEST_DIR)/cavity_farfield.py        \
     $(TEST_DIR)/chunks.py                 \
     $(TEST_DIR)/cyl_ellipsoid.py          \
+    ${DISPERSIVE_EIGENMODE_TEST}          \
     $(TEST_DIR)/dft_energy.py             \
     $(TEST_DIR)/dft_fields.py             \
     $(TEST_DIR)/faraday_rotation.py       \

python/geom.py

@@ -177,7 +177,7 @@ def __init__(self, epsilon_diag=Vector3(1, 1, 1),
                  E_chi3=None,
                  H_chi2=None,
                  H_chi3=None,
-                 valid_freq_range=None):
+                 valid_freq_range=FreqRange(min=-mp.inf, max=mp.inf)):
 
         if epsilon:
             epsilon_diag = Vector3(epsilon, epsilon, epsilon)

python/simulation.py

@@ -1264,9 +1264,9 @@ def get_field_point(self, c, pt):
         v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical)
         return self.fields.get_field_from_comp(c, v3)
 
-    def get_epsilon_point(self, pt):
+    def get_epsilon_point(self, pt, omega = 0):
         v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical)
-        return self.fields.get_eps(v3)
+        return self.fields.get_eps(v3,omega)
 
     def get_filename_prefix(self):
         if isinstance(self.filename_prefix, str):
@@ -1795,15 +1795,15 @@ def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args):
 
         return stuff
 
-    def output_component(self, c, h5file=None):
+    def output_component(self, c, h5file=None, omega=0):
         if self.fields is None:
             raise RuntimeError("Fields must be initialized before calling output_component")
 
         vol = self.fields.total_volume() if self.output_volume is None else self.output_volume
         h5 = self.output_append_h5 if h5file is None else h5file
         append = h5file is None and self.output_append_h5 is not None
 
-        self.fields.output_hdf5(c, vol, h5, append, self.output_single_precision, self.get_filename_prefix())
+        self.fields.output_hdf5(c, vol, h5, append, self.output_single_precision,self.get_filename_prefix(),omega)
 
         if h5file is None:
             nm = self.fields.h5file_name(mp.component_name(c), self.get_filename_prefix(), True)
@@ -1833,7 +1833,7 @@ def h5topng(self, rm_h5, option, *step_funcs):
         cmd = re.sub(r'\$EPS', self.last_eps_filename, opts)
         return convert_h5(rm_h5, cmd, *step_funcs)
 
-    def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None, arr=None):
+    def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None, arr=None, omega = 0):
         if component is None:
             raise ValueError("component is required")
         if isinstance(component, mp.Volume) or isinstance(component, mp.volume):
@@ -1868,9 +1868,9 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None
             arr = np.zeros(dims, dtype=np.complex128 if cmplx else np.float64)
 
         if np.iscomplexobj(arr):
-            self.fields.get_complex_array_slice(v, component, arr)
+            self.fields.get_complex_array_slice(v, component, arr, omega)
         else:
-            self.fields.get_array_slice(v, component, arr)
+            self.fields.get_array_slice(v, component, arr, omega)
 
         return arr
 
@@ -2071,8 +2071,8 @@ def run(self, *step_funcs, **kwargs):
         else:
             raise ValueError("Invalid run configuration")
 
-    def get_epsilon(self):
-        return self.get_array(component=mp.Dielectric)
+    def get_epsilon(self,omega=0):
+        return self.get_array(component=mp.Dielectric,omega=omega)
 
     def get_mu(self):
         return self.get_array(component=mp.Permeability)
@@ -2600,12 +2600,14 @@ def _output_png(sim, todo):
     return _output_png
 
 
-def output_epsilon(sim):
-    sim.output_component(mp.Dielectric)
+def output_epsilon(sim,*step_func_args,**kwargs):
+    omega = kwargs.pop('omega', 0.0)
+    sim.output_component(mp.Dielectric,omega=omega)
 
 
-def output_mu(sim):
-    sim.output_component(mp.Permeability)
+def output_mu(sim,*step_func_args,**kwargs):
+    omega = kwargs.pop('omega', 0.0)
+    sim.output_component(mp.Permeability,omega=omega)
 
 
 def output_hpwr(sim):

python/tests/dispersive_eigenmode.py

@@ -0,0 +1,168 @@
+
+# dispersive_eigenmode.py - Tests the meep eigenmode features (eigenmode source,
+# eigenmode decomposition, and get_eigenmode) with dispersive materials.
+# TODO: 
+#  * check materials with off diagonal components
+#  * check magnetic profiles
+#  * once imaginary component is supported, check that
+
+from __future__ import division
+
+import unittest
+import meep as mp
+import numpy as np
+from meep import mpb
+import h5py
+import os
+
+
+class TestDispersiveEigenmode(unittest.TestCase):
+
+    # Directly calss the C++ chi1 routine
+    def call_chi1(self,material,component,direction,omega):
+
+        sim = mp.Simulation(cell_size=mp.Vector3(1,1,1),
+                    default_material=material,
+                    resolution=10)
+
+        sim.init_sim()
+        v3 = mp.py_v3_to_vec(sim.dimensions, mp.Vector3(0,0,0), sim.is_cylindrical)
+        n = 1/np.sqrt(sim.structure.get_chi1inv(int(component),int(direction),v3,omega))
+        return n
+
+    # Pulls the "effective index" of a uniform, dispersive material
+    # (i.e. the refractive index) using meep's get_eigenmode
+    def simulate_meep(self,material,omega):
+
+        sim = mp.Simulation(cell_size=mp.Vector3(2,2,2),
+                            default_material=material,
+                            resolution=20
+                            )
+
+        direction = mp.X
+        where = mp.Volume(center=mp.Vector3(0,0,0),size=mp.Vector3(0,1,1))
+        band_num = 1
+        kpoint = mp.Vector3(2,0,0)
+        sim.init_sim()
+        em = sim.get_eigenmode(omega,direction,where,band_num,kpoint)
+        neff_meep = np.squeeze(em.k.x) / np.squeeze(em.freq)
+
+        return neff_meep
+
+    # Pulls the "effective index" of a uniform, dispersive material
+    # (i.e. the refractive index) using mpb
+    def simulate_mpb(self,material,omega):
+        ms = mpb.ModeSolver(
+            geometry_lattice=mp.Lattice(size=mp.Vector3(0,2,2)),
+            default_material=material,
+            resolution=10,
+            num_bands=1
+        )
+        k = ms.find_k(mp.NO_PARITY, omega, 1, 1, mp.Vector3(1), 1e-3, omega * 1,
+                omega * 0.1, omega * 6)
+
+        neff_mpb = k[0]/omega
+        return neff_mpb
+
+    # main test bed to check the new features
+    def compare_meep_mpb(self,material,omega,component=0,direction=0):
+        n = np.real(np.sqrt(material.epsilon(omega)[component,direction]))
+        chi1 = self.call_chi1(material,mp.Ex,mp.X,omega)
+        n_meep = self.simulate_meep(material,omega)
+        # Let's wait to check this until we enable dispersive materials in MPB...
+        #n_mpb = self.simulate_mpb(material,omega)
+
+        # Check that the chi1 value matches the refractive index
+        self.assertAlmostEqual(n,chi1, places=6)
+
+        # Check that the chi1 value matches meep's get_eigenmode
+        self.assertAlmostEqual(n,n_meep, places=6)
+
+        # Check that the chi1 value matches mpb's get_eigenmode 
+        #self.assertAlmostEqual(n,n_mpb, places=6)
+
+    def test_chi1_routine(self):
+        # This test checks the newly implemented get_chi1inv routines within the 
+        # fields and structure classes by comparing their output to the 
+        # python epsilon output.
+
+        from meep.materials import Si, Ag, LiNbO3, Au
+
+        # Check Silicon
+        w0 = Si.valid_freq_range.min
+        w1 = Si.valid_freq_range.max
+        self.compare_meep_mpb(Si,w0)
+        self.compare_meep_mpb(Si,w1)
+
+        # Check Silver
+        w0 = Ag.valid_freq_range.min
+        w1 = Ag.valid_freq_range.max
+        self.compare_meep_mpb(Ag,w0)
+        self.compare_meep_mpb(Ag,w1)
+
+        # Check Gold
+        w0 = Au.valid_freq_range.min
+        w1 = Au.valid_freq_range.max
+        self.compare_meep_mpb(Au,w0)
+        self.compare_meep_mpb(Au,w1)   
+
+        # Check Lithium Niobate (X,X)
+        w0 = LiNbO3.valid_freq_range.min
+        w1 = LiNbO3.valid_freq_range.max
+        self.compare_meep_mpb(LiNbO3,w0)
+        self.compare_meep_mpb(LiNbO3,w1)
+
+    def verify_output_and_slice(self,material,omega,component=0,direction=0):
+        filename = 'dispersive_eigenmode-eps-000000.00.h5'
+        n = np.real(np.sqrt(material.epsilon(omega)[component,direction]))
+
+        sim = mp.Simulation(cell_size=mp.Vector3(2,2,2),
+                            default_material=material,
+                            resolution=20,
+                            eps_averaging=False
+                            )
+        sim.init_sim()
+
+        # Check to make sure the get_slice routine is working with omega
+        n_slice = np.sqrt(np.min(sim.get_epsilon(omega)))
+        self.assertAlmostEqual(n,n_slice, places=6)
+
+        # Check to make sure h5 output is working with omega
+        mp.output_epsilon(sim,omega=omega)
+        n_h5 = np.sqrt(np.min(h5py.File(filename, 'r')['eps']))
+        self.assertAlmostEqual(n,n_h5, places=6)
+        os.remove(filename)
+
+    def test_get_with_dispersion(self):
+        # This test checks the get_array_slice and output_fields method
+        # with dispersive materials.
+
+        from meep.materials import Si, Ag, LiNbO3, Au
+
+        # Check Silicon
+        w0 = Si.valid_freq_range.min
+        w1 = Si.valid_freq_range.max
+        self.verify_output_and_slice(Si,w0)
+        self.verify_output_and_slice(Si,w1)
+
+        # Check Silver
+        w0 = Ag.valid_freq_range.min
+        w1 = Ag.valid_freq_range.max
+        self.verify_output_and_slice(Ag,w0)
+        self.verify_output_and_slice(Ag,w1)
+
+        # Check Gold
+        w0 = Au.valid_freq_range.min
+        w1 = Au.valid_freq_range.max
+        self.verify_output_and_slice(Au,w0)
+        self.verify_output_and_slice(Au,w1)       
+
+        # Check Lithium Niobate (X,X)
+        w0 = LiNbO3.valid_freq_range.min
+        w1 = LiNbO3.valid_freq_range.max
+        #self.verify_output_and_slice(LiNbO3,w0)
+        #self.verify_output_and_slice(LiNbO3,w1)
+
+
+if __name__ == '__main__':
+    unittest.main()

python/tests/visualization.py

@@ -31,40 +31,50 @@ def setup_sim(zDim=0):
 
     # A simple waveguide
     geometry = [mp.Block(mp.Vector3(mp.inf,1,1),
-                     center=mp.Vector3(),
-                     material=mp.Medium(epsilon=12))]
-    
+                        center=mp.Vector3(),
+                        material=mp.Medium(epsilon=12))]
+
     # Add point sources
     sources = [mp.Source(mp.ContinuousSource(frequency=0.15),
                         component=mp.Ez,
-                        center=mp.Vector3(-5,0)),
+                        center=mp.Vector3(-5,0),
+                        size=mp.Vector3(0,0,2)),
+                mp.Source(mp.ContinuousSource(frequency=0.15),
+                        component=mp.Ez,
+                        center=mp.Vector3(0,2),
+                        size=mp.Vector3(0,0,2)),
+                mp.Source(mp.ContinuousSource(frequency=0.15),
+                        component=mp.Ez,
+                        center=mp.Vector3(-1,1),
+                        size=mp.Vector3(0,0,2)),
                 mp.Source(mp.ContinuousSource(frequency=0.15),
                         component=mp.Ez,
-                        center=mp.Vector3(0,2))
+                        center=mp.Vector3(-2,-2,1),
+                        size=mp.Vector3(0,0,0)), 
                         ]
 
     # Add line sources
     sources += [mp.Source(mp.ContinuousSource(frequency=0.15),
                         component=mp.Ez,
-                        size=mp.Vector3(0,2,0),
+                        size=mp.Vector3(0,2,2),
                         center=mp.Vector3(-6,0)),
                 mp.Source(mp.ContinuousSource(frequency=0.15),
                         component=mp.Ez,
-                        size=mp.Vector3(2,0,0),
+                        size=mp.Vector3(0,2,2),
                         center=mp.Vector3(0,1))]
 
     # Add plane sources
     sources += [mp.Source(mp.ContinuousSource(frequency=0.15),
                         component=mp.Ez,
-                        size=mp.Vector3(2,2,0),
+                        size=mp.Vector3(2,2,2),
                         center=mp.Vector3(-3,0)),
                 mp.Source(mp.ContinuousSource(frequency=0.15),
                         component=mp.Ez,
-                        size=mp.Vector3(2,2,0),
+                        size=mp.Vector3(2,2,2),
                     center=mp.Vector3(0,-2))]
-    
+
     # Different pml layers
-    pml_layers = [mp.PML(2.0,mp.X),mp.PML(1.0,mp.Y,mp.Low),mp.PML(1.5,mp.Y,mp.High)]
+    pml_layers = [mp.PML(2.0,mp.X),mp.PML(1.0,mp.Y,mp.Low),mp.PML(1.5,mp.Y,mp.High),mp.PML(1.5,mp.Z)]
 
     resolution = 10
 
@@ -74,13 +84,33 @@ def setup_sim(zDim=0):
                     sources=sources,
                     resolution=resolution)
     # Line monitor
-    sim.add_flux(1,0,1,mp.FluxRegion(center=mp.Vector3(5,0,0),size=mp.Vector3(0,4), direction=mp.X))
+    sim.add_flux(1,0,1,mp.FluxRegion(center=mp.Vector3(5,0,0),size=mp.Vector3(0,4,4), direction=mp.X))
 
     # Plane monitor
-    sim.add_flux(1,0,1,mp.FluxRegion(center=mp.Vector3(2,0,0),size=mp.Vector3(4,4), direction=mp.X))
-    
+    sim.add_flux(1,0,1,mp.FluxRegion(center=mp.Vector3(2,0,0),size=mp.Vector3(4,4,4), direction=mp.X))
+
     return sim
 
+def view_sim():
+    sim = setup_sim(8)
+    xy0 = mp.Volume(center=mp.Vector3(0,0,0), size=mp.Vector3(sim.cell_size.x,sim.cell_size.y,0))
+    xy1 = mp.Volume(center=mp.Vector3(0,0,1), size=mp.Vector3(sim.cell_size.x,sim.cell_size.y,0))
+    yz0 = mp.Volume(center=mp.Vector3(0,0,0), size=mp.Vector3(0,sim.cell_size.y,sim.cell_size.z))
+    yz1 = mp.Volume(center=mp.Vector3(1,0,0), size=mp.Vector3(0,sim.cell_size.y,sim.cell_size.z))
+    xz0 = mp.Volume(center=mp.Vector3(0,0,0), size=mp.Vector3(sim.cell_size.x,0,sim.cell_size.z))
+    xz1 = mp.Volume(center=mp.Vector3(0,1,0), size=mp.Vector3(sim.cell_size.x,0,sim.cell_size.z))
+    vols = [xy0,xy1,yz0,yz1,xz0,xz1]
+    titles = ['xy0','xy1','yz0','yz1','xz0','xz1']
+    xlabel = ['x','x','y','y','x','x']
+    ylabel = ['y','y','z','z','z','z']
+    for k in range(len(vols)):
+        ax = plt.subplot(2,3,k+1)
+        sim.plot2D(ax=ax,output_plane=vols[k])
+        ax.set_xlabel(xlabel[k])
+        ax.set_ylabel(ylabel[k])
+        ax.set_title(titles[k])
+    plt.tight_layout()
+    plt.show()
 class TestVisualization(unittest.TestCase):
 
     def test_plot2D(self):