diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index ecbf2cedd..3d6c315c1 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -92,6 +92,10 @@ def _adj_src_scale(self, include_resolution=True): else: # multi frequency simulations scale = dV * iomega / adj_src_phase + # compensate for the fact that real fields take the real part of the current, + # which halves the Fourier amplitude at the positive frequency (Re[J] = (J + J*)/2) + if not self.sim.force_complex_fields: + scale *= 2 return scale def _create_time_profile(self, fwidth_frac=0.1): @@ -262,7 +266,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 +344,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): diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index e1136776b..227aab136 100644 --- a/python/tests/test_adjoint_solver.py +++ b/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): diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index f15468ebf..e503e8137 100644 --- a/src/meepgeom.cpp +++ b/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 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