more accurate Meep eigenmodes using meepgeom.cpp for MPB

[lxvm]

This issue addresses an accuracy issue (not a convergence issue) observed when calculating get_eigenmode_coefficient along an x-invariant waveguide and observing much larger-than-anticipated errors compared to manually-calculated eigenmode overlap integrals.
This continues some of the observations from #3096.

The simulation describes the geometry below and studies the variation of the power in various measurements of the eigenmode as a function of the x-dependence of the mode monitor.

The difference in powers is shown here:

This is for the Ez-polarized mode and I observed second-order convergence (see #3096), but the error may be different for the other polarization.

The code is in this notebook:

waveguide-debug.ipynb

From Meep devs discussion:

We may want to test that the material interpolation as passed to MPB is x-invariant. The hypothesis is that it isn't and so in addition to the interpolation errors observed in the manual overlap integrals, the eigenmode profile computed by Meep is introducing an additional source of error.

In the Python interface, the geometry data is stored in the Simulation.geps field that has methods to fetch the epsilon tensor at each point. We can use this to give MPB the data it wants.

May also need to refactor geom_epsilon::eff_chi1inv_row to return the full epsilon inverse matrix.

[lxvm]

In the attached notebook I found that the waveguide mode returned by the Simulation.get_eigenmode method is x-invariant and reproduces the manual overlap integral result from before that used a mode profile obtained with the MPB module.

This notebook is simpler because it doesn't call the MPB module and it also runs significantly faster because it puts all of the mode monitors into the same simulation:
waveguide-debug-eigenmode.ipynb

Since fields::get_eigenmode_coefficients internally calls fields::get_eigenmode,

meep/src/mpb.cpp

Lines 953 to 955 in 76ed307

	void *mode_data =
	get_eigenmode(flux.freq[nf], d, flux.where, eig_vol, band_num, kpoint, match_frequency,
	parity, eig_resolution, eigensolver_tol, kdom, (void **)&mdata, dp);

I would think that the error is happening downstream in the get_overlap calls or the normalization, although I don't have a way to test those from the Python interface

[stevengj]

Looking at the code, for the eigenmode source it uses add_dft_flux to create the DFT objects, but the weights look a little weird — it includes the dV and interpolation weights for the E field but not for the H field.

This is the right thing to do for computing flux because we don't want to include the weights twice when computing $E \times H$. However, for the mode-overlap integral, the E and H fields from Meep are not multiplied by each other but rather by the MPB fields. That would suggest that you should include the dV and interpolation weights in both fields for that case.

As a quick hack, can you try changing the 5th argument from false to true for H = add_dft(...) on this line:

meep/src/dft.cpp

Lines 631 to 632 in e407a5a

	H = add_dft(cH[i], where->v, freq, Nfreq, false, 1.0, H, false, std::complex<double>(1.0, 0),
	centered_grid, 0, decimation_factor);

(However, I find it hard to believe that this is the error, since it seems like we would have noticed it before if the H integral was completely wrong! But maybe a dV factor is included somewhere else in the post-processing.)

[stevengj]

Either that, or set the corresponding argument to false in the preceding E = add_dft(...) line? However, it seems like the weights might be correct as-is.

The mode overlap integral is computed on this line:

meep/src/dft.cpp

Line 1041 in e407a5a

integral += w * mode1val * dft_val;

It uses a weight w, that includes both dV and the interpolation weight, computed on this line:

meep/src/dft.cpp

Line 995 in e407a5a

double w = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2);

This might seems like it is including the weight twice for the E field, but it checks whether the chunk is already including the weight and if so divides it out on this line:

meep/src/dft.cpp

Lines 1004 to 1005 in e407a5a

	if (include_dV_and_interp_weights && dft_val != 0.0)
	dft_val /= (sqrt_dV_and_interp_weights ? sqrt(w) : w);

• If that's functioning correctly, it shouldn't matter whether you include the weights or not during the DFT calculation because it will take them into account. This is a good sanity check: change the include_dV_and_interp_weights flag in add_dft, and it shouldn't change the eigenmode coefficient at all (up to roundoff error). @lxvm, can you try this?

[lxvm]

Sure, I'll try it out and report back here when I do. Thank you for the helpful guidance!

[stevengj]

@lxvm, was there any update on this, I forget?

[lxvm]

I just tried changing the include_dV_and_interp_weights flag in add_dft line 631 of dft.cpp to true and got the same answer. I hope I rebuilt Meep right when doing this (I used the --enable-maintainer-mode flag when reconfiguring and restarted the jupyter notebook)

[stevengj]

@lxvm, can you file a separate issue just on the fact that a manual eigenmode-overlap calculation is not matching the built-in eigenmode-coefficients, even with what should be the same get_eigenmode call? In the issue, give the code for the two cases that are supposed to have the same results but don't.

(@oskooi wants to give this example to Claude to see if it can debug it.)