Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix adjoint gradient with conductivities #1823

Closed
wants to merge 0 commits into from

Conversation

mochen4
Copy link
Collaborator

@mochen4 mochen4 commented Nov 9, 2021

Fixes #1811 after PR #1801

@smartalecH smartalecH changed the title Fix adjoint gradient with connectivity Fix adjoint gradient with conductivities Nov 10, 2021
@smartalecH
Copy link
Collaborator

Would be good to add a simple test to test_adjoint_solver.py.

Maybe just borrow a simple test like this...

def test_adjoint_solver_eigenmode(self):
print("*** TESTING EIGENMODE ADJOINT ***")
## test the single frequency and multi frequency case
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.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-6)
## 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.01
self.assertClose(adj_scale,fd_grad,epsilon=tol)

...and retrofit it to include a damping term.

src/meepgeom.cpp Outdated
@@ -562,7 +562,7 @@ void epsilon_material_grid(material_data *md, double u) {
// Linearly interpolate dc epsilon values
cinterp_tensors(m1->epsilon_diag, m1->epsilon_offdiag, m2->epsilon_diag, m2->epsilon_offdiag,
&mm->epsilon_diag, &mm->epsilon_offdiag, u);

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice if you could eliminate these formatting changes from the PR. It makes review much easier and is important when comparing diffs.

Maybe check some of the default settings for your code editor -- it looks like it may be removing trailing whitespace automatically.

src/meepgeom.cpp Outdated
@@ -2612,6 +2612,46 @@ void eff_chi1inv_row_disp(meep::component c, std::complex<double> chi1inv_row[3]
case meep::NO_DIRECTION: chi1inv_row[0] = chi1inv_row[1] = chi1inv_row[2] = 0; break;
}
}
//Based on eff_chi1inv_row_disp above
void cond_inv(meep::component c, std::complex<double> condinv_row[3],
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would probably be better to refactor and use sigma_row where possible.

void geom_epsilon::sigma_row(meep::component c, double sigrow[3], const meep::vec &r) {

src/meepgeom.cpp Outdated

for (int i=0;i<3;i++) dA_du[i] = (row_1[i] - row_2[i])/(2*du);
return dA_du[dir_idx] * fields_f;
}
return dA_du[dir_idx] * fields_f * condinv_row[dir_idx]*condinv_row[dir_idx];
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't seem to be the right adjoint formulation.

For one, this would get called even during the off-diagonal calculations (where there are no sigma off-diagonal terms).

Also, shouldn't there be an additional term independent of the primary recombination step? Something like

Dᵀ_a ε dσ/du D_f + σDᵀ_a dε/du D_f

Obviously, there are more scaling terms, but the point is that there is an additional inner product because dA/du yields a sum of tensors (like in your notes). You can do all of these in a single step, but it probably shouldn't be done here.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Incorrect gradients when materials contain conductivities
2 participants