tutorial and unit test for mode sources using `DiffractedPlanewave` object #2107

oskooi · 2022-06-16T15:44:56Z

Closes #2072.

For some reason, I was not able to resolve the merge conflicts in python/tests/test_diffracted_planewave.py in #2072 during a rebase. This is why I just created this new PR.

oskooi · 2022-06-16T19:03:09Z

There is a test failure due to a segmentation fault related to garbage collection in python/tests/test_diffracted_planewave.py which I am not able to reproduce locally.

*** Process received signal ***
Signal: Segmentation fault (11)
Signal code:  (128)
Failing at address: (nil)
[ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x43090)[0x7fc847044090]
[ 1] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(+0x199d88)[0x7fc84739dd88]
[ 2] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(+0x119982)[0x7fc84731d982]
[ 3] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(+0x1993e4)[0x7fc84739d3e4]
[ 4] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(+0x203f0e)[0x7fc847407f0e]
[ 5] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(PyGC_Collect+0x55)[0x7fc847408005]
[ 6] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(Py_FinalizeEx+0x67)[0x7fc8473f99e7]
[ 7] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(Py_Exit+0xc)[0x7fc8473f9cec]
[ 8] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(+0x1fa543)[0x7fc8473fe543]
[ 9] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(PyErr_PrintEx+0x3a)[0x7fc8473fd92a]
[10] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(PyRun_SimpleFileExFlags+0x21d)[0x7fc8473fd4cd]
[11] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(+0x201811)[0x7fc847405811]
[12] /opt/hostedtoolcache/Python/3.7.13/x64/lib/libpython3.7m.so.1.0(_Py_UnixMain+0x40)[0x7fc8474054a0]
[13] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf3)[0x7fc847025083]
[14] /opt/hostedtoolcache/Python/3.7.13/x64/bin/python(_start+0x2e)[0x5635a7d9209e]
*** End of error message ***

codecov-commenter · 2022-06-19T02:27:20Z

Codecov Report

Merging #2107 (6c7abb1) into master (91d36ba) will increase coverage by 0.08%.
The diff coverage is n/a.

@@            Coverage Diff             @@
##           master    #2107      +/-   ##
==========================================
+ Coverage   73.18%   73.26%   +0.08%     
==========================================
  Files          17       17              
  Lines        4941     4941              
==========================================
+ Hits         3616     3620       +4     
+ Misses       1325     1321       -4

Impacted Files	Coverage Δ
python/simulation.py	`77.06% <0.00%> (+0.12%)`	⬆️
python/source.py	`96.21% <0.00%> (+0.75%)`	⬆️

stevengj · 2022-06-23T19:23:37Z

python/tests/test_diffracted_planewave.py

+                                    direction=mp.NO_DIRECTION,
+                                    eig_kpoint=k,


Suggested change

direction=mp.NO_DIRECTION,

eig_kpoint=k,

Should be able to delete these lines as I commented in #2072. (It's possible you were having problems before that are now fixed by #2114.)

If it's not working, make sure the MPB/Meep verbosity is > 1 and look at the NEW KPOINT and MPB solved for frequency verbose outputs — you should see that it is setting the parallel component of k by the boundary conditions, and varying the perpendicular component of k.

Looks like there is a potential bug because printing out NEW KPOINT shows that its value is clearly incorrect:

NEW KPOINT: 3, 2, 0

The y component (i.e., the periodic direction of the monitor monitor) should be 1.333... (i.e., 4/3) to match the k_point of the Simulation object. The reason it is not this value is due to these lines:

meep/src/mpb.cpp

Lines 426 to 433 in 405b7e6

for (int i = 0; i < 3; ++i) {

k[i] = dot_product(R[i], kcart); // convert k to reciprocal basis

// G = inverse of R transpose, via cross-product formula

cross_product(G[i], R[(i + 1) % 3], R[(i + 2) % 3]);

double GdotR = dot_product(G[i], R[i]);

for (int j = 0; j < 3; ++j)

G[i][j] /= GdotR;

}

In this example, R = (1,1.5,1) and kcart = (0,1.33333,0) which results in k = (0, 2.0, 0).

The fix seems to be to just replace k[i] = dot_product(R[i], kcart) with k[i] = kcart[i]?

