WIP: mode-solver improvements (#333)

* mpb match-freq should only search k in direction d

* adjust for bloch-periodic boundaries in get_eigenmode, so that it works for both sources and mode decomposition

* allow get_eigenmode to return NULL if mode was not found, clean up a memory leak in mode coefficients

* whoops

* add parity argument to get_eigenmode_coefficients (fises #294)

* add missing get_eigenmode_coefficients parameters, fix defaults

doc/docs/Python_Tutorials/Mode_Decomposition.md

@@ -26,25 +26,25 @@ def main(args):
     Lt = args.Lt      # taper length
 
     Si = mp.Medium(epsilon=12.0)
-    
+
     dair = 3.0
     dpml = 3.0
-    
+
     sx = dpml + Lw + Lt + Lw + dpml
     sy = dpml + dair + w2 + dair + dpml
     cell_size = mp.Vector3(sx, sy, 0)
 
-    geometry = [ mp.Block(material=Si, center=mp.Vector3(0,0,0), size=mp.Vector3(mp.inf,w1,mp.inf)), 
+    geometry = [ mp.Block(material=Si, center=mp.Vector3(0,0,0), size=mp.Vector3(mp.inf,w1,mp.inf)),
                  mp.Block(material=Si, center=mp.Vector3(0.5*sx-0.5*(Lt+Lw+dpml),0,0), size=mp.Vector3(Lt+Lw+dpml,w2,mp.inf)) ]
 
     # form linear taper
 
     hh = w2
     ww = 2*Lt
-    
+
     # taper angle (CCW, relative to +X axis)
     rot_theta = math.atan(0.5*(w2-w1)/Lt)
-    
+
     pvec = mp.Vector3(-0.5*sx+dpml+Lw,0.5*w1,0)
     cvec = mp.Vector3(-0.5*sx+dpml+Lw+0.5*ww,0.5*hh+0.5*w1,0)
     rvec = cvec-pvec
@@ -69,7 +69,7 @@ def main(args):
 
     # mode frequency
     fcen = 0.15
-    
+
     sources = [ mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=0.5*fcen),
                                    component=mp.Ez,
                                    size=mp.Vector3(0,sy-2*dpml,0),
@@ -78,9 +78,9 @@ def main(args):
                                    eig_parity=mp.ODD_Z,
                                    eig_kpoint=mp.Vector3(0.4,0,0),
                                    eig_resolution=32) ]
-    
+
     symmetries = [ mp.Mirror(mp.Y,+1) ]
-    
+
     sim = mp.Simulation(resolution=resolution,
                         cell_size=cell_size,
                         boundary_layers=boundary_layers,
@@ -90,21 +90,21 @@ def main(args):
 
     xm = -0.5*sx + dpml + 0.5*Lw; # x-coordinate of monitor
     mflux = sim.add_eigenmode(fcen, 0, 1, mp.FluxRegion(center=mp.Vector3(xm,0), size=mp.Vector3(0,sy-2*dpml)));
-        
+
     sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(xm,0), 1e-10))
 
     bands = np.array([1],dtype=np.int32) # indices of modes for which to compute expansion coefficients
     num_bands = 1;
     alpha = np.zeros(2*num_bands,dtype=np.complex128) # preallocate array to store coefficients
     vgrp = np.zeros(num_bands,dtype=np.float64)       # also store mode group velocities
     mvol = mp.volume(mp.vec(xm,-0.5*sy+dpml),mp.vec(xm,+0.5*sy-dpml))
-    sim.fields.get_eigenmode_coefficients(mflux, mp.X, mvol, bands, alpha, vgrp)
+    sim.fields.get_eigenmode_coefficients(mflux, mp.X, mvol, bands, mp.ODD_Z, alpha, vgrp)
 
     alpha0Plus  = alpha[2*0 + 0]; # coefficient of forward-traveling fundamental mode
     alpha0Minus = alpha[2*0 + 1]; # coefficient of backward-traveling fundamental mode
 
     print("refl:, {}, {:.8f}".format(Lt, abs(alpha0Minus)**2))
-    
+
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()
     parser.add_argument('-Lt',  type=float, default=3.0, help='taper length (default: 3.0)')

python/Makefile.am

@@ -1,4 +1,5 @@
 HPPFILES=                                    \
+ $(top_srcdir)/src/meep.hpp                  \
  $(top_srcdir)/src/meep_internals.hpp        \
  $(top_srcdir)/src/bicgstab.hpp              \
  $(top_srcdir)/src/meep/vec.hpp              \

python/simulation.py

@@ -1024,13 +1024,16 @@ def get_dft_array(self, dft_obj, component, num_freq):
         else:
             raise ValueError("Invalid type of dft object: {}".format(dft_obj))
 
-    def get_eigenmode_coefficients(self, flux, bands, kpoint_func=None):
+    def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY,
+                                   eig_vol=None, eig_resolution=0, eig_tolerance=1e-7, kpoint_func=None):
         if self.fields is None:
             raise ValueError("Fields must be initialized before calling get_eigenmode_coefficients")
+        if eig_vol is None:
+            eig_vol = flux.where
 
         num_bands = len(bands)
         coeffs = np.zeros(2 * num_bands, dtype=np.complex128)
-        self.fields.get_eigenmode_coefficients(flux, np.array(bands, dtype=np.intc), coeffs, None, kpoint_func)
+        self.fields.get_eigenmode_coefficients(flux, eig_vol, np.array(bands, dtype=np.intc), eig_parity, eig_resolution, eig_tolerance, coeffs, None, kpoint_func)
 
         return coeffs
 

src/meep.hpp

@@ -1441,10 +1441,10 @@ class fields {
   // values of field components at arbitrary points in space.
   // call destroy_eigenmode_data() to deallocate it when finished.
   void *get_eigenmode(double omega_src, direction d, const volume where,
-	              const volume eig_vol, int band_num,
-		      const vec &kpoint, bool match_frequency=true,
-                      int parity=0, double resolution=0.0,
-                      double eigensolver_tol=1.0e-7, bool verbose=false);
+                      const volume eig_vol, int band_num,
+                      const vec &kpoint, bool match_frequency,
+                      int parity, double resolution,
+                      double eigensolver_tol, bool verbose=false);
 
   void add_eigenmode_source(component c, const src_time &src,
 	  		    direction d, const volume &where,
@@ -1457,11 +1457,22 @@ class fields {
 			    std::complex<double> A(const vec &) = 0);
 
   void get_eigenmode_coefficients(dft_flux flux,
-                                  int *bands, int num_bands,
+                                  const volume &eig_vol,
+                                  int *bands, int num_bands, int parity,
+                  			          double eig_resolution, double eigensolver_tol,
                                   std::complex<double> *coeffs,
                                   double *vgrp, kpoint_func user_kpoint_func=0,
                                   void *user_kpoint_data=0);
 
+  void get_eigenmode_coefficients(dft_flux flux,
+                                  int *bands, int num_bands, int parity,
+                                  std::complex<double> *coeffs,
+                                  double *vgrp, kpoint_func user_kpoint_func=0,
+                                  void *user_kpoint_data=0) {
+      get_eigenmode_coefficients(flux, flux.where, bands, num_bands, parity,
+          0.0, 1e-7, coeffs, vgrp, user_kpoint_func, user_kpoint_data);
+  }
+
   // initialize.cpp:
   void initialize_field(component, std::complex<double> f(const vec &));
   void initialize_with_nth_te(int n);