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

Apply slight offset to Er source in tutorial on computing extraction efficiency in cylindrical coordinates #2707

Merged
merged 1 commit into from Nov 8, 2023
Merged

Apply slight offset to Er source in tutorial on computing extraction efficiency in cylindrical coordinates #2707

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
109 changes: 51 additions & 58 deletions doc/docs/Python_Tutorials/Local_Density_of_States.md
View file
Open in desktop
Expand Up @@ -329,45 +329,46 @@ To demonstrate this feature of the LDOS, we will compute the extraction efficien

The simulation setup is shown in the figures below for 3D Cartesian (cross section in $xz$) and cylindrical coordinates. (Note that this single-dipole calculation differs from the somewhat related flux calculation in [Tutorials/Custom Source/Stochastic Dipole Emission in Light Emitting Diodes](Custom_Source.md#stochastic-dipole-emission-in-light-emitting-diodes) which involved modeling a *collection* of dipoles.) In this example, the point-dipole source is positioned at $r=0$ which involves a single simulation. Nonaxisymmetric dipoles positioned at $r>0$, however, are more challenging because they involve multiple simulations. For a demonstration, see [Cylindrical Coordinates/Nonaxisymmetric Dipole Sources](Cylindrical_Coordinates.md#nonaxisymmetric-dipole-sources).

Note: because of a [bug](https://github.com/NanoComp/meep/issues/2704) for an $E_r$ point source at $r = 0$ and $m = \pm 1$ simulation, it is necessary to slightly offset the source to $r = 1.5\Delta r$. This incurs a small error which decreases linearly with resolution.

![](../images/dipole_extraction_eff_3D.png#center)

![](../images/dipole_extraction_eff_cyl.png#center)

The total emitted power obtained from the LDOS terms of the formula above must be multiplied by $\Delta V$, the volume of the voxel. In cylindrical coordinates, $\Delta V = \Delta r \times \Delta z \times 2 \pi r$. Meep implements an $r = 0$ source at $r = 0.5 \Delta r$, corresponding to the smallest-$r$ `Er` Yee grid point. This means that for a source at $r = 0$, $\Delta V = \pi / resolution^3$ since $\Delta r = \Delta z = 1 / resolution$. In 3D, $\Delta V = \Delta x \times \Delta y \times \Delta z = 1 / resolution^3$ for every voxel in the cell.
The total emitted power obtained from the LDOS terms of the formula above must be multiplied by $\Delta V$, the volume of the voxel. In cylindrical coordinates, $\Delta V = \Delta r \times \Delta z \times 2 \pi r$. Meep implements an $r = 0$ source at $r = 0.5 \Delta r$, corresponding to the smallest-$r$ $E_r$ Yee grid point. This means that for a source at $r = 0$, $\Delta V = \pi / resolution^3$ since $\Delta r = \Delta z = 1 / resolution$. In 3D, $\Delta V = \Delta x \times \Delta y \times \Delta z = 1 / resolution^3$ for every voxel in the cell.

As shown in the figure below, the results from the two coordinate systems have good agreement.

The simulation script is [examples/extraction_eff_ldos.py](https://github.com/NanoComp/meep/blob/master/python/examples/extraction_eff_ldos.py).

```py
import numpy as np
import meep as mp
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import meep as mp
import numpy as np


resolution = 80 # pixels/μm
dpml = 0.5 # thickness of PML
dair = 1.0 # thickness of air padding
L = 6.0 # length of non-PML region
n = 2.4 # refractive index of surrounding medium
wvl = 1.0 # wavelength (in vacuum)
dpml = 0.5 # thickness of PML
dair = 1.0 # thickness of air padding
L = 6.0 # length of non-PML region
n = 2.4 # refractive index of surrounding medium
wvl = 1.0 # wavelength (in vacuum)

fcen = 1 / wvl # center frequency of source/monitor
fcen = 1 / wvl # center frequency of source/monitor

# runtime termination criteria
tol = 1e-8


def extraction_eff_cyl(dmat: float, h: float) -> float:
"""Computes the extraction efficiency of a point dipole embedded
within a dielectric layer above a lossless ground plane in
cylindrical coordinates.
"""Computes the extraction efficiency in cylindrical coordinates.

Args:
dmat: thickness of dielectric layer.
h: height of dipole above ground plane as fraction of dmat.

Args:
dmat: thickness of dielectric layer.
h: height of dipole above ground plane as fraction of dmat.
Returns:
The extraction efficiency of the dipole within the dielecric layer.
"""
sr = L + dpml
sz = dmat + dair + dpml
Expand All @@ -379,7 +380,14 @@ def extraction_eff_cyl(dmat: float, h: float) -> float:
]

src_cmpt = mp.Er
src_pt = mp.Vector3(0, 0, -0.5 * sz + h * dmat)

# Because (1) Er is not defined at r=0 on the Yee grid, and (2) there
# seems to be a bug in the interpolation of an Er point source at r=0,
# the source is placed at r=~Δr (just outside the first voxel).
# This incurs a small error which decreases linearly with resolution.
# Ref: https://github.com/NanoComp/meep/issues/2704
src_pt = mp.Vector3(1.5 / resolution, 0, -0.5 * sz + h * dmat)

sources = [
mp.Source(
src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
Expand Down Expand Up @@ -422,37 +430,38 @@ def extraction_eff_cyl(dmat: float, h: float) -> float:

sim.run(
mp.dft_ldos(fcen, 0, 1),
until_after_sources=mp.stop_when_fields_decayed(
20, src_cmpt, src_pt, tol
),
until_after_sources=mp.stop_when_fields_decayed(20, src_cmpt, src_pt, tol),
)

out_flux = mp.get_fluxes(flux_air)[0]
dV = np.pi / (resolution**3)
if src_pt.x == 0:
dV = np.pi / (resolution**3)
else:
dV = 2 * np.pi * src_pt.x / (resolution**2)
total_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV
ext_eff = out_flux / total_flux
print(f"extraction efficiency (cyl):, "
f"{dmat:.4f}, {h:.4f}, {ext_eff:.6f}")
print(f"extraction efficiency (cyl):, " f"{dmat:.4f}, {h:.4f}, {ext_eff:.6f}")

return ext_eff


def extraction_eff_3D(dmat: float, h: float) -> float:
"""Computes the extraction efficiency of a point dipole embedded
within a dielectric layer above a lossless ground plane in
3D Cartesian coordinates.
"""Computes the extraction efficiency in 3D Cartesian coordinates.

Args:
dmat: thickness of dielectric layer.
h: height of dipole above ground plane as fraction of dmat.

Args:
dmat: thickness of dielectric layer.
h: height of dipole above ground plane as fraction of dmat.
Returns:
The extraction efficiency of the dipole within the dielecric layer.
"""
sxy = L + 2 * dpml
sz = dmat + dair + dpml
cell_size = mp.Vector3(sxy, sxy, sz)

symmetries = [
mp.Mirror(direction=mp.X, phase=-1),
mp.Mirror(direction=mp.Y)
mp.Mirror(direction=mp.Y),
]

boundary_layers = [
Expand Down Expand Up @@ -497,73 +506,57 @@ def extraction_eff_3D(dmat: float, h: float) -> float:
size=mp.Vector3(L, L, 0),
),
mp.FluxRegion(
center=mp.Vector3(
0.5 * L, 0, 0.5 * sz - dpml - 0.5 * dair
),
center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml - 0.5 * dair),
size=mp.Vector3(0, L, dair),
),
mp.FluxRegion(
center=mp.Vector3(
-0.5 * L, 0, 0.5 * sz - dpml - 0.5 * dair
),
center=mp.Vector3(-0.5 * L, 0, 0.5 * sz - dpml - 0.5 * dair),
size=mp.Vector3(0, L, dair),
weight=-1.0,
),
mp.FluxRegion(
center=mp.Vector3(
0, 0.5 * L, 0.5 * sz - dpml - 0.5 * dair
),
center=mp.Vector3(0, 0.5 * L, 0.5 * sz - dpml - 0.5 * dair),
size=mp.Vector3(L, 0, dair),
),
mp.FluxRegion(
center=mp.Vector3(
0, -0.5 * L, 0.5 * sz - dpml - 0.5 * dair
),
center=mp.Vector3(0, -0.5 * L, 0.5 * sz - dpml - 0.5 * dair),
size=mp.Vector3(L, 0, dair),
weight=-1.0,
),
)

sim.run(
mp.dft_ldos(fcen, 0, 1),
until_after_sources=mp.stop_when_fields_decayed(
20, src_cmpt, src_pt, tol
),
until_after_sources=mp.stop_when_fields_decayed(20, src_cmpt, src_pt, tol),
)

out_flux = mp.get_fluxes(flux_air)[0]
dV = 1 / (resolution**3)
total_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV
ext_eff = out_flux / total_flux
print(f"extraction efficiency (3D):, "
f"{dmat:.4f}, {h:.4f}, {ext_eff:.6f}")
print(f"extraction efficiency (3D):, {dmat:.4f}, {h:.4f}, {ext_eff:.6f}")

return ext_eff


if __name__ == "__main__":
layer_thickness = 0.7 * wvl / n
dipole_height = np.linspace(0.1,0.9,21)
dipole_height = np.linspace(0.1, 0.9, 21)

exteff_cyl = np.zeros(len(dipole_height))
exteff_3D = np.zeros(len(dipole_height))
for j in range(len(dipole_height)):
exteff_cyl[j] = extraction_eff_cyl(layer_thickness, dipole_height[j])
exteff_3D[j] = extraction_eff_3D(layer_thickness, dipole_height[j])

plt.plot(dipole_height,exteff_cyl,'bo-',label='cylindrical')
plt.plot(dipole_height,exteff_3D,'ro-',label='3D Cartesian')
plt.xlabel(f"height of dipole above ground plane "
f"(fraction of layer thickness)")
plt.plot(dipole_height, exteff_cyl, "bo-", label="cylindrical")
plt.plot(dipole_height, exteff_3D, "ro-", label="3D Cartesian")
plt.xlabel("height of dipole above ground plane (fraction of layer thickness)")
plt.ylabel("extraction efficiency")
plt.legend()

if mp.am_master():
plt.savefig(
'extraction_eff_vs_dipole_height.png',
dpi=150,
bbox_inches='tight'
)
plt.savefig("extraction_eff_vs_dipole_height.png", dpi=150, bbox_inches="tight")
```


Expand Down
52 changes: 33 additions & 19 deletions python/examples/extraction_eff_ldos.py
View file
Open in desktop
@@ -1,14 +1,13 @@
# Verifies that the extraction efficiency of a point dipole in a
# dielectric layer above a lossless ground plane computed in
# cylindrical and 3D Cartesian coordinates agree.
"""Computes the extraction efficiency in 3D and cylindrical coordinates.

import numpy as np
import meep as mp
import matplotlib
Verifies that the extraction efficiency of a point dipole in a dielectric
layer above a lossless metallic ground plane computed in two different
coordinate systems agree.
"""

matplotlib.use("agg")
import matplotlib.pyplot as plt

import meep as mp
import numpy as np

resolution = 80 # pixels/μm
dpml = 0.5 # thickness of PML
Expand All @@ -24,13 +23,14 @@


def extraction_eff_cyl(dmat: float, h: float) -> float:
"""Computes the extraction efficiency of a point dipole embedded
within a dielectric layer above a lossless ground plane in
cylindrical coordinates.
"""Computes the extraction efficiency in cylindrical coordinates.

Args:
dmat: thickness of dielectric layer.
h: height of dipole above ground plane as fraction of dmat.

Returns:
The extraction efficiency of the dipole within the dielecric layer.
"""
sr = L + dpml
sz = dmat + dair + dpml
Expand All @@ -42,7 +42,14 @@ def extraction_eff_cyl(dmat: float, h: float) -> float:
]

src_cmpt = mp.Er
src_pt = mp.Vector3(0, 0, -0.5 * sz + h * dmat)

# Because (1) Er is not defined at r=0 on the Yee grid, and (2) there
# seems to be a bug in the interpolation of an Er point source at r=0,
# the source is placed at r=~Δr (just outside the first voxel).
# This incurs a small error which decreases linearly with resolution.
# Ref: https://github.com/NanoComp/meep/issues/2704
src_pt = mp.Vector3(1.5 / resolution, 0, -0.5 * sz + h * dmat)

sources = [
mp.Source(
src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
Expand Down Expand Up @@ -89,7 +96,10 @@ def extraction_eff_cyl(dmat: float, h: float) -> float:
)

out_flux = mp.get_fluxes(flux_air)[0]
dV = np.pi / (resolution**3)
if src_pt.x == 0:
dV = np.pi / (resolution**3)
else:
dV = 2 * np.pi * src_pt.x / (resolution**2)
total_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV
ext_eff = out_flux / total_flux
print(f"extraction efficiency (cyl):, " f"{dmat:.4f}, {h:.4f}, {ext_eff:.6f}")
Expand All @@ -98,19 +108,23 @@ def extraction_eff_cyl(dmat: float, h: float) -> float:


def extraction_eff_3D(dmat: float, h: float) -> float:
"""Computes the extraction efficiency of a point dipole embedded
within a dielectric layer above a lossless ground plane in
3D Cartesian coordinates.
"""Computes the extraction efficiency in 3D Cartesian coordinates.

Args:
dmat: thickness of dielectric layer.
h: height of dipole above ground plane as fraction of dmat.

Returns:
The extraction efficiency of the dipole within the dielecric layer.
"""
sxy = L + 2 * dpml
sz = dmat + dair + dpml
cell_size = mp.Vector3(sxy, sxy, sz)

symmetries = [mp.Mirror(direction=mp.X, phase=-1), mp.Mirror(direction=mp.Y)]
symmetries = [
mp.Mirror(direction=mp.X, phase=-1),
mp.Mirror(direction=mp.Y),
]

boundary_layers = [
mp.PML(dpml, direction=mp.X),
Expand Down Expand Up @@ -182,7 +196,7 @@ def extraction_eff_3D(dmat: float, h: float) -> float:
dV = 1 / (resolution**3)
total_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV
ext_eff = out_flux / total_flux
print(f"extraction efficiency (3D):, " f"{dmat:.4f}, {h:.4f}, {ext_eff:.6f}")
print(f"extraction efficiency (3D):, {dmat:.4f}, {h:.4f}, {ext_eff:.6f}")

return ext_eff

Expand All @@ -199,7 +213,7 @@ def extraction_eff_3D(dmat: float, h: float) -> float:

plt.plot(dipole_height, exteff_cyl, "bo-", label="cylindrical")
plt.plot(dipole_height, exteff_3D, "ro-", label="3D Cartesian")
plt.xlabel(f"height of dipole above ground plane " f"(fraction of layer thickness)")
plt.xlabel("height of dipole above ground plane (fraction of layer thickness)")
plt.ylabel("extraction efficiency")
plt.legend()

Expand Down