oskooi
commented
1 month ago I have put together a demonstration of this calculation which seems to work. The example involves surrounding the disc with $n = 2.4$ from the previous tutorial with an ITO coating of thickness 120 nm. However, the map of the absorbed power density is not that interesting: there is a prominent dark spot but is otherwise featureless.
import math
from typing import Tuple
import matplotlib.pyplot as plt
from meep.materials import SiO2, ITO
import meep as mp
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
DEBUG_OUTPUT = False
RESOLUTION_UM = 100
WAVELENGTH_UM = 0.55
N_DISC = 2.4
DISC_RADIUS_UM = 1.2
DISC_THICKNESS_UM = 0.7 * WAVELENGTH_UM / N_DISC
COATING_UM = 0.12
PML_UM = 1.0
PADDING_UM = 1.0
NUM_DIPOLES = 11
def plot_absorbed_power_density_map(
absorbed_power_density: np.ndarray,
filename: str
) -> None:
r_coord = np.linspace(
0, DISC_RADIUS_UM + COATING_UM, absorbed_power_density.shape[1]
)
z_coord = np.linspace(
0, DISC_THICKNESS_UM + COATING_UM, absorbed_power_density.shape[0]
)
print(f"shapes:, {r_coord.shape} (r_coord), {z_coord.shape} (z_coord), {absorbed_power_density.shape} (absorbed_power_density)")
fig, ax = plt.subplots()
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
im = ax.pcolormesh(
r_coord,
z_coord,
absorbed_power_density,
cmap="hot_r",
shading="gouraud",
)
ax.set_aspect("equal")
fig.colorbar(im, cax=cax)
if mp.am_master():
fig.savefig(filename, dpi=150, bbox_inches="tight")
def dipole_in_disc(zpos: float, rpos_um: float, m: int) -> Tuple[float, np.ndarray]:
"""Computes the total flux and absorbed power density of a dipole in a disc.
Args:
zpos: height of dipole above ground plane as fraction of disc thickness.
rpos_um: radial position of dipole.
m: angular φ dependence of the fields exp(imφ).
Returns:
The map of the absorbed power density in the disc with coating.
"""
frequency = 1 / WAVELENGTH_UM # center frequency of source/monitor
# Runtime termination criteria.
field_decay_threshold = 1e-8
size_r = DISC_RADIUS_UM + PADDING_UM + PML_UM
size_z = DISC_THICKNESS_UM + COATING_UM + PADDING_UM + PML_UM
cell_size = mp.Vector3(size_r, 0, size_z)
boundary_layers = [
mp.PML(PML_UM, direction=mp.R),
mp.PML(PML_UM, direction=mp.Z, side=mp.High),
]
src_pt = mp.Vector3(rpos_um, 0, -0.5 * size_z + zpos * DISC_THICKNESS_UM)
sources = [
mp.Source(
src=mp.GaussianSource(frequency, fwidth=0.05 * frequency),
component=mp.Er,
center=src_pt,
)
]
geometry = [
mp.Block(
material=mp.Medium(index=N_DISC),
center=mp.Vector3(
0.5 * DISC_RADIUS_UM, 0, -0.5 * size_z + 0.5 * DISC_THICKNESS_UM
),
size=mp.Vector3(DISC_RADIUS_UM, mp.inf, DISC_THICKNESS_UM)
),
mp.Block(
material=ITO,
center=mp.Vector3(
0.5 * (DISC_RADIUS_UM + COATING_UM),
0,
-0.5 * size_z + DISC_THICKNESS_UM + 0.5 * COATING_UM
),
size=mp.Vector3(DISC_RADIUS_UM + COATING_UM, mp.inf, COATING_UM)
),
mp.Block(
material=ITO,
center=mp.Vector3(
DISC_RADIUS_UM + 0.5 * COATING_UM,
0,
-0.5 * size_z + 0.5 * DISC_THICKNESS_UM
),
size=mp.Vector3(COATING_UM, mp.inf, DISC_THICKNESS_UM)
)
]
sim = mp.Simulation(
resolution=RESOLUTION_UM,
cell_size=cell_size,
dimensions=mp.CYLINDRICAL,
m=m,
boundary_layers=boundary_layers,
sources=sources,
geometry=geometry,
force_complex_fields=True,
)
dft_fields = sim.add_dft_fields(
[mp.Dr, mp.Dp, mp.Dz, mp.Er, mp.Ep, mp.Ez],
frequency,
0,
1,
center=mp.Vector3(
0.5 * (DISC_RADIUS_UM + COATING_UM),
0,
-0.5 * size_z + 0.5 * (DISC_THICKNESS_UM + COATING_UM)
),
size=mp.Vector3(
DISC_RADIUS_UM + COATING_UM,
0,
DISC_THICKNESS_UM + COATING_UM
),
yee_grid=False
)
sim.run(
mp.dft_ldos(frequency, 0, 1),
until_after_sources=mp.stop_when_fields_decayed(
25.0, mp.Er, src_pt, field_decay_threshold
),
)
delta_vol = 2 * np.pi * rpos_um / (RESOLUTION_UM**2)
dipole_power = (-np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) *
delta_vol)
dft_dr = sim.get_dft_array(dft_fields, mp.Dr, 0)
dft_dp = sim.get_dft_array(dft_fields, mp.Dp, 0)
dft_dz = sim.get_dft_array(dft_fields, mp.Dz, 0)
dft_er = sim.get_dft_array(dft_fields, mp.Er, 0)
dft_ep = sim.get_dft_array(dft_fields, mp.Ep, 0)
dft_ez = sim.get_dft_array(dft_fields, mp.Ez, 0)
absorbed_power_density = 2 * np.pi * frequency * np.imag(np.conj(dft_er) * dft_dr +
np.conj(dft_ep) * dft_dp +
np.conj(dft_ez) * dft_dz)
if DEBUG_OUTPUT:
fig, ax = plt.subplots()
sim.plot2D(ax=ax)
if mp.am_master():
fig.savefig("disc_coating_layout.png", dpi=150, bbox_inches="tight")
print(f"absorbed_power_density:, {absorbed_power_density.shape}")
plot_absorbed_power_density_map(absorbed_power_density, "absorbed_power_density")
return dipole_power, absorbed_power_density
if __name__ == "__main__":
dipole_height = 0.5
dipole_rpos_um = np.linspace(0, DISC_RADIUS_UM, NUM_DIPOLES)
delta_rpos_um = DISC_RADIUS_UM / (NUM_DIPOLES - 1)
num_m_per_dipole_rpos = np.zeros(NUM_DIPOLES - 1)
# 1. Er source at r = 0 requires a single simulation with m = ±1.
# An Er source at r = 0 needs to be slighty offset due to a bug.
# https://github.com/NanoComp/meep/issues/2704
dipole_rpos_um[0] = 1.5 / RESOLUTION_UM
m = -1
dipole_flux, dipole_absorbed_power = dipole_in_disc(
dipole_height,
dipole_rpos_um[0],
m,
)
flux_total = dipole_flux * dipole_rpos_um[0] * delta_rpos_um
absorbed_power_total = (
dipole_absorbed_power * dipole_rpos_um[0] * delta_rpos_um
)
print(f"dipole:, {dipole_rpos_um[0]:.4f}, {dipole_flux:.6f}")
# 2. Er source at r > 0 requires Fourier-series expansion of φ.
# Threshold flux to determine when to truncate expansion.
flux_decay_threshold = 1e-2
for i, rpos_um in enumerate(dipole_rpos_um[1:]):
dipole_absorbed_power_total = 0
dipole_flux_total = 0
dipole_flux_max = 0
m = 0
while True:
dipole_flux, dipole_absorbed_power = dipole_in_disc(
dipole_height, rpos_um, m
)
dipole_flux_total += dipole_flux * (1 if m == 0 else 2)
dipole_absorbed_power_total += dipole_absorbed_power * (
1 if m == 0 else 2
)
print(f"dipole:, {rpos_um:.4f}, {m}, {dipole_flux:.6f}")
if dipole_flux > dipole_flux_max:
dipole_flux_max = dipole_flux
if (
m > 0
and (dipole_flux / dipole_flux_max)
< flux_decay_threshold
):
num_m_per_dipole_rpos[i] = m
break
else:
m += 1
dipole_position_scale_factor = 0.5 * (dipole_rpos_um[0] / rpos_um) ** 2
flux_total += (
dipole_flux_total * dipole_position_scale_factor * rpos_um * delta_rpos_um
)
absorbed_power_total += (
dipole_absorbed_power_total
* dipole_position_scale_factor
* rpos_um
* delta_rpos_um
)
if mp.am_master():
np.savez(
"disc_absorbed_power.npz",
flux_total=flux_total,
absorbed_power_total=absorbed_power_total,
num_m_per_dipole_rpos=num_m_per_dipole_rpos
)
plot_absorbed_power_density_map(absorbed_power_total, "disc_absorbed_power")
stevengj
Might be nice to have an analogue of the "absorbed power density map of a lossy cylinder" tutorial but in cylindrical coordinates.
Basically the same calculation, except that the fields are sums over $m$, e.g. $\vec{E} = \sum_m \vec{E}_m(r,z) e^{im\phi}$. When computing the absorbed power $\sim \omega \Im[\vec{E}^ \cdot \vec{D}]$ at a particular $(r,\phi,z)$, you need to perform the sum before taking the dot product. However, if you are computing power averaged/integrated over all $\phi$, then the cross terms vanish and you can compute $\sum_m \omega \Im[\vec{E}_m^ \cdot \vec{D}_m]$, the sum of the dot products.
For example, you could do something similar to the "extraction efficiency of a collection of dipoles in a disc", but put an absorbing layer above the disk. Since the averaged absorbed power over an axisymmetric ensemble of dipoles is axisymmetric, you can use the simplified formula of computing the $\phi$-averaged power density for each dipole. Then you would sum over the different dipole radii using the same weighting factor as in the tutorial.
(Note that in general you want
yee_grid=Falseso that all field components are computed at the same location in order to combine them into a single number $\vec{E}^* \cdot \vec{D}$ at that location.)