Fix factor of 2 in adjoint gradients when not using complex fields

python/adjoint/objective.py

@@ -92,6 +92,8 @@ def _adj_src_scale(self, include_resolution=True):
         else:
             # multi frequency simulations
             scale = dV * iomega / adj_src_phase
+        if not self.sim.force_complex_fields:
+            scale *= 2
         return scale
 
     def _create_time_profile(self, fwidth_frac=0.1):
@@ -262,7 +264,7 @@ def place_adjoint_source(self, dJ):
                 for yi in range(y_dim):
                     for xi in range(x_dim):
                         '''We only need to add a current source if the
-                        jacobian is nonzero for all frequencies at 
+                        jacobian is nonzero for all frequencies at
                         that particular point. Otherwise, the fitting
                         algorithm is going to fail.
                         '''
@@ -340,7 +342,7 @@ def place_adjoint_source(self, dJ):
                     time_src.frequency,
                     self._frequencies,
                     scale,
-                    dt,
+                    self.sim.fields.dt,
                 )
                 (num_basis, num_pts) = src.nodes.shape
                 for basis_i in range(num_basis):

python/tests/test_adjoint_solver.py

@@ -53,13 +53,13 @@
                               eig_parity=eig_parity)]
 
 
-def forward_simulation(design_params,mon_type,frequencies=None):
+def forward_simulation(design_params,mon_type,frequencies=None, use_complex=False):
     matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                               mp.air,
                               silicon,
                               weights=design_params.reshape(Nx,Ny),
                               grid_type='U_MEAN')
-            
+
     matgrid_geometry = [mp.Block(center=mp.Vector3(),
                                  size=mp.Vector3(design_shape.x,design_shape.y,0),
                                  material=matgrid)]
@@ -70,7 +70,8 @@ def forward_simulation(design_params,mon_type,frequencies=None):
                         cell_size=cell_size,
                         boundary_layers=boundary_layers,
                         sources=sources,
-                        geometry=geometry)
+                        geometry=geometry,
+                        force_complex_fields=use_complex)
     if not frequencies:
         frequencies = [fcen]
 
@@ -107,7 +108,7 @@ 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, use_complex=False):
     matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                               mp.air,
                               silicon,
@@ -127,7 +128,8 @@ def adjoint_solver(design_params, mon_type, frequencies=None):
                         cell_size=cell_size,
                         boundary_layers=boundary_layers,
                         sources=sources,
-                        geometry=geometry)
+                        geometry=geometry,
+                        force_complex_fields=use_complex)
     if not frequencies:
         frequencies = [fcen]
 
@@ -182,14 +184,14 @@ def test_adjoint_solver_DFT_fields(self):
 
             ## compute unperturbed S12
             S12_unperturbed = forward_simulation(p, MonitorObject.DFT, frequencies)
-            
+
             ## compare objective results
             print("|Ez|^2 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
             self.assertVectorsClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)
 
             ## compute perturbed S12
             S12_perturbed = forward_simulation(p+dp, MonitorObject.DFT, frequencies)
-            
+
             ## compare gradients
             if adjsol_grad.ndim < 2:
                 adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
@@ -203,28 +205,29 @@ def test_adjoint_solver_DFT_fields(self):
     def test_adjoint_solver_eigenmode(self):
         print("*** TESTING EIGENMODE ADJOINT FEATURES ***")
         ## test the single frequency and multi frequency cases
-        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)
-
-            ## compute unperturbed S12
-            S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies)
-
-            ## compare objective results
-            print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
-            self.assertVectorsClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)
-
-            ## compute perturbed S12
-            S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, 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.04 if mp.is_single_precision() else 0.03
-            self.assertVectorsClose(adj_scale,fd_grad,epsilon=tol)
+        for use_complex in [True, False]:
+            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, use_complex)
+
+                ## compute unperturbed S12
+                S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies, use_complex)
+
+                ## compare objective results
+                print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
+                self.assertVectorsClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)
+
+                ## compute perturbed S12
+                S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies, use_complex)
+
+                ## 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.1 if mp.is_single_precision() else 0.03
+                self.assertVectorsClose(adj_scale,fd_grad,epsilon=tol)
 
 
     def test_gradient_backpropagation(self):

src/meepgeom.cpp

@@ -408,7 +408,7 @@ meep::vec material_grid_grad(vector3 p, material_data *md, const geometric_objec
 #undef D
 
   // [du_dx,du_dy,du_dz] is the gradient ∇u with respect to the transformed coordinate
-  // r1 of the matgrid_val function but what we want is the gradient of u(g(r2)) with 
+  // r1 of the matgrid_val function but what we want is the gradient of u(g(r2)) with
   // respect to r2 where g(r2) is the to_geom_object_coords function (in libctl/utils/geom.c).
   // computing this quantity involves using the chain rule and thus the vector-Jacobian product
   // ∇u J where J is the Jacobian matrix of g.
@@ -2579,9 +2579,9 @@ double get_material_gradient(
   if ((mm->E_susceptibilities.size() == 0) && mm->D_conductivity_diag.x == 0 &&
       mm->D_conductivity_diag.y == 0 && mm->D_conductivity_diag.z == 0){
         switch (field_dir){
-          case meep::Ex: return 2 * (m2->epsilon_diag.x - m1->epsilon_diag.x) * (fields_a * fields_f).real();
-          case meep::Ey: return 2 * (m2->epsilon_diag.y - m1->epsilon_diag.y) * (fields_a * fields_f).real();
-          case meep::Ez: return 2 * (m2->epsilon_diag.z - m1->epsilon_diag.z) * (fields_a * fields_f).real();
+          case meep::Ex: return (m2->epsilon_diag.x - m1->epsilon_diag.x) * (fields_a * fields_f).real();
+          case meep::Ey: return (m2->epsilon_diag.y - m1->epsilon_diag.y) * (fields_a * fields_f).real();
+          case meep::Ez: return (m2->epsilon_diag.z - m1->epsilon_diag.z) * (fields_a * fields_f).real();
           default: meep::abort("Invalid field component specified in gradient.");
         }
       }
@@ -2612,7 +2612,7 @@ double get_material_gradient(
     meep::abort("Invalid adjoint field component");
 
   std::complex<double> result = fields_a * dA_du[3 * dir_idx + dir_idx] * fields_f;
-  return 2 * result.real();
+  return result.real();
 }
 
 /* add the weights from linear_interpolate (see the linear_interpolate