tutorial for LDOS and extraction efficiency of an LED (#2215)

* tutorial for LDOS and extraction efficiency

* tweaks and fixes

* Update Local_Density_of_States.md

Co-authored-by: Steven G. Johnson <stevenj@mit.edu>

doc/docs/Python_User_Interface.md

@@ -88,20 +88,20 @@ control various parameters of the Meep computation.
 def __init__(self,
              cell_size: Union[meep.geom.Vector3, Tuple[float, ...], NoneType] = None,
              resolution: float = None,
-             geometry: Union[List[meep.geom.GeometricObject], NoneType] = None,
-             sources: Union[List[meep.source.Source], NoneType] = None,
+             geometry: Optional[List[meep.geom.GeometricObject]] = None,
+             sources: Optional[List[meep.source.Source]] = None,
              eps_averaging: bool = True,
              dimensions: int = 3,
-             boundary_layers: Union[List[meep.simulation.PML], NoneType] = None,
-             symmetries: Union[List[meep.simulation.Symmetry], NoneType] = None,
+             boundary_layers: Optional[List[meep.simulation.PML]] = None,
+             symmetries: Optional[List[meep.simulation.Symmetry]] = None,
              force_complex_fields: bool = False,
              default_material: meep.geom.Medium = Medium(),
              m: float = 0,
              k_point: Union[meep.geom.Vector3, Tuple[float, ...], bool] = False,
              kz_2d: str = 'complex',
-             extra_materials: Union[List[meep.geom.Medium], NoneType] = None,
-             material_function: Union[Callable[[Union[meep.geom.Vector3, Tuple[float, ...]]], meep.geom.Medium], NoneType] = None,
-             epsilon_func: Union[Callable[[Union[meep.geom.Vector3, Tuple[float, ...]]], float], NoneType] = None,
+             extra_materials: Optional[List[meep.geom.Medium]] = None,
+             material_function: Optional[Callable[[Union[meep.geom.Vector3, Tuple[float, ...]]], meep.geom.Medium]] = None,
+             epsilon_func: Optional[Callable[[Union[meep.geom.Vector3, Tuple[float, ...]]], float]] = None,
              epsilon_input_file: str = '',
              progress_interval: float = 4,
              subpixel_tol: float = 0.0001,
@@ -112,8 +112,8 @@ def __init__(self,
              num_chunks: int = 0,
              Courant: float = 0.5,
              accurate_fields_near_cylorigin: bool = False,
-             filename_prefix: Union[str, NoneType] = None,
-             output_volume: Union[meep.simulation.Volume, NoneType] = None,
+             filename_prefix: Optional[str] = None,
+             output_volume: Optional[meep.simulation.Volume] = None,
              output_single_precision: bool = False,
              geometry_center: Union[meep.geom.Vector3, Tuple[float, ...]] = Vector3<0.0, 0.0, 0.0>,
              force_all_components: bool = False,
@@ -2085,7 +2085,8 @@ arbitrarily spaced frequencies. One can also pass in an `Ldos` object as
 
 The resulting spectrum is outputted as comma-delimited text, prefixed by `ldos:,`, and
 is also stored in the `ldos_data` variable of the `Simulation` object after the `run`
-is complete.
+is complete. The Fourier-transformed electric field and current source are stored in
+the `ldos_Fdata` and `ldos_Jdata` of the `Simulation` object, respectively.
 
 </div>
 
@@ -6240,7 +6241,7 @@ def __init__(self,
              side: int = -1,
              R_asymptotic: float = 1e-15,
              mean_stretch: float = 1.0,
-             pml_profile: Callable[[float], float] = <function <lambda> at 0x7f3ad518a310>):
+             pml_profile: Callable[[float], float] = <function <lambda> at 0x7f3a89740ca0>):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -6979,7 +6980,7 @@ def __init__(self,
              size: Union[meep.geom.Vector3, Tuple[float, ...]] = Vector3<0.0, 0.0, 0.0>,
              direction: int = -1,
              weight: float = 1.0,
-             volume: Union[meep.simulation.Volume, NoneType] = None):
+             volume: Optional[meep.simulation.Volume] = None):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -7395,7 +7396,7 @@ def __init__(self,
              pt: Union[meep.geom.Vector3, Tuple[float, ...]] = None,
              fcen: float = None,
              df: float = None,
-             mxbands: Union[int, NoneType] = None):
+             mxbands: Optional[int] = None):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -7859,7 +7860,8 @@ arbitrarily spaced frequencies. One can also pass in an `Ldos` object as
 
 The resulting spectrum is outputted as comma-delimited text, prefixed by `ldos:,`, and
 is also stored in the `ldos_data` variable of the `Simulation` object after the `run`
