tutorial and unit test for mode sources using DiffractedPlanewave object

[oskooi]

Closes #2107.

Unfortunately, I was not able to (again) fix the merge conflicts in #2107 and so just created another new PR.

I think the discussion in #2107 is now resolved by #2107 (comment).

This PR also includes some minor fixes (i.e., proper formatting) and additional comments to run_binary_grating_diffraction which is part of test_diffracted_planewave in python/tests/test_diffracted_planewave.py.

[codecov-commenter]

Codecov Report

Merging #2226 (82734d1) into master (a93a158) will increase coverage by 0.02%.
The diff coverage is 87.50%.

❗ Current head 82734d1 differs from pull request most recent head 4aa5e0b. Consider uploading reports for the commit 4aa5e0b to get more accurate results

@@            Coverage Diff             @@
##           master    #2226      +/-   ##
==========================================
+ Coverage   73.20%   73.23%   +0.02%     
==========================================
  Files          17       17              
  Lines        4904     4916      +12     
==========================================
+ Hits         3590     3600      +10     
- Misses       1314     1316       +2     

Impacted Files	Coverage Δ
python/simulation.py	76.85% <70.00%> (+0.01%)	⬆️
python/adjoint/objective.py	93.03% <100.00%> (+0.03%)	⬆️
python/adjoint/optimization_problem.py	57.29% <100.00%> (+0.22%)	⬆️

[oskooi]

The test_mode_source submodule (added as part of this PR) of tests_diffracted_planewave.py which involves using a DiffractedPlanewave object as the eig_band parameter of EigenModeSource is failing CI due to a segmentation fault. For some reason, the CI error occurs only when building Meep with MPI. The serial version passes.

This test passes on my local machine which builds Meep using MPI. Since fields::get_eigenmode has not been updated recently, it's unclear why MPI of the CI build is causing this test to fail.

[oskooi]

Looks like there was an issue with Python's garbage collection which can be resolved simply by putting Simulation.reset_meep() before the end of the function in order to manually trigger the garbage collection.

[mawc2019]

The segmentation fault may go away if import numpy as np is put after import meep as mp, i.e., using

import meep as mp
import numpy as np

instead of the opposite. In my experience, this sequence issue alone does not cause error, but may cause error in the presence of EigenModeSource. In the small example as follows, if the places of import numpy as np and import meep as mp are exchanged, the segmentation fault disappears. Alternatively, if EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),center=src_pt,size=mp.Vector3(0,sy)) is replaced with Source(src=mp.GaussianSource(fcen,fwidth=df),component=mp.Ez,center=src_pt,size=mp.Vector3(0,sy)), the segmentation fault also disappears.

import numpy as np
import meep as mp

resolution = 20
sx, sy = 4.0, 4.0
cell_size = mp.Vector3(sx,sy)

fcen, df = 1, 0.2
src_pt = mp.Vector3(-1)
sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),center=src_pt,size=mp.Vector3(0,sy))]
sim = mp.Simulation(resolution=resolution,cell_size=cell_size,sources=sources)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(), 1e-5))

[stevengj]

(Note that for any direction with zero thickness, the Cartesian and reciprocal bases are the same in that direction.)

[stevengj]

If you are worried about the initial guess as commented in #333, you could always change the initial guess manually. It shouldn't make any different in the mode that it converges to, it should only affect the number of Newton iterations.

[oskooi]

If you are worried about the initial guess as commented in #333, you could always change the initial guess manually. It shouldn't make any different in the mode that it converges to, it should only affect the number of Newton iterations.

This statement is only partially true in this example, unfortunately. No matter what I manually set as the initial guess for k inside fields::get_eigenmode, the Newton solver always finds a mode for which the $x$ component of the final k is $n\omega=(1.5)(2)=3$ and its $y$ component is the same as that of the initial k. This is why the correct mode is never found: the $y$ component (in the parallel direction) is simply being ignored by MPB's eigensolver routine. It's not obvious why this happens though.

Here are two examples of different initial ks and the output from MPB.

Example 1

initial k:

  k[0] = 2.68742 * R[0][0]; k[1] = 1.33333 * R[1][1]; k[2] = 0 * R[2][2];
  update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);

output of Newton solver and MPB's eigensolver:

MPB solved for frequency_1(2.68742,1.33333,0) = 1.79161 after 30 iters
MPB solved for frequency_1(3.31258,1.33333,0) = 2.20839 after 1 iters
MPB solved for frequency_1(3,1.33333,0) = 2 after 1 iters

Example 2

initial k:

  k[0] = 1.65873 * R[0][0]; k[1] = -0.89742 * R[1][1]; k[2] = 0 * R[2][2];
  update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);

output of Newton solver and MPB's eigensolver:

