Anisotropic material gradient support

python/adjoint/optimization_problem.py

@@ -13,35 +13,41 @@ def __init__(
             design_parameters,
             volume=None,
             size=None,
-            center=mp.Vector3(),
-            MaterialGrid=None,
+            center=mp.Vector3()
     ):
         self.volume = volume if volume else mp.Volume(center=center, size=size)
         self.size = self.volume.size
         self.center = self.volume.center
         self.design_parameters = design_parameters
         self.num_design_params = design_parameters.num_params
-        self.MaterialGrid = MaterialGrid
 
     def update_design_parameters(self, design_parameters):
         self.design_parameters.update_weights(design_parameters)
+
+    def update_beta(self,beta):
+        self.design_parameters.beta=beta
 
     def get_gradient(self, sim, fields_a, fields_f, frequencies):
+        num_freqs = np.array(frequencies).size
+        shapes = []
         for c in range(3):
+            if (sim.is_cylindrical or sim.dimensions == mp.CYLINDRICAL):
+                # roll r, z, and phi
+                fields_a[c] = np.transpose(fields_a[c],(0,2,3,1))
+                fields_f[c] = np.transpose(fields_f[c],(0,2,3,1))
+            shapes.append(fields_a[c].shape)
             fields_a[c] = fields_a[c].flatten(order='C')
             fields_f[c] = fields_f[c].flatten(order='C')
+        shapes = np.asarray(shapes).flatten(order='C')
         fields_a = np.concatenate(fields_a)
         fields_f = np.concatenate(fields_f)
-        num_freqs = np.array(frequencies).size
-
+
         grad = np.zeros((num_freqs, self.num_design_params))  # preallocate
-
         geom_list = sim.geometry
         f = sim.fields
         vol = sim._fit_volume_to_simulation(self.volume)
         # compute the gradient
-        sim_is_cylindrical = (sim.dimensions == mp.CYLINDRICAL) or sim.is_cylindrical
-        mp._get_gradient(grad,fields_a,fields_f,vol,np.array(frequencies),geom_list,f, sim_is_cylindrical)
+        mp._get_gradient(grad,fields_a,fields_f,sim.gv,vol.swigobj,np.array(frequencies),sim.geps,shapes)
 
         return np.squeeze(grad).T
 
@@ -76,9 +82,9 @@ def __init__(
     ):
 
         self.sim = simulation
-        self.components = [mp.Ex,mp.Ey,mp.Ez]
+        self.components = [mp.Dx,mp.Dy,mp.Dz]
         if self.sim.is_cylindrical or self.sim.dimensions == mp.CYLINDRICAL:
-            self.components = [mp.Er,mp.Ep,mp.Ez]
+            self.components = [mp.Dr,mp.Dp,mp.Dz]
 
         if isinstance(objective_functions, list):
             self.objective_functions = objective_functions
@@ -327,6 +333,7 @@ def adjoint_run(self):
                 for ic, c in enumerate(self.components):
                     for f in range(self.nf):
                         if (self.sim.is_cylindrical or self.sim.dimensions == mp.CYLINDRICAL):
+                            # FIXME probably shouldn't go here...
                             # Addtional factor of 2 for cyldrical coordinate
                             self.a_E[ar][nb][ic][f, :, :, :] = 2 * atleast_3d(self.sim.get_dft_array(dgm, c, f))
                         else:
