update tests and add comments to code

python/adjoint/filters.py

@@ -7,6 +7,8 @@
 import meep as mp
 from scipy import special
 from scipy import signal
+import skfmm
+from scipy import interpolate
 
 def _proper_pad(x,n):
     '''
@@ -854,3 +856,19 @@ def gray_indicator(x):
     density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218.
     '''
     return npa.mean(4 * x.flatten() * (1 - x.flatten())) * 100
+
+def make_sdf(data):
+    '''
+    Assume the input, data, is the desired output shape
+    (i.e. 1d, 2d, or 3d) and that it's values are between
+    0 and 1.
+    '''
+    # create signed distance function
+    sd = skfmm.distance(data- 0.5 , dx = 1)
+
+    # interpolate zero-levelset onto 0.5-levelset
+    x = [np.min(sd.flatten()), 0, np.max(sd.flatten())]
+    y = [0, 0.5, 1]
+    f = interpolate.interp1d(x, y, kind='linear')
+
+    return f(sd)

python/requirements.txt

@@ -7,3 +7,4 @@ numpy
 parameterized
 pytest
 scipy
+scikit-fmm

python/tests/test_material_grid.py

@@ -13,7 +13,6 @@
 except:
     import adjoint as mpa
 import numpy as np
-from scipy.ndimage import gaussian_filter
 import unittest
 
 rad = 0.301943
@@ -92,8 +91,7 @@ def compute_transmittance(matgrid_symmetry=False):
 
         return tran
 
-
-def compute_resonant_mode_2d(res,default_mat=False):
+def compute_resonant_mode_2d(res,radius=rad,default_mat=False,cylinder=False,cylinder_matgrid=True,do_averaging=True):
         fcen = 0.3
         df = 0.2*fcen
         sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
@@ -110,24 +108,37 @@ def compute_resonant_mode_2d(res,default_mat=False):
         x = np.linspace(-0.5*matgrid_size.x,0.5*matgrid_size.x,Nx)
         y = np.linspace(-0.5*matgrid_size.y,0.5*matgrid_size.y,Ny)
         xv, yv = np.meshgrid(x,y)
-        weights = np.sqrt(np.square(xv) + np.square(yv)) < rad
-        filtered_weights = gaussian_filter(weights,
-                                           sigma=3.0,
-                                           output=np.double)
-
+        weights = mpa.make_sdf(np.sqrt(np.square(xv) + np.square(yv)) < radius)
+        if cylinder:
+            weights = 0.0*weights+1.0          
+
         matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                                   mp.air,
                                   Si,
-                                  weights=filtered_weights,
-                                  do_averaging=True,
-                                  beta=1000,
+                                  weights=weights,
+                                  beta=np.inf,
                                   eta=0.5)