-is complete.
+is complete. The Fourier-transformed electric field and current source are stored in
+the `ldos_Fdata` and `ldos_Jdata` of the `Simulation` object, respectively.
 
 </div>
 

python/examples/extraction_eff_ldos.py

@@ -0,0 +1,206 @@
+# 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.
+
+import numpy as np
+import meep as mp
+import matplotlib
+
+matplotlib.use("agg")
+import matplotlib.pyplot as plt
+
+
+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)
+
+fcen = 1 / wvl  # center frequency of source/monitor
+
+# source properties (cylindrical)
+df = 0.05 * fcen
+cutoff = 10.0
+src = mp.GaussianSource(fcen, fwidth=df, cutoff=cutoff)
+
+# 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.
+
+    Args:
+      dmat: thickness of dielectric layer.
+      h: height of dipole above ground plane as fraction of dmat.
+    """
+    sr = L + dpml
+    sz = dmat + dair + dpml
+    cell_size = mp.Vector3(sr, 0, sz)
+
+    boundary_layers = [
+        mp.PML(dpml, direction=mp.R),
+        mp.PML(dpml, direction=mp.Z, side=mp.High),
+    ]
+
+    src_cmpt = mp.Er
+    src_pt = mp.Vector3(0, 0, -0.5 * sz + h * dmat)
+    sources = [mp.Source(src=src, component=src_cmpt, center=src_pt)]
+
+    geometry = [
+        mp.Block(
+            material=mp.Medium(index=n),
+            center=mp.Vector3(0, 0, -0.5 * sz + 0.5 * dmat),
+            size=mp.Vector3(mp.inf, mp.inf, dmat),
+        )
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        dimensions=mp.CYLINDRICAL,
+        m=-1,
+        boundary_layers=boundary_layers,
+        sources=sources,
+        geometry=geometry,
+    )
+
+    flux_air = sim.add_flux(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
+            size=mp.Vector3(L, 0, 0),
+        ),
+        mp.FluxRegion(
+            center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * dair),
+            size=mp.Vector3(0, 0, dair),
+        ),
+    )
+
+    sim.run(
+        mp.dft_ldos(fcen, 0, 1),
+        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)
+    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}")
+
+    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.
+
+    Args:
+      dmat: thickness of dielectric layer.
+      h: height of dipole above ground plane as fraction of dmat.
+    """
+    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)]
+
+    boundary_layers = [
+        mp.PML(dpml, direction=mp.X),
+        mp.PML(dpml, direction=mp.Y),
+        mp.PML(dpml, direction=mp.Z, side=mp.High),
+    ]
+
+    src_cmpt = mp.Ex
+    src_pt = mp.Vector3(0, 0, -0.5 * sz + h * dmat)
+    sources = [
+        mp.Source(
+            src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
+            component=src_cmpt,
+            center=src_pt,
+        )
+    ]
+
+    geometry = [
+        mp.Block(
+            material=mp.Medium(index=n),
+            center=mp.Vector3(0, 0, -0.5 * sz + 0.5 * dmat),
+            size=mp.Vector3(mp.inf, mp.inf, dmat),
+        )
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        boundary_layers=boundary_layers,
+        sources=sources,
+        geometry=geometry,
+        symmetries=symmetries,
+    )
+
+    flux_air = sim.add_flux(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(0, 0, 0.5 * sz - dpml),
+            size=mp.Vector3(L, L, 0),
+        ),
+        mp.FluxRegion(
+            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),
+            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),
+            size=mp.Vector3(L, 0, dair),
+        ),
+        mp.FluxRegion(
+            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),
+    )
+
+    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}")
+
+    return ext_eff
+
+
+if __name__ == "__main__":
+    layer_thickness = 0.7 * wvl / n
+    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.ylabel("extraction efficiency")
+    plt.legend()
+
+    if mp.am_master():
+        plt.savefig("extraction_eff_vs_dipole_height.png", dpi=150, bbox_inches="tight")

python/simulation.py

@@ -5825,7 +5825,8 @@ def dft_ldos(*args, **kwargs):
 
     The resulting spectrum is outputted as comma-delimited text, prefixed by `ldos:,`, and
     is also stored in the `ldos_data` variable of the `Simulation` object after the `run`
-    is complete.
+    is complete. The Fourier-transformed electric field and current source are stored in
+    the `ldos_Fdata` and `ldos_Jdata` of the `Simulation` object, respectively.
     """
     ldos = kwargs.get("ldos", None)
     if ldos is None: