Approximate material grid smoothing (β=∞)

libpympb/pympb.cpp

@@ -376,12 +376,9 @@ int mode_solver::mean_epsilon(symmetric_matrix *meps, symmetric_matrix *meps_inv
 
   bool o1_is_var = o1 && meep_geom::is_variable(o1->material);
   bool o2_is_var = o2 && meep_geom::is_variable(o2->material);
-  bool default_is_var_or_file =
-      meep_geom::is_variable(default_material) || meep_geom::is_file(default_material);
+  bool default_is_var = meep_geom::is_variable(default_material);
 
-  if (o1_is_var || o2_is_var ||
-      (default_is_var_or_file &&
-       (!o1 || !o2 || meep_geom::is_file(o1->material) || meep_geom::is_file(o2->material)))) {
+  if (o1_is_var || o2_is_var || (default_is_var && (!o1 || !o2))) {
     return 0; /* arbitrary material functions are non-analyzable */
   }
 

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/adjoint/utils.py

@@ -19,7 +19,7 @@
 _Z_AXIS = 1
 
 # default finite difference step size when calculating Aᵤ
-FD_DEFAULT = 1e-3
+FD_DEFAULT = 1e-6
 
 class DesignRegion(object):
     def __init__(

python/meep.i

@@ -1504,8 +1504,6 @@ void _get_gradient(PyObject *grad, double scalegrad,
 %ignore material_type_equal;
 %ignore is_variable;
 %ignore is_variable;
-%ignore is_file;
-%ignore is_file;
 %ignore is_medium;
 %ignore is_medium;
 %ignore is_metal;
@@ -1986,7 +1984,7 @@ meep_geom::geom_epsilon* _set_materials(meep::structure * s,
         geps = existing_geps;
     } else {
         geps = meep_geom::make_geom_epsilon(s, &gobj_list, center, _ensure_periodicity, _default_material,
-                                                extra_materials);
+                                                extra_materials, use_anisotropic_averaging);
     }
     if (set_materials) {
         meep_geom::set_materials_from_geom_epsilon(s, geps, use_anisotropic_averaging, tol,

python/requirements.txt

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

python/tests/test_material_grid.py

@@ -1,12 +1,24 @@
+'''
+Checks that a material grid works with symmetries, and that the
+new subpixel smoothing feature is accurate.
+
+TODO:
+    Create a 3D test that works well with the new smoothing
+'''
+
 import meep as mp
+from meep import mpb
 try:
     import meep.adjoint as mpa
 except:
     import adjoint as mpa
 import numpy as np
-from scipy.ndimage import gaussian_filter
 import unittest
 
+rad = 0.301943
+k_point = mp.Vector3(0.3892,0.1597,0)
+Si = mp.Medium(index=3.5)
+cell_size = mp.Vector3(1,1,0)
 
 def compute_transmittance(matgrid_symmetry=False):
         resolution = 25
@@ -24,13 +36,16 @@ def compute_transmittance(matgrid_symmetry=False):
         rng = np.random.RandomState(2069588)
 
         w = rng.rand(Nx,Ny)
+        # for subpixel smoothing, the underlying design grid must be smoothly varying
+        w = mpa.tanh_projection(rng.rand(Nx,Ny),1e5,0.5)
+        w = mpa.conic_filter(w,0.25,matgrid_size.x,matgrid_size.y,matgrid_resolution)
         weights = 0.5 * (w + np.fliplr(w)) if not matgrid_symmetry else w
 
         matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                                   mp.air,
                                   mp.Medium(index=3.5),
                                   weights=weights,
-                                  do_averaging=False,
+                                  beta=np.inf,
                                   grid_type='U_MEAN')
 
         geometry = [mp.Block(center=mp.Vector3(),
@@ -76,20 +91,13 @@ def compute_transmittance(matgrid_symmetry=False):
 
         return tran
 
-
-def compute_resonant_mode_2d(res,default_mat=False):
-        cell_size = mp.Vector3(1,1,0)
-
-        rad = 0.301943
-
+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),
                              component=mp.Hz,
                              center=mp.Vector3(-0.1057,0.2094,0))]
 
-        k_point = mp.Vector3(0.3892,0.1597,0)
-
         matgrid_size = mp.Vector3(1,1,0)
         matgrid_resolution = 1200
 
@@ -100,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,
-                                  mp.Medium(index=3.5),
-                                  weights=filtered_weights,
-                                  do_averaging=True,
-                                  beta=1000,
+                                  Si,
+                                  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 [],
@@ -126,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:
@@ -137,109 +158,70 @@ def compute_resonant_mode_2d(res,default_mat=False):
 
         return freq
 
+def compute_resonant_mode_mpb(resolution=512):
+    geometry = [mp.Cylinder(rad, material=Si)]
+    geometry_lattice = mp.Lattice(cell_size)
 
-def compute_resonant_mode_3d(use_matgrid=True):
-    resolution = 25
-
-    wvl = 1.27
-    fcen = 1/wvl
-    df = 0.02*fcen
-
-    nSi = 3.45
-    Si = mp.Medium(index=nSi)
-    nSiO2 = 1.45
-    SiO2 = mp.Medium(index=nSiO2)
-
-    s = 1.0
-    cell_size = mp.Vector3(s,s,s)
-
-    rad = 0.34  # radius of sphere
-
-    if use_matgrid:
-        matgrid_resolution = 2*resolution
-        N = int(s*matgrid_resolution)
-        coord = np.linspace(-0.5*s,0.5*s,N)
-        xv, yv, zv = np.meshgrid(coord,coord,coord)
-
-        weights = np.sqrt(np.square(xv) + np.square(yv) + np.square(zv)) < rad
-        filtered_weights = gaussian_filter(weights,
-                                           sigma=4/resolution,
-                                           output=np.double)
-
-        matgrid = mp.MaterialGrid(mp.Vector3(N,N,N),
-                                  SiO2,
-                                  Si,
-                                  weights=filtered_weights,
-                                  do_averaging=True,
-                                  beta=1000,
-                                  eta=0.5)
-
-        geometry = [mp.Block(center=mp.Vector3(),
-                             size=cell_size,
-                             material=matgrid)]
-    else:
-        geometry = [mp.Sphere(center=mp.Vector3(),
-                              radius=rad,
-                              material=Si)]
-
-    sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=df),
-                         size=mp.Vector3(),
-                         center=mp.Vector3(0.13,0.25,0.06),
-                         component=mp.Ez)]
+    ms = mpb.ModeSolver(num_bands=1,
+                        k_points=[mp.cartesian_to_reciprocal(k_point, geometry_lattice)],
+                        geometry=geometry,
+                        geometry_lattice=geometry_lattice,
+                        resolution=resolution)
 
