Skip to content

Commit

Permalink
Fix FourierFields adjoint for D and B fields in isotropic media (#2095)
Browse files Browse the repository at this point in the history 
* deal with D and B fields

* add unit test

* adjust tolerance in unit test

* adjust tolerance further

* adjust tolerance further
  • Loading branch information
@mawc2019
mawc2019 authored Jun 23, 2022
1 parent 5aac79a commit 405b7e6
Show file tree
Hide file tree
Showing 2 changed files with 32 additions and 20 deletions.
35 changes: 18 additions & 17 deletions python/tests/test_adjoint_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -213,23 +213,23 @@ 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),
yee_grid=False)

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):
Expand All @@ -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,
Expand All @@ -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)
Expand Down Expand Up @@ -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):
Expand Down
17 changes: 14 additions & 3 deletions src/dft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1468,7 +1468,10 @@ std::vector<struct sourcedata> dft_fields::fourier_sourcedata(const volume &wher
std::vector<ptrdiff_t> idx_arr;
std::vector<std::complex<double> > amp_arr;
std::complex<double> EH0 = std::complex<double>(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) {
Expand All @@ -1491,7 +1494,11 @@ std::vector<struct sourcedata> 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);
}
Expand All @@ -1503,7 +1510,11 @@ std::vector<struct sourcedata> 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);
}
Expand Down

0 comments on commit 405b7e6

Please sign in to comment.