antenna_freespace_cyl.py

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


parser = argparse.ArgumentParser()
parser.add_argument(
    '-res',
    type=float,
    help='resolution (pixels/μm)',
    default=100.,
)
args = parser.parse_args()

resolution = args.res  # pixels/um                                                                                                                          

s = 5.
dpml_r = 4.0
dpml_z = 2.0
cell_size = mp.Vector3(s+dpml_r,0,s+2*dpml_z)

pml_layers = [mp.PML(dpml_r,direction=mp.R),
              mp.PML(dpml_z,direction=mp.Z),]

fcen = 1.

sources = [mp.Source(
    src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
    center=mp.Vector3(),
    component=mp.Er)]

sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    dimensions=mp.CYLINDRICAL,
    m=-1,
    sources=sources,
    boundary_layers=pml_layers)

flux_box = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5*s,0,0.5*s),
        size=mp.Vector3(s,0,0)),
    mp.FluxRegion(
        center=mp.Vector3(s,0,0),
        size=mp.Vector3(0,0,s)),
    mp.FluxRegion(
        center=mp.Vector3(0.5*s,0,-0.5*s),
        size=mp.Vector3(s,0,0),
        weight=-1.0))

sim.run(mp.dft_ldos(fcen,0,1),
        until_after_sources=40)
        # until_after_sources=mp.stop_when_energy_decayed(                                                                                                  
        #     5.,                                                                                                                                           
        #     1e-7)) 
      
flux = mp.get_fluxes(flux_box)[0]

dV = 1/(resolution**2)
ldos1 = -np.real(sim.ldos_Fdata[0]*np.conj(sim.ldos_Jdata[0]))*dV

dV = 0.5/(resolution**3)
ldos2 = -np.real(sim.ldos_Fdata[0]*np.conj(sim.ldos_Jdata[0]))*dV

dV = 2*np.pi/(resolution**2)
ldos3 = -np.real(sim.ldos_Fdata[0]*np.conj(sim.ldos_Jdata[0]))*dV

dV = np.pi/(resolution**3)
ldos4 = -np.real(sim.ldos_Fdata[0]*np.conj(sim.ldos_Jdata[0]))*dV

print(f"flux:, {flux:.6f}, {ldos1:.6f}, {ldos2:.6f}, {ldos3:.6f}, {ldos4:.6f}")

sim.plot2D()
if mp.am_master():
    plt.savefig('cyl_layout.png',dpi=150,bbox_inches='tight')

import meep as mp
import numpy as np


resolution = 30  # pixels/um                                                                                                                          

s = 5.
dpml_r = 4.
dpml_z = 2.
cell_size = mp.Vector3(s+dpml_r,0,s+2*dpml_z)

pml_layers = [mp.PML(dpml_r,direction=mp.R),
              mp.PML(dpml_z,direction=mp.Z),]

fcen = 1.

sources = [mp.Source(
    src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
    center=mp.Vector3(),
    component=mp.Er)]

sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    dimensions=mp.CYLINDRICAL,
    m=-1,
    sources=sources,
    boundary_layers=pml_layers)

flux_plus_z = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5*s,0,0.5*s),
        size=mp.Vector3(s,0,0)))

flux_plus_r = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(s,0,0),
        size=mp.Vector3(0,0,s)))

flux_minus_z = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5*s,0,-0.5*s),
        size=mp.Vector3(s,0,0)))

for t in [100, 200, 400]:
    sim.run(until_after_sources=t)

    flux_plusz = mp.get_fluxes(flux_plus_z)[0]
    flux_plusr = mp.get_fluxes(flux_plus_r)[0]
    flux_minusz = mp.get_fluxes(flux_minus_z)[0]
    print(f"flux:, {sim.meep_time()}, {flux_plusz:.6f}, {flux_plusr:.6f}, {flux_minusz:.6f}")

t = 100

flux:, 150.0, -0.059667, 0.025234, 0.059667

t = 200

flux:, 250.0, 0.085347, 0.025234, -0.085347

t = 400

flux:, 450.0, -0.044804, 0.025234, 0.044804

If we shrink the size of the flux monitors in the $z$ direction to exclude $r=0$ then the flux is a constant value independent of the run time:

flux_plus_z = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5*s,0,0.5*s),
        size=mp.Vector3(0.95*s,0,0)))

flux_plus_r = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(s,0,0),
        size=mp.Vector3(0,0,s)))

flux_minus_z = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5*s,0,-0.5*s),
        size=mp.Vector3(0.95*s,0,0)))

 150.0, +0.022526, +0.025275, -0.022526
 250.0, +0.022526, +0.025275, -0.022526
 450.0, +0.022526, +0.025275, -0.022526

This suggests there is likely something wrong with the E_r or E_ɸ fields at $r=0$ and $z \neq 0$.

(Based on this, it also seems that using dV = np.pi/(resolution**3) to compute the total emitted power based on the LDOS via -np.real(sim.ldos_Fdata[0]*np.conj(sim.ldos_Jdata[0]))*dV matches the Poynting flux value.)

flux_plus_z = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
	center=mp.Vector3(0.01,0,0.5*s),
	size=mp.Vector3(0.02,0,0)))

for t in [100, 200, 400]:
    sim.run(until_after_sources=t)

    flux_plusz = mp.get_fluxes(flux_plus_z)[0]
    print(f"flux:, {sim.meep_time()}, {flux_plusz:.6f}")

flux:, 150.0, -0.038800
flux:, 250.0, 0.029357
flux:, 450.0, -0.031814