@@ -486,7 +493,7 @@ def calculate_fd_gradient(
 
         return fd_gradient, fd_gradient_idx
 
-    def update_design(self, rho_vector):
+    def update_design(self, rho_vector, beta=None):
         """Update the design permittivity function.
 
         rho_vector ....... a list of numpy arrays that maps to each design region
@@ -497,6 +504,8 @@ def update_design(self, rho_vector):
                     "Each vector of design variables must contain only one dimension."
                 )
             b.update_design_parameters(rho_vector[bi])
+            if beta:
+                b.update_beta(beta)
 
         self.sim.reset_meep()
         self.current_state = "INIT"

python/meep.i

@@ -844,7 +844,9 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
 //--------------------------------------------------
 
 %inline %{
-void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObject *grid_volume, PyObject *frequencies, PyObject *py_geom_list, PyObject *f, bool sim_is_cylindrical) {
+void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, 
+    meep::grid_volume *grid_volume, meep::volume *where, PyObject *frequencies, 
+    meep_geom::geom_epsilon *geps, PyObject *fields_shapes) {
     // clean the gradient array
     PyArrayObject *pao_grad = (PyArrayObject *)grad;
     if (!PyArray_Check(pao_grad)) meep::abort("grad parameter must be numpy array.");
@@ -867,16 +869,15 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
     if (PyArray_NDIM(pao_fields_f) !=1) {meep::abort("Numpy forward fields array must have 1 dimension.");}
     std::complex<double> *fields_f_c = (std::complex<double> *)PyArray_DATA(pao_fields_f);
 
+    // clean shapes array
+    PyArrayObject *pao_fields_shapes = (PyArrayObject *)fields_shapes;
+    if (!PyArray_Check(pao_fields_shapes)) meep::abort("fields shape parameter must be numpy array.");
+    if (!PyArray_ISCARRAY(pao_fields_shapes)) meep::abort("Numpy fields shape array must be C-style contiguous.");
+    size_t *fields_shapes_c = (size_t *)PyArray_DATA(pao_fields_shapes);
+
     // scalegrad not currently used
     double scalegrad = 1.0;
 
-    // clean the volume object
-    void* where;
-
-    PyObject *swigobj = PyObject_GetAttrString(grid_volume, "swigobj");
-    SWIG_ConvertPtr(swigobj,&where,0,NULL);
-    const meep::volume* where_vol = (const meep::volume*)where;
-
     // clean the frequencies array
     PyArrayObject *pao_freqs = (PyArrayObject *)frequencies;
     if (!PyArray_Check(pao_freqs)) meep::abort("frequencies parameter must be numpy array.");
@@ -885,24 +886,8 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
     npy_intp nf = PyArray_DIMS(pao_freqs)[0];
     if (PyArray_DIMS(pao_grad)[0] != nf) meep::abort("Numpy grad array is allocated for %td frequencies; it should be allocated for %td.",PyArray_DIMS(pao_grad)[0],nf);
 
-    // prepare a geometry_tree
-    // TODO eventually it would be nice to store the geom tree within the structure object so we don't have to recreate it here
-    geometric_object_list *l;
-    l = new geometric_object_list();
-    if (!py_list_to_gobj_list(py_geom_list,l)) meep::abort("Unable to convert geometry tree.");
-    geom_box_tree geometry_tree = calculate_tree(*where_vol,*l);
-
-    // clean the fields pointer
-    void* f_v;
-    SWIG_ConvertPtr(f,&f_v,NULL,NULL);
-    meep::fields* f_c = (meep::fields *)f_v;
-
     // calculate the gradient
-    meep_geom::material_grids_addgradient(grad_c,ng,fields_a_c,fields_f_c,frequencies_c,nf,scalegrad,*where_vol,geometry_tree,f_c, sim_is_cylindrical);
-
-    destroy_geom_box_tree(geometry_tree);
-    delete l;
-    Py_DECREF(swigobj);
+    meep_geom::material_grids_addgradient(grad_c,ng,fields_a_c,fields_f_c,fields_shapes_c,frequencies_c,scalegrad,*grid_volume,*where,geps);
 }
 %}
 
@@ -1413,11 +1398,6 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
 // For some reason SWIG needs the namespaced version too
 %apply material_type_list { meep_geom::material_type_list };
 
-// Typemaps for geom_epsilon object
-%typemap(out) geom_epsilon {
-    $result = PyCapsule_New($1);
-}
-
 // Typemap suite for kpoint_func
 
 %typecheck(SWIG_TYPECHECK_POINTER) (meep::kpoint_func user_kpoint_func, void *user_kpoint_data) {
@@ -2014,12 +1994,11 @@ meep_geom::geom_epsilon* _set_materials(meep::structure * s,
                     const meep::binary_partition *my_bp) {
 
     meep_geom::geom_epsilon *geps;
-    if (existing_geps) {
-        geps = existing_geps;
-    } else {
-        geps = meep_geom::make_geom_epsilon(s, &gobj_list, center, _ensure_periodicity, _default_material,
+    // FIXME we need to properly garbage collect the geps
+    //if (existing_geps)
+    //    delete existing_geps;
+    geps = meep_geom::make_geom_epsilon(s, &gobj_list, center, _ensure_periodicity, _default_material,
                                                 extra_materials);
-    }
     if (set_materials) {
         meep_geom::set_materials_from_geom_epsilon(s, geps, center, use_anisotropic_averaging, tol,
                                              maxeval,alist);

python/simulation.py

@@ -1725,7 +1725,7 @@ def _init_structure(self, k=False):
             self.extra_materials,
             self.split_chunks_evenly,
             True,
-            None,
+            self.geps,
             False,
             None
         )
@@ -1903,7 +1903,7 @@ def set_materials(self, geometry=None, default_material=None):
             self.extra_materials,
             self.split_chunks_evenly,
             True,
-            None,
+            self.geps,
             False,
             None
         )

python/tests/test_adjoint_solver.py

@@ -15,6 +15,7 @@
 resolution = 30
 
 silicon = mp.Medium(epsilon=12)
+sapphire = mp.Medium(epsilon_diag=(10.225,10.225,9.95),epsilon_offdiag=(-0.825,-0.55*np.sqrt(3/2),0.55*np.sqrt(3/2)))
 
 sxy = 5.0
 cell_size = mp.Vector3(sxy,sxy,0)
@@ -60,10 +61,10 @@
 k_point = mp.Vector3(0.23,-0.38)
 
 
-def forward_simulation(design_params, mon_type, frequencies=None):
+def forward_simulation(design_params, mon_type, frequencies=None, mat2=silicon):
     matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                               mp.air,
-                              silicon,
+                              mat2,
                               weights=design_params.reshape(Nx,Ny))
 
     matgrid_geometry = [mp.Block(center=mp.Vector3(),
@@ -114,10 +115,10 @@ def forward_simulation(design_params, mon_type, frequencies=None):
         return Ez2
 
 
-def adjoint_solver(design_params, mon_type, frequencies=None):
+def adjoint_solver(design_params, mon_type, frequencies=None, mat2=silicon):
     matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                               mp.air,
-                              silicon,
+                              mat2,
                               weights=np.ones((Nx,Ny)))
 
     matgrid_region = mpa.DesignRegion(matgrid,
@@ -399,6 +400,31 @@ def test_complex_fields(self):
             tol = 0.005 if mp.is_single_precision() else 0.0008
             self.assertClose(adj_scale,fd_grad,epsilon=tol)
 
+    def test_offdiagonal(self):
+        print("*** TESTING OFFDIAGONAL COMPONENTS ***")
+
+        ## test the single frequency and multi frequency case
+        for frequencies in [[fcen], [1/1.58, fcen, 1/1.53]]:
+            ## compute gradient using adjoint solver
+            adjsol_obj, adjsol_grad = adjoint_solver(p, MonitorObject.EIGENMODE, frequencies, sapphire)
 
+            ## compute unperturbed S12
+            S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies, sapphire)
+
+            ## compare objective results
+            print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
+            self.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-6)
+
+            ## compute perturbed S12
+            S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies, sapphire)
+
+            ## compare gradients
+            if adjsol_grad.ndim < 2:
+                adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
+            adj_scale = (dp[None,:]@adjsol_grad).flatten()
+            fd_grad = S12_perturbed-S12_unperturbed
+            print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
+            tol = 0.04
+            self.assertClose(adj_scale,fd_grad,epsilon=tol)
 if __name__ == '__main__':
     unittest.main()

src/loop_in_chunks.cpp

@@ -229,7 +229,7 @@ void chunkloop_field_components::update_values(ptrdiff_t idx) {
    component of equal_shift (which should be either -2, 0, or +2).
    (equal_shift is there to prevent us from counting edge points twice.) */
 
-static ivec vec2diel_floor(const vec &pt, double a, const ivec &equal_shift) {
+ivec fields::vec2diel_floor(const vec &pt, double a, const ivec &equal_shift) {
   ivec ipt(pt.dim);
   LOOP_OVER_DIRECTIONS(pt.dim, d) {
     ipt.set_direction(d, 1 + 2 * int(floor(pt.in_direction(d) * a - .5)));
@@ -238,7 +238,7 @@ static ivec vec2diel_floor(const vec &pt, double a, const ivec &equal_shift) {
   }
   return ipt;
 }
-static ivec vec2diel_ceil(const vec &pt, double a, const ivec &equal_shift) {
+ivec fields::vec2diel_ceil(const vec &pt, double a, const ivec &equal_shift) {
   ivec ipt(pt.dim);
   LOOP_OVER_DIRECTIONS(pt.dim, d) {
     ipt.set_direction(d, 1 + 2 * int(ceil(pt.in_direction(d) * a - .5)));

src/meep.hpp

@@ -1931,6 +1931,8 @@ class fields {
   int is_phasing();
 
   // loop_in_chunks.cpp
+  static ivec vec2diel_floor(const vec &pt, double a, const ivec &equal_shift);
+  static ivec vec2diel_ceil(const vec &pt, double a, const ivec &equal_shift);
   void loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data, const volume &where,
                       component cgrid = Centered, bool use_symmetry = true,
                       bool snap_unit_dims = false);

src/meep/vec.hpp

@@ -789,6 +789,12 @@ class ivec {
     return result;
   };
 
+  ivec operator+(int s) const {
+    ivec result = *this;
+    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] += s;
+    return result;
+  };
+
   ivec operator+=(const ivec &a) {
     LOOP_OVER_DIRECTIONS(dim, d) t[d] += a.t[d];
     return *this;
@@ -847,6 +853,19 @@ class ivec {
     return result;
   };
 
+  // element-wise product
+  ivec operator*(const ivec &a) const {
+    ivec result = *this;
+    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] *= a.t[d];
+    return result;
+  };
+
+  ivec operator/(int s) const {
+    ivec result = *this;
+    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] /= s;
+    return result;
+  };
+
   vec operator*(double s) const {
     vec result(dim);
     LOOP_OVER_DIRECTIONS(dim, d) result.set_direction(d, t[d] * s);