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

tutorial and unit test for mode sources using DiffractedPlanewave object #2226

Open
wants to merge 5 commits into
base: master
Choose a base branch
from

Conversation

oskooi
Copy link
Collaborator

@oskooi oskooi commented Sep 6, 2022

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
Copy link

codecov-commenter commented Sep 6, 2022

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
Copy link
Collaborator Author

oskooi commented Sep 7, 2022

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
Copy link
Collaborator Author

oskooi commented Sep 7, 2022

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
Copy link
Contributor

mawc2019 commented Sep 7, 2022

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
Copy link
Collaborator

stevengj commented Sep 22, 2022

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

@stevengj
Copy link
Collaborator

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
Copy link
Collaborator Author

oskooi commented Sep 26, 2022

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
Copy link
Collaborator Author

oskooi commented Sep 26, 2022

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

k_point_and_diffracted_orders

k_point_and_diffracted_orders (1)

k_point_and_diffracted_orders (2)

@stevengj
Copy link
Collaborator

stevengj commented Sep 29, 2022

(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
Copy link
Collaborator Author

oskooi commented Oct 2, 2022

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);

planewave_oblique

  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)
  1. 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
Copy link
Collaborator Author

oskooi commented Oct 8, 2022

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.

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

Successfully merging this pull request may close these issues.

4 participants