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.

Sign up for GitHub

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

Jump to bottom

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

Closed
wants to merge 2 commits into from
Closed

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fe094a0
tutorial and unit test for mode sources using DiffractedPlanewave object
oskooi Jun 16, 2022
6c7abb1
remove sim.reset_meep() at the end of run_mode_source
oskooi Jun 19, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 18 additions & 9 deletions doc/docs/Python_Tutorials/Eigenmode_Source.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -243,9 +243,9 @@ cell_size = mp.Vector3(14,10,0)
pml_layers = [mp.PML(thickness=2,direction=mp.X)]

# rotation angle (in degrees) of planewave, counter clockwise (CCW) around z-axis
rot_angle = np.radians(0)
rot_angle = np.radians(23.2)

fsrc = 1.0 # frequency of planewave (wavelength = 1/fsrc)
fsrc = 1.0 # frequency (wavelength = 1/fsrc)

n = 1.5 # refractive index of homogeneous material
default_material = mp.Medium(index=n)
Expand All @@ -258,8 +258,7 @@ sources = [mp.EigenModeSource(src=mp.ContinuousSource(fsrc),
direction=mp.AUTOMATIC if rot_angle == 0 else mp.NO_DIRECTION,
eig_kpoint=k_point,
eig_band=1,
eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
eig_match_freq=True)]
eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z)]

sim = mp.Simulation(cell_size=cell_size,
resolution=resolution,
Expand All @@ -269,16 +268,14 @@ sim = mp.Simulation(cell_size=cell_size,
default_material=default_material,
symmetries=[mp.Mirror(mp.Y)] if rot_angle == 0 else [])

sim.run(until=100)
sim.run(until=20)

nonpml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(10,10,0))

sim.plot2D(fields=mp.Ez,
output_plane=nonpml_vol)

