From bc70d311380472c2a756702e6e79f48711a2da66 Mon Sep 17 00:00:00 2001 From: mawc2019 Date: Mon, 6 Jun 2022 23:12:10 -0400 Subject: [PATCH 1/5] deal with D and B fields --- src/dft.cpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) 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); } From a816bc9edd275cf1b117a82c274e2d4653cf19b0 Mon Sep 17 00:00:00 2001 From: mawc2019 Date: Thu, 16 Jun 2022 21:26:43 -0400 Subject: [PATCH 2/5] add unit test --- python/tests/test_adjoint_solver.py | 32 ++++++++++++++--------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 39fdbfab4..59d1e2a43 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,21 +505,21 @@ 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)) + self.assertClose(adjsol_obj,Dz2_unperturbed,epsilon=1e-6) - ## 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 self.assertClose(adj_scale,fd_grad,epsilon=tol) From 22cd1bce691fb93046af5505f69c5a31a6b1317c Mon Sep 17 00:00:00 2001 From: mawc2019 Date: Fri, 17 Jun 2022 00:13:02 -0400 Subject: [PATCH 3/5] adjust tolerance in unit test --- python/tests/test_adjoint_solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 59d1e2a43..719bcd739 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -510,7 +510,7 @@ def test_complex_fields(self): ## compare objective results print("Dz2 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,Dz2_unperturbed)) - self.assertClose(adjsol_obj,Dz2_unperturbed,epsilon=1e-6) + self.assertClose(adjsol_obj,Dz2_unperturbed,epsilon=1e-4 if mp.is_single_precision() else 1e-6) ## compute perturbed |Dz|^2 Dz2_perturbed = forward_simulation_complex_fields(p+dp, frequencies) From aaad41544e60cb8917f20950fb028eec062229dd Mon Sep 17 00:00:00 2001 From: mawc2019 Date: Sun, 19 Jun 2022 21:24:08 -0400 Subject: [PATCH 4/5] adjust tolerance further --- python/tests/test_adjoint_solver.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 719bcd739..29b3b6ff7 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -510,7 +510,8 @@ def test_complex_fields(self): ## compare objective results print("Dz2 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,Dz2_unperturbed)) - self.assertClose(adjsol_obj,Dz2_unperturbed,epsilon=1e-4 if mp.is_single_precision() else 1e-6) + tol = 1e-5 if mp.is_single_precision() else 1e-6 + self.assertClose(adjsol_obj,Dz2_unperturbed,epsilon=tol) ## compute perturbed |Dz|^2 Dz2_perturbed = forward_simulation_complex_fields(p+dp, frequencies) @@ -521,7 +522,7 @@ def test_complex_fields(self): adj_scale = (dp[None,:]@adjsol_grad).flatten() 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.055 if mp.is_single_precision() else 0.002 self.assertClose(adj_scale,fd_grad,epsilon=tol) def test_damping(self): From 50e6f3fea6a5362c176a8371dcbe39b6f2ad592b Mon Sep 17 00:00:00 2001 From: mawc2019 Date: Sun, 19 Jun 2022 23:19:19 -0400 Subject: [PATCH 5/5] adjust tolerance further --- python/tests/test_adjoint_solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 29b3b6ff7..53a97d123 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -522,7 +522,7 @@ def test_complex_fields(self): adj_scale = (dp[None,:]@adjsol_grad).flatten() fd_grad = Dz2_perturbed-Dz2_unperturbed print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad)) - tol = 0.055 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):