diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 39fdbfab4..53a97d123 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -204,7 +204,7 @@ def forward_simulation_complex_fields(design_params, frequencies=None): material=matgrid)] sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, + cell_size=cell_size,default_material=silicon, k_point=k_point, boundary_layers=pml_x, sources=pt_source, @@ -213,7 +213,7 @@ def forward_simulation_complex_fields(design_params, frequencies=None): if not frequencies: frequencies = [fcen] - mode = sim.add_dft_fields([mp.Ez], + mode = sim.add_dft_fields([mp.Dz], frequencies, center=mp.Vector3(0.9), size=mp.Vector3(0.2,0.5), @@ -221,15 +221,15 @@ def forward_simulation_complex_fields(design_params, frequencies=None): sim.run(until_after_sources=mp.stop_when_dft_decayed()) - Ez2 = [] + Dz2 = [] for f in range(len(frequencies)): - Ez_dft = sim.get_dft_array(mode, mp.Ez, f) - Ez2.append(np.power(np.abs(Ez_dft[3,9]),2)) - Ez2 = np.array(Ez2) + Dz_dft = sim.get_dft_array(mode, mp.Dz, f) + Dz2.append(np.power(np.abs(Dz_dft[3,9]),2)) + Dz2 = np.array(Dz2) sim.reset_meep() - return Ez2 + return Dz2 def adjoint_solver_complex_fields(design_params, frequencies=None): @@ -249,7 +249,7 @@ def adjoint_solver_complex_fields(design_params, frequencies=None): material=matgrid)] sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, + cell_size=cell_size,default_material=silicon, k_point=k_point, boundary_layers=pml_x, sources=pt_source, @@ -261,7 +261,7 @@ def adjoint_solver_complex_fields(design_params, frequencies=None): obj_list = [mpa.FourierFields(sim, mp.Volume(center=mp.Vector3(0.9), size=mp.Vector3(0.2,0.5)), - mp.Ez)] + mp.Dz)] def J(dft_mon): return npa.power(npa.abs(dft_mon[:,3,9]),2) @@ -505,23 +505,24 @@ def test_complex_fields(self): ## compute gradient using adjoint solver adjsol_obj, adjsol_grad = adjoint_solver_complex_fields(p, frequencies) - ## compute unperturbed |Ez|^2 - Ez2_unperturbed = forward_simulation_complex_fields(p, frequencies) + ## compute unperturbed |Dz|^2 + Dz2_unperturbed = forward_simulation_complex_fields(p, frequencies) ## compare objective results - print("Ez2 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,Ez2_unperturbed)) - self.assertClose(adjsol_obj,Ez2_unperturbed,epsilon=1e-6) + print("Dz2 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,Dz2_unperturbed)) + tol = 1e-5 if mp.is_single_precision() else 1e-6 + self.assertClose(adjsol_obj,Dz2_unperturbed,epsilon=tol) - ## compute perturbed |Ez|^2 - Ez2_perturbed = forward_simulation_complex_fields(p+dp, frequencies) + ## compute perturbed |Dz|^2 + Dz2_perturbed = forward_simulation_complex_fields(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 = Ez2_perturbed-Ez2_unperturbed + fd_grad = Dz2_perturbed-Dz2_unperturbed print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad)) - tol = 0.018 if mp.is_single_precision() else 0.002 + tol = 0.06 if mp.is_single_precision() else 0.002 self.assertClose(adj_scale,fd_grad,epsilon=tol) def test_damping(self): diff --git a/src/dft.cpp b/src/dft.cpp index bc18c2fe2..d7780e090 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -1468,7 +1468,10 @@ std::vector dft_fields::fourier_sourcedata(const volume &wher std::vector idx_arr; std::vector > amp_arr; std::complex EH0 = std::complex(0,0); - sourcedata temp_struct = {component(f->c), idx_arr, f->fc->chunk_idx, amp_arr}; + component c = component(f->c); + direction cd = component_direction(c); + sourcedata temp_struct = {c, idx_arr, f->fc->chunk_idx, amp_arr}; + int position_array[3] = {0, 0, 0}; // array indicating the position of a point relative to the minimum corner of the monitor LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) { @@ -1491,7 +1494,11 @@ std::vector dft_fields::fourier_sourcedata(const volume &wher temp_struct.idx_arr.push_back(idx); for (size_t i = 0; i < Nfreq; ++i) { EH0 = dJ_weight*dJ[reduced_grid_size*i+idx_1d]; - if (is_E_or_D(temp_struct.near_fd_comp)) EH0 *= -1; + + if (is_electric(c)) EH0 *= -1; + if (is_D(c) && f->fc->s->chi1inv[c - Dx + Ex][cd]) EH0 /= -f->fc->s->chi1inv[c - Dx + Ex][cd][idx]; + if (is_B(c) && f->fc->s->chi1inv[c - Bx + Hx][cd]) EH0 /= f->fc->s->chi1inv[c - Bx + Hx][cd][idx]; + EH0 /= f->S.multiplicity(ix0); temp_struct.amp_arr.push_back(EH0); } @@ -1503,7 +1510,11 @@ std::vector dft_fields::fourier_sourcedata(const volume &wher temp_struct.idx_arr.push_back(site_ind[j]); for (size_t i = 0; i < Nfreq; ++i) { EH0 = dJ_weight*dJ[reduced_grid_size*i+idx_1d]*0.25; // split the amplitude of the adjoint source into four parts - if (is_E_or_D(temp_struct.near_fd_comp)) EH0 *= -1; + + if (is_electric(c)) EH0 *= -1; + if (is_D(c) && f->fc->s->chi1inv[c - Dx + Ex][cd]) EH0 /= -f->fc->s->chi1inv[c - Dx + Ex][cd][idx]; + if (is_B(c) && f->fc->s->chi1inv[c - Bx + Hx][cd]) EH0 /= f->fc->s->chi1inv[c - Bx + Hx][cd][idx]; + EH0 /= f->S.multiplicity(ix0); temp_struct.amp_arr.push_back(EH0); }