if mp.am_master():
plt.axis('off')
plt.savefig('pw.png',bbox_inches='tight')
plt.axis('off')
plt.show()
```

Note that the line source spans the *entire* length of the cell. (Planewave sources extending into the PML region must include `is_integrated=True`.) This example involves a continuous-wave (CW) time profile. For a pulse profile, the oblique planewave is incident at a given angle for only a *single* frequency component of the source as described in [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](../Python_Tutorials/Basics.md#angular-reflectance-spectrum-of-a-planar-interface). This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency $\omega$, you will need to do separate simulations involving different values of $\vec{k}$ (`k_point`) since each set of $(\vec{k},\omega)$ specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707).
Expand All @@ -288,3 +285,15 @@ Shown below are the steady-state field profiles generated by the planewave for t
<center>
![](../images/eigenmode_planewave.png)
</center>

An alternative method to launching a planewave in homogeneous media which provides more control over specifying its polarization, particularly in 3D, is to use a `DiffractedPlanewave` object for the `eig_band` property of the `EigenModeSource` instead of a band number and parity as in the previous example which used 1 and `mp.ODD_Z`, respectively. A planewave with wavevector $\vec{k}$ can be defined using the *zeroth* diffraction order combined with the `k_point` of the `Simulation` object set to $\vec{k}$. As a demonstration, the previous example which involved an $E_z$-polarized source (equivalent to the $\mathcal{S}$ polarization given the plane of incidence of $xy$), can also be defined using:

```py
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))]
```
25 changes: 16 additions & 9 deletions python/examples/oblique-planewave.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,17 @@
import numpy as np
import matplotlib.pyplot as plt


resolution = 50 # pixels/μm

cell_size = mp.Vector3(14,10,0)

pml_layers = [mp.PML(thickness=2,direction=mp.X)]

# rotation angle (in degrees) of planewave, counter clockwise (CCW) around z-axis
rot_angle = np.radians(0)
rot_angle = np.radians(23.2)

fsrc = 1.0 # frequency of planewave (wavelength = 1/fsrc)
fsrc = 1.0 # frequency (wavelength = 1/fsrc)

n = 1.5 # refractive index of homogeneous material
default_material = mp.Medium(index=n)
Expand All @@ -24,8 +25,16 @@
direction=mp.AUTOMATIC if rot_angle == 0 else mp.NO_DIRECTION,
eig_kpoint=k_point,
eig_band=1,
eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
eig_match_freq=True)]
eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z)]

## equivalent definition
# 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))]

sim = mp.Simulation(cell_size=cell_size,
resolution=resolution,
Expand All @@ -35,13 +44,11 @@
default_material=default_material,
symmetries=[mp.Mirror(mp.Y)] if rot_angle == 0 else [])

sim.run(until=100)
sim.run(until=20)

nonpml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(10,10,0))

sim.plot2D(fields=mp.Ez,
output_plane=nonpml_vol)

if mp.am_master():
plt.axis('off')
plt.savefig('pw.png',bbox_inches='tight')
plt.axis('off')
plt.show()
106 changes: 99 additions & 7 deletions python/tests/test_diffracted_planewave.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,12 @@
import numpy as np


# Computes the mode coefficient of the transmitted orders of
# a binary grating given an incident planewave and verifies
# that the results are the same when using either a band number
# or `DiffractedPlanewave` object in `get_eigenmode_coefficients`.

class TestDiffractedPlanewave(unittest.TestCase):

@classmethod
def setUp(cls):
cls.resolution = 50 # pixels/um
cls.resolution = 50 # pixels/μm

cls.dpml = 1.0 # PML thickness
cls.dsub = 3.0 # substrate thickness
cls.dpad = 3.0 # length of padding between grating and PML
Expand All @@ -29,6 +25,12 @@ def setUp(cls):


def run_binary_grating_diffraction(self,gp,gh,gdc,theta):
"""Computes the mode coefficient of the transmitted orders of
a binary grating given an incident planewave using the two
approaches of band number and `DiffractedPlanewave` object in
`get_eigenmode_coefficients` and verifies that they produce
the same result.
"""
sx = self.dpml+self.dsub+gh+self.dpad+self.dpml
sy = gp
cell_size = mp.Vector3(sx,sy,0)
Expand Down Expand Up @@ -88,7 +90,7 @@ def _pw_amp(x):
m_minus = int(np.ceil((-self.fcen-k.y)*gp))

if theta == 0:
orders = range(m_plus)
orders = range(m_plus+1)
else:
# ordering of the modes computed by MPB is according to *decreasing*
# values of kx (i.e., closest to propagation direction of 0° or +x)
Expand Down Expand Up @@ -137,12 +139,102 @@ def _pw_amp(x):

print("PASSED.")


def test_diffracted_planewave(self):
self.run_binary_grating_diffraction(2.6,0.4,0.6,0)
self.run_binary_grating_diffraction(2.6,0.4,0.6,13.4)

# self.run_binary_grating_diffraction(10.0,0.5,0.5,0)
# self.run_binary_grating_diffraction(10.0,0.5,0.5,10.7)


def run_mode_source(self,gp,gh,gdc,m,diffpw):
"""Computes the transmitted flux of a
binary grating given an incident planewave
specified by the diffraction order `m` in the
y direction. The incident planewave is defined
using a mode source with either a band number
or `DiffractedPlanewave` object specified by
the boolean flag `diffpw`.
"""
sx = self.dpml+self.dsub+gh+self.dpad+self.dpml
sy = gp
cell_size = mp.Vector3(sx,sy,0)

ky = m/gp
theta = math.asin(ky/(self.fcen*self.ng))

# k (in source medium) with correct length (plane of incidence: XY)
k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(z=1), theta)

eig_parity = mp.ODD_Z
if theta == 0:
k = mp.Vector3()
eig_parity += mp.EVEN_Y
symmetries = [mp.Mirror(direction=mp.Y)]
else:
symmetries = []

src_pt = mp.Vector3(-0.5*sx+self.dpml,0,0)
if diffpw:
# the *zeroth* diffraction order defines a planewave with a
# wavevector equal to the `k_point` of the `Simulation` object
sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=0.1*self.fcen),
center=src_pt,
size=mp.Vector3(0,sy,0),
eig_band=mp.DiffractedPlanewave((0,0,0),
mp.Vector3(0,1,0),
1,
0))]
else:
sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=0.1*self.fcen),
center=src_pt,
size=mp.Vector3(0,sy,0),
direction=mp.NO_DIRECTION,
eig_kpoint=k,
Comment on lines +193 to +194
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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]?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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].

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

eig_band=1,
eig_parity=eig_parity)]

geometry = [mp.Block(material=self.glass,
size=mp.Vector3(self.dpml+self.dsub,mp.inf,mp.inf),
center=mp.Vector3(-0.5*sx+0.5*(self.dpml+self.dsub),0,0)),
mp.Block(material=self.glass,
size=mp.Vector3(gh,gdc*gp,mp.inf),
center=mp.Vector3(-0.5*sx+self.dpml+self.dsub+0.5*gh,0,0))]

sim = mp.Simulation(resolution=self.resolution,
cell_size=cell_size,
boundary_layers=self.pml_layers,
k_point=k,
geometry=geometry,
sources=sources,
symmetries=symmetries)

tran_pt = mp.Vector3(0.5*sx-self.dpml,0,0)
tran_flux = sim.add_flux(self.fcen,
0,
1,
mp.FluxRegion(center=tran_pt,
size=mp.Vector3(0,sy,0)))

sim.run(until_after_sources=mp.stop_when_fields_decayed(20,mp.Ez,src_pt,1e-6))

tran = mp.get_fluxes(tran_flux)[0]

return tran


def test_mode_source(self):
tran_bn = self.run_mode_source(1.5,0.5,0.3,2,False)
tran_dp = self.run_mode_source(1.5,0.5,0.3,2,True)
print("mode-source:, "
"{:.5f} (band number), "
"{:.5f} (diffraction order)".format(tran_bn,
tran_dp))
self.assertAlmostEqual(tran_bn,
tran_dp,
places=3 if mp.is_single_precision() else 4)


if __name__ == '__main__':
unittest.main()