Fix adjoint gradient with conductivities (NanoComp#1830)

* damp_fix

* increase run time

Co-authored-by: Mo Chen <mochen@Mos-MacBook-Pro.local>

python/tests/test_adjoint_solver.py

@@ -261,7 +261,86 @@ def J(dft_mon):
     sim.reset_meep()
 
     return f, dJ_du
+
+def forward_simulation_damping(design_params, frequencies=None, mat2=silicon):
+    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
+                              mp.air,
+                              mat2,
+                              weights=design_params.reshape(Nx,Ny),
+                              damping = 3.14*fcen)
+
+    matgrid_geometry = [mp.Block(center=mp.Vector3(),
+                                 size=mp.Vector3(design_region_size.x,design_region_size.y,0),
+                                 material=matgrid)]
+
+    geometry = waveguide_geometry + matgrid_geometry
+
+    sim = mp.Simulation(resolution=resolution,
+                        cell_size=cell_size,
+                        boundary_layers=pml_xy,
+                        sources=wvg_source,
+                        geometry=geometry)
+
+    if not frequencies:
+        frequencies = [fcen]
+
+    mode = sim.add_mode_monitor(frequencies,
+                                    mp.ModeRegion(center=mp.Vector3(0.5*sxy-dpml-0.1),
+                                                  size=mp.Vector3(0,sxy-2*dpml,0)),
+                                    yee_grid=True,
+                                    eig_parity=eig_parity)
+
+    sim.run(until_after_sources=mp.stop_when_dft_decayed())
+
+
+    coeff = sim.get_eigenmode_coefficients(mode,[1],eig_parity).alpha[0,:,0]
+    S12 = np.power(np.abs(coeff),2)
+    sim.reset_meep()
+    return S12
+
+def adjoint_solver_damping(design_params, frequencies=None, mat2=silicon):
+    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
+                              mp.air,
+                              mat2,
+                              weights=np.ones((Nx,Ny)),
+                              damping = 3.14*fcen)
+    matgrid_region = mpa.DesignRegion(matgrid,
+                                      volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_region_size.x, design_region_size.y, 0)))
+
+    matgrid_geometry = [mp.Block(center=matgrid_region.center,
+                                 size=matgrid_region.size,
+                                 material=matgrid)]
 
+    geometry = waveguide_geometry + matgrid_geometry
+
+    sim = mp.Simulation(resolution=resolution,
+                        cell_size=cell_size,
+                        boundary_layers=pml_xy,
+                        sources=wvg_source,
+                        geometry=geometry)
+
+    if not frequencies:
+        frequencies = [fcen]
+
+    obj_list = [mpa.EigenmodeCoefficient(sim, mp.Volume(center=mp.Vector3(0.5*sxy-dpml-0.1),
+                                        size=mp.Vector3(0,sxy-2*dpml,0)), 1, eig_parity=eig_parity)]
+
+    def J(mode_mon):
+        return npa.power(npa.abs(mode_mon),2)
+
+
+    opt = mpa.OptimizationProblem(
+        simulation=sim,
+        objective_functions=J,
+        objective_arguments=obj_list,
+        design_regions=[matgrid_region],
+        frequencies=frequencies,
+        minimum_run_time=150)
+
+    f, dJ_du = opt([design_params])
+
+    sim.reset_meep()
+    return f, dJ_du
 
 def mapping(x,filter_radius,eta,beta):
     filtered_field = mpa.conic_filter(x,
@@ -403,6 +482,32 @@ def test_complex_fields(self):
             print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
             tol = 0.005 if mp.is_single_precision() else 0.0008
             self.assertClose(adj_scale,fd_grad,epsilon=tol)
+
+    def test_damping(self):
+        print("*** TESTING CONDUCTIVITIES ***")
+
+        for frequencies in [[1/1.58, fcen, 1/1.53]]:
+            ## compute gradient using adjoint solver
+            adjsol_obj, adjsol_grad = adjoint_solver_damping(p, frequencies)
+
+            ## compute unperturbed S12
+            S12_unperturbed = forward_simulation_damping(p, frequencies)
+
+            ## 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_damping(p+dp, frequencies)
+
+            ## 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.06 if mp.is_single_precision() else 0.03
+            self.assertClose(adj_scale,fd_grad,epsilon=tol)
 
     def test_offdiagonal(self):
         print("*** TESTING OFFDIAGONAL COMPONENTS ***")

src/meepgeom.cpp

@@ -2572,7 +2572,7 @@ void eff_chi1inv_row_disp(meep::component c, std::complex<double> chi1inv_row[3]
     vector3 dummy;
     dummy.x = dummy.y = dummy.z = 0.0;
     double conductivityCur = vec_to_value(mm->D_conductivity_diag, dummy, i);
-    a = std::complex<double>(1.0, conductivityCur / (freq));
+    a = std::complex<double>(1.0, conductivityCur / (2*meep::pi*freq));
 
     // compute lorentzian component
     b = cvec_to_value(mm->epsilon_diag, mm->epsilon_offdiag, i);
@@ -2613,6 +2613,23 @@ void eff_chi1inv_row_disp(meep::component c, std::complex<double> chi1inv_row[3]
     }
 }
 
+std::complex<double> cond_cmp(meep::component c, const meep::vec &r, double freq, geom_epsilon *geps) {
+  // locate the proper material
+  material_type md;
+  geps->get_material_pt(md, r);
+  const medium_struct *mm = &(md->medium);
+
+  // get the row we care about
+      switch (component_direction(c)) {
+      case meep::X:
+      case meep::R: return std::complex<double>(1.0, mm->D_conductivity_diag.x / (2*meep::pi*freq));
+      case meep::Y:
+      case meep::P: return std::complex<double>(1.0, mm->D_conductivity_diag.y / (2*meep::pi*freq));
+      case meep::Z: return std::complex<double>(1.0, mm->D_conductivity_diag.z / (2*meep::pi*freq));
+      case meep::NO_DIRECTION: meep::abort("Invalid adjoint field component");
+    }
+}
+
 std::complex<double> get_material_gradient(
     const meep::vec &r,               // current point
     const meep::component adjoint_c,  // adjoint field component
@@ -2686,7 +2703,7 @@ std::complex<double> get_material_gradient(
     u[idx] = orig;
 
     for (int i=0;i<3;i++) dA_du[i] = (row_1[i] - row_2[i])/(2*du);
-    return dA_du[dir_idx] * fields_f;
+    return dA_du[dir_idx] * fields_f  * cond_cmp(forward_c,r,freq,geps);
   }  
 }
 
@@ -2912,6 +2929,11 @@ void material_grids_addgradient(double *v, size_t ng, std::complex<double> *fiel
         meep::ivec ip = gv.iloc(adjoint_c,idx);
         meep::vec p   = gv.loc(adjoint_c,idx);
         std::complex<double> adj = GET_FIELDS(fields_a,ci_adjoint,f_i,idx_fields);
+
+        material_type md;
+        geps->get_material_pt(md, p);
+        if (!md->trivial) adj *= cond_cmp(adjoint_c,p,frequencies[f_i], geps);
+
         double cyl_scale;
         int ci_forward = 0;
         FOR_MY_COMPONENTS(forward_c) {