-    k_point = mp.Vector3(0.23,-0.17,0.35)
-
-    sim = mp.Simulation(resolution=resolution,
-                        cell_size=cell_size,
-                        sources=sources,
-                        default_material=SiO2,
-                        k_point=k_point,
-                        geometry=geometry)
-
-    h = mp.Harminv(mp.Ez, mp.Vector3(-0.2684,0.1185,0.0187), fcen, df)
-
-    sim.run(mp.after_sources(h),
-            until_after_sources=200)
-
-    try:
-        for m in h.modes:
-            print("harminv:, {}, {}, {}".format(resolution,m.freq,m.Q))
-        freq = h.modes[0].freq
-    except:
-        raise RuntimeError("No resonant modes found.")
-
-    return freq
 
+    ms.run_te()
+    return ms.freqs[0]
 
 class TestMaterialGrid(unittest.TestCase):
 
     def test_subpixel_smoothing(self):
-        # "exact" frequency computed using MaterialGrid at resolution = 300
-        freq_ref = 0.29826813873225283
-
-        res = [25, 50]
+        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()
+        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_matgrid_3d(self):
-        freq_matgrid = compute_resonant_mode_3d(True)
-        freq_geomobj = compute_resonant_mode_3d(False)
-        self.assertAlmostEqual(freq_matgrid, freq_geomobj, places=2)
-
-
     def test_symmetry(self):
             tran_nosym = compute_transmittance(False)
             tran_sym = compute_transmittance(True)

python/typemap_utils.cpp

@@ -529,8 +529,12 @@ static int pymaterial_grid_to_material_grid(PyObject *po, material_data *md) {
   if (!PyArray_ISCARRAY(pao)) {
     meep::abort("Numpy array weights must be C-style contiguous.");
   }
-  md->weights = new double[PyArray_SIZE(pao)];
-  memcpy(md->weights, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
+  md->weights.clear();
+  for (size_t i=0;i<PyArray_SIZE(pao);i++){
+    md->weights.push_back(((double *)PyArray_DATA(pao))[i]);
+  }
+  //md->weights = new double[PyArray_SIZE(pao)];
+  //memcpy(md->weights, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
 
   // if needed, combine sus structs to main object
   PyObject *py_e_sus_m1 = PyObject_GetAttrString(po_medium1, "E_susceptibilities");

src/material_data.cpp

@@ -76,10 +76,10 @@ void material_data::copy_from(const material_data &from) {
   }
 
   grid_size = from.grid_size;
-  if (from.weights) {
+  if (from.weights.size() != 0) {
     size_t N = from.grid_size.x * from.grid_size.y * from.grid_size.z;
-    weights = new double[N];
-    std::copy_n(from.weights, N, weights);
+    weights.clear();
+    weights.insert(weights.end(), from.weights.begin(), from.weights.end());
   }
 
   medium_1 = from.medium_1;