Incorrect gradients when materials contain conductivities

[smartalecH]

As described in #1788, the adjoint solver does not produce accurate gradients when materials contain conductivities (including artificial damping terms). This is why the following lines are commented out on master:

meep/src/meepgeom.cpp

Lines 603 to 606 in 7019c11

	// assign interpolated, nondimensionalized conductivity term
	// TODO: dampen the lorentzians to improve stability
	// mm->D_conductivity_diag.x = mm->D_conductivity_diag.y = mm->D_conductivity_diag.z = u*(1-u) *
	// omega_mean;

(Note that even if these lines were active, the actual conductivities themselves still need to be interpolated, as described in our paper -- but that's not why the gradients were failing, as the artificial damping terms introduced in #1788 also fail.)

When calculating Aᵤ (the derivative of the system matrix w.r.t. the design parameters), do a finite difference, performed here:

meep/src/meepgeom.cpp

Lines 2566 to 2592 in 7019c11

	/* For now we do a finite difference approach to estimate the
	gradient of the system matrix `A` since it's fairly accurate,
	cheap, and easy to generalize. */
	std::complex<double> dA_du_0[9] = {std::complex<double>(0, 0)};
	epsilon_material_grid(md, u - du);
	get_material_tensor(mm, freq, dA_du_0);
	std::complex<double> dA_du_1[9] = {std::complex<double>(0, 0)};
	epsilon_material_grid(md, u + du);
	get_material_tensor(mm, freq, dA_du_1);
	std::complex<double> dA_du[9] = {std::complex<double>(0, 0)};
	for (int i = 0; i < 9; i++)
	dA_du[i] = (dA_du_1[i] - dA_du_0[i]) / (2 * du);
	int dir_idx = 0;
	if (field_dir == meep::Ex || field_dir == meep::Er)
	dir_idx = 0;
	else if (field_dir == meep::Ey || field_dir == meep::Ep)
	dir_idx = 1;
	else if (field_dir == meep::Ez)
	dir_idx = 2;
	else
	meep::abort("Invalid adjoint field component");
	std::complex<double> result = fields_a * dA_du[3 * dir_idx + dir_idx] * fields_f;
	return result.real();

The code that evaluates ε(ω,u) is found here:

meep/src/meepgeom.cpp

Lines 2506 to 2533 in 7019c11

	void get_material_tensor(const medium_struct *mm, double freq,
	std::complex<double> *tensor) {
	/*
	Note that the current implementation assumes that any dispersion
	is either a lorentzian or drude profile.
	*/
	for (int i = 0; i < 9; i++) {
	std::complex<double> a = std::complex<double>(1, 0);
	std::complex<double> b = std::complex<double>(0, 0);
	// compute first part containing conductivity
	vector3 dummy;
	dummy.x = dummy.y = dummy.z = 0.0;
	double conductivityCur = vec_to_value(mm->D_conductivity_diag, dummy, i);
	a += std::complex<double>(0.0, conductivityCur / freq);
	// compute lorentzian component
	b = cvec_to_value(mm->epsilon_diag, mm->epsilon_offdiag, i);
	for (const auto &mm_susc: mm->E_susceptibilities) {
	meep::lorentzian_susceptibility sus = meep::lorentzian_susceptibility(
	mm_susc.frequency, mm_susc.gamma, mm_susc.drude);
	double sigma = vec_to_value(mm_susc.sigma_diag, mm_susc.sigma_offdiag, i);
	b += sus.chi1(freq, sigma);
	}
	// elementwise multiply
	tensor[i] = a * b;
	}
	}

What's interesting is gradients for dispersive materials without any global conductivities (e.g. LN in the material database) work just fine -- meaning at least part of the above code is working well...

Here is a test script that can be used during debugging.

(Here is another script that uses LN for the material grid and produces accurate gradients)

(cc @mochen4)

[stevengj]

Note that conductivity is handled in a very different place in the code from other dispersive materials — in the step_curl operation, rather than in update_edhb. Not sure if this would affect the adjoint formula?

[mochen4]

It seems that when damping is on, dA_du isn't calculated from finite difference but from the trivial case here:

meep/src/meepgeom.cpp

Lines 2553 to 2555 in 7019c11

	// trivial case
	if ((mm->E_susceptibilities.size() == 0) && mm->D_conductivity_diag.x == 0 &&
	mm->D_conductivity_diag.y == 0 && mm->D_conductivity_diag.z == 0){

Besides check additionally md->damping == 0 in the if statement above, shouldn't we also use 2*meep::pi*freq instead of freq here?

meep/src/meepgeom.cpp

Line 2519 in 7019c11

a += std::complex<double>(0.0, conductivityCur / freq);

Making these two changes seems to fix the issue.

[smartalecH]

shouldn't we also use 2meep::pifreq instead of freq here?

Yes that's a bug.

Making these two changes seems to fix the issue.

What do the results look like when you make those two changes?

[mochen4]

What do the results look like when you make those two changes?

With your test script above, the outputs are:
adjoint solver: [-2.46963186e-05], FD: [-2.46378424e-05]|| rel_error = [0.00237343]

I also ran the script a few times without fixing the random seed, and the rel_error is consistently around 1e-3.