MPB solved for frequency_1(1.65873,-0.89742,0) = 1.11647 after 25 iters
MPB solved for frequency_1(4.33806,-0.89742,0) = 2.89613 after 1 iters
MPB solved for frequency_1(2.99197,-0.89742,0) = 2.00057 after 1 iters
MPB solved for frequency_1(2.99111,-0.89742,0) = 2 after 1 iters

[oskooi]

As a reference, here is a figure which illustrates the problem as well as a list of the main parameters and the expected result.

[stevengj]

(Note that R[0] should be (1,0,0), R[1] should be [0,Λ,0], and R[2] should be (0,0,1).)

You might try a few different ky values and see what kx values you get. It should be some function of the form kx = sqrt(a^2 - (b*ky)^2), and it might be revealing to learn what a and b are.

[oskooi]

There seems to be an inconsistency in the output from MPB related to the mode's wavevector compared to the actual mode that is found. (This inconsistency makes it difficult to debug this PR since we are relying on the output from MPB to determine whether it has found the correct mode or not.)

As a demonstration, this inconsistency can be seen in the output of python/examples/oblique-planewave.py in this PR which involves launching a planewave at an angle of 23.2° (rotated counterclockwise from the $+x$ axis) using an EigenModeSource in two different ways: (1) eig_band=1, direction=mp.NO_DIRECTION, eig_kpoint=k_point, eig_parity=mp.ODD_Z and (2) eig_band=mp.DiffractedPlanewave(...). The mode has k_point = (1.37870, 0.59091, 0.0) and $\omega=1.0$ in a homogenous medium with $n=1.5$. Note that the magnitude of k_point is $n\omega=1.5$.

The correct mode is launched in both cases as verified by the field snapshot (shown below). However, in the case of (1), the output from MPB is: MPB solved for frequency_1(1.63197,0.590913,0) = 1 after 68 iters. That is, MPB reports that the mode's wavevector is (1.63197,0.590913,0) rather than (1.37870, 0.59091, 0.0). The magnitude of MPB's reported wavevector is 1.73566 which is clearly not the expected result of 1.5.

Given the non-orthogonal lattice and reciprocal lattice basis vectors, perhaps the wavevector is simply not being converted to Cartesian coordinates correctly in the output?

meep/src/mpb.cpp

Lines 573 to 576 in 82734d1

	master_printf("MPB solved for frequency_%d(%g,%g,%g) "
	"= %g after %d iters\n",
	band_num, G[0][0] * k[0], G[1][1] * k[1], G[2][2] * k[2],
	sqrt(eigvals[band_num - 1]), num_iters);

1. EigenModeSource with eig_band=1

sources = [
    mp.EigenModeSource(
        src=mp.ContinuousSource(fsrc),
        center=mp.Vector3(),
        size=mp.Vector3(y=10),
        direction=mp.NO_DIRECTION,
        eig_kpoint=k_point,
        eig_parity=mp.ODD_Z,
    )
]

KPOINT: 1.5, 5.90913, 0
MPB solved for frequency_1(1.63197,0.590913,0) = 1 after 68 iters

For reference, the lattice R and reciprocal lattice G basis vectors are:

R[0] = (0.919135, 0.393942, 0)
R[1] = (0, 10, 0)
R[2] = (0, 0, 1)

G[0] = (1.08798, 0, 0)
G[1] = (-0.0428601, 0.1, 0)
G[2] = (0, 0, 1)

2. EigenModeSource with eig_band=mp.DiffractedPlanewave(...)

sources = [
    mp.EigenModeSource(
        src=mp.ContinuousSource(fsrc),
        center=mp.Vector3(),
        size=mp.Vector3(y=10),
        eig_band=mp.DiffractedPlanewave((0, 0, 0), mp.Vector3(0, 1, 0), 1, 0),
    )
]

KPOINT: 0, 5.90913, 0
NEW KPOINT: 1.5, 5.90913, 0

[oskooi]

Some additional information. It seems that while #2114 fixed one bug in fields::get_eigenmoderelated to oblique modes, it uncovered a new and separate bug which had somehow gone unnoticed. Prior to #2114, the example in #2054 (comment) which involves specifying an oblique mode using kpoint_func=lambda *not_used: kdiff in get_eigenmode_coefficients worked fine and the adjoint gradient agreed with the finite-difference result. However, this example started failing with #2114.

The problematic issue in the unit test of this PR as well as the example in #2054 (comment) involves launching an input planewave source with a given $\vec{k}$ (which requires specifying the k_point parameter of the Simulation object) and then trying to compute the mode coefficient for a different wavevector $\vec{k'}$. The key is to do the mode decomposition without specifying the kpoint_func parameter of get_eigenmode_coefficients.