tutorial for absorbed power density in cylindrical coordinates

[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=False so 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.)

[oskooi]

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]

You can't see much here because you just see the singularity in the absorption from the emitter that is adjacen to the absorber.

Probably better to just replicate the tutorial with an absorbing cylinder and a planewave, but in cylindrical coordinates (where it will be a 1d simultation). This way you can also validate the results by comparing them to the full 2d simulation.