-
-        geometry = [mp.Block(center=mp.Vector3(),
+
+        # Use a cylinder object, not a block object
+        if cylinder:
+            # within the cylinder, use a (uniform) material grid
+            if cylinder_matgrid:
+                geometry = [mp.Cylinder(center=mp.Vector3(),
+                             radius=radius,
+                             material=matgrid)]
+            # within the cylinder, just use a normal medium
+            else:
+                geometry = [mp.Cylinder(center=mp.Vector3(),
+                             radius=radius,
+                             material=Si)]
+        # use a block object
+        else:
+            geometry = [mp.Block(center=mp.Vector3(),
                              size=mp.Vector3(matgrid_size.x,matgrid_size.y,0),
                              material=matgrid)]
 
         sim = mp.Simulation(resolution=res,
+                            eps_averaging=do_averaging,
                             cell_size=cell_size,
                             default_material=matgrid if default_mat else mp.Medium(),
                             geometry=geometry if not default_mat else [],
@@ -136,7 +147,7 @@ def compute_resonant_mode_2d(res,default_mat=False):
 
         h = mp.Harminv(mp.Hz, mp.Vector3(0.3718,-0.2076), fcen, df)
         sim.run(mp.after_sources(h),
-                until_after_sources=200)
+                until_after_sources=300)
 
         try:
             for m in h.modes:
@@ -147,13 +158,12 @@ def compute_resonant_mode_2d(res,default_mat=False):
 
         return freq
 
-
-def compute_resonant_mode_mpb(resolution=1024):
+def compute_resonant_mode_mpb(resolution=512):
     geometry = [mp.Cylinder(rad, material=Si)]
     geometry_lattice = mp.Lattice(cell_size)
 
     ms = mpb.ModeSolver(num_bands=1,
-                        k_points=[k_point],
+                        k_points=[mp.cartesian_to_reciprocal(k_point, geometry_lattice)],
                         geometry=geometry,
                         geometry_lattice=geometry_lattice,
                         resolution=resolution)
@@ -165,23 +175,51 @@ def compute_resonant_mode_mpb(resolution=1024):
 class TestMaterialGrid(unittest.TestCase):
 
     def test_subpixel_smoothing(self):
+        res = 25
+
+        def subpixel_test_matrix(radius,do_averaging):
+            freq_cylinder = compute_resonant_mode_2d(res, radius, default_mat=False, cylinder=True, cylinder_matgrid=False, do_averaging=do_averaging)
+            freq_matgrid = compute_resonant_mode_2d(res, radius, default_mat=False, cylinder=False, do_averaging=do_averaging)
+            freq_matgrid_cylinder = compute_resonant_mode_2d(res, radius, default_mat=False, cylinder=True, cylinder_matgrid=True, do_averaging=do_averaging)
+            return [freq_cylinder,freq_matgrid,freq_matgrid_cylinder]
+
+        # when smoothing is off, all three tests should be identical to machine precision
+        no_smoothing = subpixel_test_matrix(rad,False)
+        self.assertEqual(no_smoothing[0], no_smoothing[1])
+        self.assertEqual(no_smoothing[1], no_smoothing[2])
+
+        # when we slightly perturb the radius, the results should be the same as before.
+        no_smoothing_perturbed = subpixel_test_matrix(rad+0.01/res,False)
+        self.assertEqual(no_smoothing, no_smoothing_perturbed)
+
+        # when smoothing is on, the simple material results should be different from the matgrid results
+        smoothing = subpixel_test_matrix(rad,True)
+        self.assertNotEqual(smoothing[0], smoothing[1])
+        self.assertAlmostEqual(smoothing[1], smoothing[2], 3)
+
+        # when we slighty perturb the radius, the results should all be different from before.
+        smoothing_perturbed = subpixel_test_matrix(rad+0.01/res,True)
+        self.assertNotEqual(smoothing, smoothing_perturbed)
+
         # "exact" frequency computed using MPB
         freq_ref = compute_resonant_mode_mpb()
-
-        res = [25, 50]
+        print(freq_ref)
+        res = [50, 100]
         freq_matgrid = []
         for r in res:
-            freq_matgrid.append(compute_resonant_mode_2d(r))
+            freq_matgrid.append(compute_resonant_mode_2d(r,cylinder=False))
             # verify that the frequency of the resonant mode is
             # approximately equal to the reference value
             self.assertAlmostEqual(freq_ref, freq_matgrid[-1], 2)
+        print("results:    ",freq_ref,freq_matgrid)
 
         # verify that the relative error is decreasing with increasing resolution
         # and is better than linear convergence because of subpixel smoothing
         self.assertLess(abs(freq_matgrid[1]-freq_ref)*(res[1]/res[0]),
                         abs(freq_matgrid[0]-freq_ref))
 
-        freq_matgrid_default_mat = compute_resonant_mode_2d(res[0], True)
+        # ensure that a material grid as a default material works
+        freq_matgrid_default_mat = compute_resonant_mode_2d(res[0], rad, True)
         self.assertAlmostEqual(freq_matgrid[0], freq_matgrid_default_mat)
 
     def test_symmetry(self):

src/meepgeom.cpp

@@ -362,8 +362,14 @@ cvector3 to_geom_object_coords_VJP(cvector3 v, const geometric_object *o) {
     /* case geometric_object::CYLINDER:
        NOT YET IMPLEMENTED */
     case geometric_object::BLOCK: {
-      vector3 v_real = cvector3_re(v);
-      vector3 v_imag = cvector3_im(v);
+      /* In order to leverage the underlying libctl infrastructure
+      *and* the dual library that makes computing derivatives so easy,
+      me perform a bit of a trick: we use the *complex* libctl vector,
+      storing the real and dual parts of the dual library and *manually*
+      propagate the dual through existing libraries as needed.
+      */
+      vector3 v_real = cvector3_re(v); // real part
+      vector3 v_imag = cvector3_im(v); // "dual" part
 
       vector3 size = o->subclass.block_data->size;
       if (size.x != 0.0) {v_real.x /= size.x; v_imag.x /= size.x;}
@@ -379,6 +385,9 @@ cvector3 to_geom_object_coords_VJP(cvector3 v, const geometric_object *o) {
 }
 
 cvector3 material_grid_grad(vector3 p, material_data *md, const geometric_object *o) {
+  /* computes the actual spatial gradient at point `p` 
+  for the specified material grid `md`. */
+
   if (!is_material_grid(md)) {meep::abort("Invalid material grid detected.\n"); }
 
   cvector3 gradient = cvector_zero();
@@ -468,6 +477,14 @@ void map_lattice_coordinates(double &px, double &py, double &pz) {
 }
 
 cvector3 matgrid_grad(vector3 p, geom_box_tree tp, int oi, material_data *md) {
+  /* loops through all the material grids at a current point
+  and computes the final *spatial* gradient (w.r.t. x,y,z)
+  after all appropriate transformations (e.g. due to 
+  overlapping grids). Calls the helper function, `material_grid_grad`,
+  which is what actually computes the spatial gradient.
+  */
+
+  // check for proper overlapping grids
   if (md->material_grid_kinds == material_data::U_MIN ||
       md->material_grid_kinds == material_data::U_PROD)
     meep::abort("%s:%i:matgrid_grad does not support overlapping grids with U_MIN or U_PROD\n",__FILE__,__LINE__);
@@ -493,6 +510,7 @@ cvector3 matgrid_grad(vector3 p, geom_box_tree tp, int oi, material_data *md) {
     ++matgrid_val_count;
   }
 
+  // compensate for overlapping grids
   if (md->material_grid_kinds == material_data::U_MEAN)
     gradient = cvector3_scale(1.0/matgrid_val_count,gradient);
 
@@ -1250,12 +1268,11 @@ duals::duald get_material_grid_fill(meep::ndim dim, duals::duald d, double r, du
     return -1.0; // garbage fill
   } else {
     if (dim == meep::D1)
-      rel_fill = (r-d)/(2*r);
-    else if (dim == meep::D2 || dim == meep::Dcyl){
-      rel_fill = (1/(r*r*meep::pi)) * (r*r*acos(d/r)-d*sqrt(r*r-d*d));
-    }
+      return (r-d)/(2*r);
+    else if (dim == meep::D2 || dim == meep::Dcyl)
+      return (1/(r*r*meep::pi)) * (r*r*acos(d/r)-d*sqrt(r*r-d*d));
     else if (dim == meep::D3)
-      rel_fill = (((r-d)*(r-d))/(4*meep::pi*r*r*r))*(2*r+d);
+      return (((r-d)*(r-d))/(4*meep::pi*r*r*r))*(2*r+d);
   }
 
   return rel_fill;
@@ -1593,20 +1610,7 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3]
   material =
       (material_type)material_of_unshifted_point_in_tree_inobject(p, restricted_tree, &inobject);
   material_data *md = material;
-  meep::vec gradient(zero_vec(v.dim));
-  double uval = 0;
-  // TODO cleanup and remove
-  if (md->which_subclass == material_data::MATERIAL_GRID) {
-    geom_box_tree tp;
-    int oi;
-    tp = geom_tree_search(p, restricted_tree, &oi);
-    //gradient = matgrid_grad(p, tp, oi, md);
-    //uval = matgrid_val(p, tp, oi, md)+this->u_p;
-  }
-  else {
-    gradient = normal_vector(meep::type(c), v);
-  }
-
+  meep::vec gradient = normal_vector(meep::type(c), v);
   get_material_pt(material, v.center());
   material_epsmu(meep::type(c), material, &chi1p1, &chi1p1_inv);
   material_gc(material);
@@ -1634,73 +1638,47 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3]
   integer errflag;
   double meps, minveps;
 
-  if (md->which_subclass == material_data::MATERIAL_GRID) {
-    number xmin[1], xmax[1];
-    matgrid_volavg mgva;
-    mgva.dim = v.dim;
-    mgva.ugrad_abs = meep::abs(gradient);
-    mgva.uval = uval;
-    mgva.rad = v.diameter()/2;
-    mgva.beta = md->beta;
-    mgva.eta = md->eta;
-    mgva.eps1 = (md->medium_1.epsilon_diag.x+md->medium_1.epsilon_diag.y+md->medium_1.epsilon_diag.z)/3;
-    mgva.eps2 = (md->medium_2.epsilon_diag.x+md->medium_2.epsilon_diag.y+md->medium_2.epsilon_diag.z)/3;
-    xmin[0] = -v.diameter()/2;
-    xmax[0] = v.diameter()/2;
-#ifdef CTL_HAS_COMPLEX_INTEGRATION
-    cnumber ret = cadaptive_integration(matgrid_ceps_func, xmin, xmax, 1, (void *)&mgva, 0, tol, maxeval,
-                                        &esterr, &errflag);
-    meps = ret.re;
-    minveps = ret.im;
-#else
-    meps = adaptive_integration(matgrid_eps_func, xmin, xmax, 1, (void *)&mgva, 0, tol, maxeval, &esterr,
-                                &errflag);
-    minveps = adaptive_integration(matgrid_inveps_func, xmin, xmax, 1, (void *)&mgva, 0, tol, maxeval, &esterr,
-                                   &errflag);
-#endif
+  integer n;
+  number xmin[3], xmax[3];
+  vector3 gvmin, gvmax;
+  gvmin = vec_to_vector3(v.get_min_corner());
+  gvmax = vec_to_vector3(v.get_max_corner());
+  xmin[0] = gvmin.x;
+  xmax[0] = gvmax.x;
+  if (dim == meep::Dcyl) {
+    xmin[1] = gvmin.z;
+    xmin[2] = gvmin.y;
+    xmax[1] = gvmax.z;
+    xmax[2] = gvmax.y;
   }
   else {
-    integer n;
-    number xmin[3], xmax[3];
-    vector3 gvmin, gvmax;
-    gvmin = vec_to_vector3(v.get_min_corner());
-    gvmax = vec_to_vector3(v.get_max_corner());
-    xmin[0] = gvmin.x;
-    xmax[0] = gvmax.x;
-    if (dim == meep::Dcyl) {
-      xmin[1] = gvmin.z;
-      xmin[2] = gvmin.y;
-      xmax[1] = gvmax.z;
-      xmax[2] = gvmax.y;
-    }
-    else {
-      xmin[1] = gvmin.y;
-      xmin[2] = gvmin.z;
-      xmax[1] = gvmax.y;
-      xmax[2] = gvmax.z;
-    }
-    if (xmin[2] == xmax[2])
-      n = xmin[1] == xmax[1] ? 1 : 2;
-    else
-      n = 3;
-    double vol = 1;
-    for (int i = 0; i < n; ++i)
-      vol *= xmax[i] - xmin[i];
-    if (dim == meep::Dcyl) vol *= (xmin[0] + xmax[0]) * 0.5;
-    eps_ever_negative = 0;
-    func_ft = meep::type(c);
+    xmin[1] = gvmin.y;
+    xmin[2] = gvmin.z;
+    xmax[1] = gvmax.y;
+    xmax[2] = gvmax.z;
+  }
+  if (xmin[2] == xmax[2])
+    n = xmin[1] == xmax[1] ? 1 : 2;
+  else
+    n = 3;
+  double vol = 1;
+  for (int i = 0; i < n; ++i)
+    vol *= xmax[i] - xmin[i];
+  if (dim == meep::Dcyl) vol *= (xmin[0] + xmax[0]) * 0.5;
+  eps_ever_negative = 0;
+  func_ft = meep::type(c);
 #ifdef CTL_HAS_COMPLEX_INTEGRATION
-    cnumber ret = cadaptive_integration(ceps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval,
-                                        &esterr, &errflag);
-    meps = ret.re / vol;
-    minveps = ret.im / vol;
+  cnumber ret = cadaptive_integration(ceps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval,
+                                      &esterr, &errflag);
+  meps = ret.re / vol;
+  minveps = ret.im / vol;
 #else
-    meps = adaptive_integration(eps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr,
-                                &errflag) / vol;
-    minveps = adaptive_integration(inveps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr,
-                                   &errflag) / vol;
+  meps = adaptive_integration(eps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr,
+                              &errflag) / vol;
+  minveps = adaptive_integration(inveps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr,
+                                  &errflag) / vol;
 #endif
-  }
+
   if (eps_ever_negative) // averaging negative eps causes instability
     minveps = 1.0 / (meps = eps(v.center()));
   {