No, that code is correct — MPB expects the k point in the reciprocal-lattice basis, so it has to convert from the Cartesian basis of Meep. (It does the same thing with the DiffractedPlanewave object.)

Instead, something weird is going on. You said the output is:

MPB solved for frequency_1(3,1.33333,0) = 2 after 29 iters

which means that it is taking the nonzero ky into account, but somehow it is getting the wrong kx — somehow, it is acting as though it is passing ky==0 to MPB, even though we are clearly printing that it is passing a nonzero k[1].

I think I have finally gotten to the bottom of this issue. TLDR: the current behavior is a feature and not a bug.

Currently, to launch oblique planewaves using an EigenModeSource with eig_band set to an integer (rather than an DiffractedPlanewave object), the user must provide its wavevector k (a meep.Vector3 object) via direction=mp.NO_DIRECTION and eig_kpoint=k. This is because even though internally Meep sets the component of k parallel to the line source (in the directions with periodic boundaries) to the components of the k_point of the Simulation object, the correct perpendicular component $\left|k_\perp\right| = (n\omega)^2 - \left|k_\parallel\right|$ is only computed automatically for the case when eig_band is a DiffractedPlanewave object.

It would be nice to extend this feature to when eig_band is an integer but that can be a separate PR since it is not directly related to this PR.

No. A DiffractedPlanewave object only enables the analytical calculation of k. Otherwise, it uses a Newton solver to find the k where ω(k) = specified frequency. This of course should be the same thing (but slower, and with more ambiguity about polarization and more difficulty in choosing a specific diffraction order).

The wavevector of the oblique planewave in Cartesian coordinates should be: k = (2.68742, 1.33333, 0.00000). However, the Newton solver shows:

MPB solved for frequency_1(3,1.33333,0) = 2 after 29 iters update_maxwell_data_k:, k=(3.00000, 2.00000, 0.00000) update_maxwell_data_k:, G[0]=(1.00000, 0.00000, 0.00000) update_maxwell_data_k:, G[1]=(0.00000, 0.66667, 0.00000) update_maxwell_data_k:, G[2]=(0.00000, 0.00000, 1.00000)

Note that the $y$ component (parallel direction) of the MPB mode (1.33333) is correct but the $x$ component (perpendicular direction) (3) is not 2.68742.

$|k|=n\omega=(1.5)(2)=3$ but clearly the wavevector from the Newton solver is wrong since $|3,1.3333,0)| > 3$.

Turns out that there was a bug in the output of MPB's mode wavevector which is fixed by #2253. (The mode itself was correct.)

With the bug fix, the output from run_mode_source of test_diffracted_planewave.py in this PR is:

MPB solved for frequency_1(2.68742,1.33333,0) = 2 after 23 iters

This means k = (2.68742,1.33333,0) which is equivalent to the k_point of the Simulation class object as expected.

oskooi · 2022-09-06T05:14:00Z

Replaced by #2226.

tutorial and unit test for mode sources using DiffractedPlanewave object

fe094a0

oskooi added enhancement documentation labels Jun 16, 2022

remove sim.reset_meep() at the end of run_mode_source

6c7abb1

stevengj reviewed Jun 23, 2022

View reviewed changes

This was referenced Sep 6, 2022

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

Closed

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

Open

oskooi closed this Sep 6, 2022

oskooi deleted the diffpw_mode_source2 branch September 6, 2022 05:14

oskooi mentioned this pull request Sep 20, 2022

WIP: mode-solver improvements #333

Merged

	for (int i = 0; i < 3; ++i) {
	k[i] = dot_product(R[i], kcart); // convert k to reciprocal basis
	// G = inverse of R transpose, via cross-product formula
	cross_product(G[i], R[(i + 1) % 3], R[(i + 2) % 3]);
	double GdotR = dot_product(G[i], R[i]);
	for (int j = 0; j < 3; ++j)
	G[i][j] /= GdotR;
	}

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

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

Uh oh!

Conversation

oskooi commented Jun 16, 2022 • edited by stevengj Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

oskooi commented Jun 16, 2022

Uh oh!

codecov-commenter commented Jun 19, 2022

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

oskooi commented Sep 6, 2022

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

tutorial and unit test for mode sources using `DiffractedPlanewave` object #2107

tutorial and unit test for mode sources using `DiffractedPlanewave` object #2107

oskooi commented Jun 16, 2022 •

edited by stevengj

Loading