Skip to content

Commit

Permalink
update Er source in tutorial on computing extraction efficiency in cy…
Browse files Browse the repository at this point in the history 
…lindrical coordinates (#2707)
  • Loading branch information
@oskooi
oskooi authored Nov 8, 2023
1 parent ebf6ea6 commit 91ccc5b
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 77 deletions.
109 changes: 51 additions & 58 deletions doc/docs/Python_Tutorials/Local_Density_of_States.md
View file
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 91ccc5b

Please sign in to comment.