Hi,
I am trying to check some optical plots in the following paper:
https://www.nature.com/articles/s41598-020-68704-w
Here is a smaller/quicker-to-simulate version of the problem. I am
interested in silicon at long wavelength (up to 1.2 μm), and have fit
Lorentzians to Schinke et al data to get the following:
eps01=1.000000000000757e+00
cSi_sig1=1.742452232381217e+00
cSi_frq1=2.268393842625499e-03
cSi_gam1=1.000000000000000e-06
cSi_sig2=4.413478307846660e+00
cSi_frq2=5.563486517688901e-03
cSi_gam2=1.039533704134854e-05
cSi_sig3=4.346562778556655e-04
cSi_frq3=1.015803616840379e-03
cSi_gam3=4.436686697711304e-05
cSi_sig4=2.667407380976183e+00
cSi_frq4=3.019419346733225e-03
cSi_gam4=7.381034557205803e-06
cSi_sig5=1.874745406362289e+00
cSi_frq5=2.999546106465388e-03
cSi_gam5=2.153817716389790e-06
cSi_sig6=-4.082635925775152e-02
cSi_frq6=2.261035517057193e-03
cSi_gam6=5.061060071803692e-04
cSi_susc = [mp.LorentzianSusceptibility(frequency=cSi_frq1*1000,
gamma=cSi_gam1*2000, sigma=cSi_sig1),
mp.LorentzianSusceptibility(frequency=cSi_frq2*1000,
gamma=cSi_gam2*2000, sigma=cSi_sig2),
mp.LorentzianSusceptibility(frequency=cSi_frq3*1000,
gamma=cSi_gam3*2000, sigma=cSi_sig3),
mp.LorentzianSusceptibility(frequency=cSi_frq4*1000,
gamma=cSi_gam4*2000, sigma=cSi_sig4),
mp.LorentzianSusceptibility(frequency=cSi_frq5*1000,
gamma=cSi_gam5*2000, sigma=cSi_sig5),
mp.LorentzianSusceptibility(frequency=cSi_frq6*1000,
gamma=cSi_gam6*2000, sigma=cSi_sig6)]
cSi = mp.Medium(epsilon=eps01, E_susceptibilities=cSi_susc,
valid_freq_range=cSi_range)
I am interested in the absorption of an inverted pyramid photonic
crystal structure in the 1-1.2 μm wavelength range. If I start with just
a silicon slab:
zSi = dpml+2*monpad+dsi/2
siFilm = mp.Block(size=mp.Vector3(sx,sy,dsi),
center=mp.Vector3(0,0,zSi),
material=cSi)
geometry = [siFilm]
I get that my code converges and the absorption looks basically how I
expect it to. If I add a pyramid of air using the sidewall_angle argument:
zPyrAir = zSi + dsi/2
vertices = [mp.Vector3(sx*0.5,sy*0.5,zPyrAir),
mp.Vector3(sx*0.5,sy*-0.5,zPyrAir),
mp.Vector3(sx*-0.5,sy*-0.5,zPyrAir),
mp.Vector3(sx*-0.5,sy*0.5,zPyrAir)]
pyr_air = mp.Prism(vertices,
# remove small number so that vertices at top
don't overlap (otherwise crashes)
height=pyr_height_h-1e-10,
axis=mp.Vector3(0,0,-1),
sidewall_angle=-t_angle*pi/180,
material=mp.Medium(epsilon=1))
geometry = [siFilm, pyr_air]
I get that the fields diverge (albeit very slowly), at least so it
seems. Boundary conditions are k_point=Vector3(0,0,0) with PML in the Z
direction. My guess has been that the pyramid reflects light such that
it hits the PML at a glancing angle and this gives reflective artefacts.
I've hence played with PML parameters (thickness, profile, R_asymptotic)
and the distance between sources, monitors and structure, but so far
haven't gotten anything to work. Is perhaps the absorption too small for
MEEP to handle (i.e. leading to significant roundoff errors)? Is there a
workaround? What are some typical values to try if the defaults do not work?
I've attached code as well which contains all the details.
For what it's worth, I also find that with this structure I can only add
one symmetry (as discussed earlier on this mailing list), although my
results are slightly different (<1%) and the run time improvement is not
substantial.
Philip
"""
..moduleauthor:: jph
..date:: 2020-07-28
One minimal example to test MEEP long wavelength absorption and symmetry:
1 micron silicon slab. I am mostly checking for convergence, not necessarily
correctness
Sent version
"""
import sys
import meep as mp
from math import pi, tan, sin
import numpy as np
from matplotlib import pyplot as plt
cSi_range = mp.FreqRange(min=1/1.2, max=1/1.)
eps01=1.000000000000757e+00
cSi_sig1=1.742452232381217e+00
cSi_frq1=2.268393842625499e-03
cSi_gam1=1.000000000000000e-06
cSi_sig2=4.413478307846660e+00
cSi_frq2=5.563486517688901e-03
cSi_gam2=1.039533704134854e-05
cSi_sig3=4.346562778556655e-04
cSi_frq3=1.015803616840379e-03
cSi_gam3=4.436686697711304e-05
cSi_sig4=2.667407380976183e+00
cSi_frq4=3.019419346733225e-03
cSi_gam4=7.381034557205803e-06
cSi_sig5=1.874745406362289e+00
cSi_frq5=2.999546106465388e-03
cSi_gam5=2.153817716389790e-06
cSi_sig6=-4.082635925775152e-02
cSi_frq6=2.261035517057193e-03
cSi_gam6=5.061060071803692e-04
cSi_susc = [mp.LorentzianSusceptibility(frequency=cSi_frq1*1000, gamma=cSi_gam1*2000, sigma=cSi_sig1),
mp.LorentzianSusceptibility(frequency=cSi_frq2*1000, gamma=cSi_gam2*2000, sigma=cSi_sig2),
mp.LorentzianSusceptibility(frequency=cSi_frq3*1000, gamma=cSi_gam3*2000, sigma=cSi_sig3),
mp.LorentzianSusceptibility(frequency=cSi_frq4*1000, gamma=cSi_gam4*2000, sigma=cSi_sig4),
mp.LorentzianSusceptibility(frequency=cSi_frq5*1000, gamma=cSi_gam5*2000, sigma=cSi_sig5),
mp.LorentzianSusceptibility(frequency=cSi_frq6*1000, gamma=cSi_gam6*2000, sigma=cSi_sig6)]
cSi = mp.Medium(epsilon=eps01, E_susceptibilities=cSi_susc, valid_freq_range=cSi_range)
if __name__ == '__main__':
### INPUT PARAMETERS ###
is_pyr = int(sys.argv[1])
wvl_min = 1.0
wvl_max = 1.2
max_time = 2000
decay_factor = 1e-5
if is_pyr:
resolution = 50
outdir = 'test_simple_pyr'
else:
resolution = 50
outdir = 'test_slab'
dpml = 2.0
# PYRAMID #
pyr_base_L = 1
# pyr_height_h = 1
t_angle = 35.3
tilt_angle = 90-t_angle
pyr_height_h = tan(tilt_angle*pi/180)*pyr_base_L/2
dsi = 1 # 1 micron thick silicon
### DERIVED PARAMETERS ###
nfreq = int(1000*(wvl_max-wvl_min)+1) # number of frequency samples
sourcepad = wvl_max
monpad = wvl_max
freq_min = 1/wvl_max
freq_max = 1/wvl_min
freq_centre = (freq_max+freq_min)/2
freq_width = freq_max-freq_min
### DOMAIN ###
sx = pyr_base_L
sy = pyr_base_L
sz = 2*dpml + 2*sourcepad + monpad + dsi + 2*monpad
domain = mp.Vector3(sx, sy, sz)
### SOURCE DEFINITION ###
sources = [mp.Source(mp.GaussianSource(frequency=freq_centre, width=freq_width),
component=mp.Ex,
center=mp.Vector3(0,0,sz-dpml-sourcepad),
size=mp.Vector3(sx,sy,0))]
### SIMULATION PARAMETERS ###
k = mp.Vector3(0,0,0)
# it seems symmetry works only to an approximation
# and can only have one of them with PBCs
symm = [] # [mp.Mirror(mp.Y, phase=1)] #, mp.Mirror(mp.X, phase=-1)] # right now, no symmetries used
pml_layers = [mp.PML(dpml, direction=mp.Z, side=mp.ALL)]
### NORMALISATION SIM ###
sim = mp.Simulation(cell_size=domain,
boundary_layers=pml_layers,
geometry=[],
resolution=resolution,
k_point=k,
geometry_center=mp.Vector3(0,0,+0.5*sz),
symmetries=symm,
sources=sources)
### DEFINE FLUX REGIONS ###
refl_flux_region = mp.FluxRegion(center=mp.Vector3(0,0,sz-dpml-sourcepad-monpad),
size=mp.Vector3(sx,sy,0))
refl = sim.add_flux(freq_centre, freq_width, nfreq, refl_flux_region)
tran_flux_region = mp.FluxRegion(center=mp.Vector3(0,0,dpml+monpad),
size=mp.Vector3(sx,sy,0))
tran = sim.add_flux(freq_centre, freq_width, nfreq, tran_flux_region)
### RUN NORM SIMULATION ###
checkpt = mp.Vector3(0,0,dpml+monpad) # check at centre of tran flux
sim.run(until_after_sources=mp.stop_when_fields_decayed(25, mp.Ex, checkpt, 1e-10))
### GET FLUXES FOR SIMULATION ###
refl_flux = mp.get_fluxes(refl)
tran_flux = mp.get_fluxes(tran)
flux_freq = mp.get_flux_freqs(refl)
refl_norm_data = sim.get_flux_data(refl)
tran_norm_flux = mp.get_fluxes(tran)
### DEFINE GEOMETRY FOR SIMULATION ###
sim.reset_meep() # first must reset
zSi = dpml+2*monpad+dsi/2
siFilm = mp.Block(size=mp.Vector3(sx,sy,dsi),
center=mp.Vector3(0,0,zSi),
material=cSi)
if is_pyr:
zPyrAir = zSi + dsi/2
vertices = [mp.Vector3(sx*0.5,sy*0.5,zPyrAir),
mp.Vector3(sx*0.5,sy*-0.5,zPyrAir),
mp.Vector3(sx*-0.5,sy*-0.5,zPyrAir),
mp.Vector3(sx*-0.5,sy*0.5,zPyrAir)]
pyr_air = mp.Prism(vertices,
# remove small number so that vertices at top don't overlap (otherwise crashes!!)
height=pyr_height_h-1e-10,
axis=mp.Vector3(0,0,-1),
sidewall_angle=-t_angle*pi/180,
material=mp.Medium(epsilon=1))
geometry = [siFilm, pyr_air]
else:
geometry = [siFilm]
### DEFINE SIMULATION ###
sim = mp.Simulation(cell_size=domain,
boundary_layers=pml_layers,
geometry=geometry,
resolution=resolution,
k_point=k,
geometry_center=mp.Vector3(0,0,+0.5*sz),
symmetries=symm,
sources=sources)
### DEFINE FLUXES ###
refl_flux_region = mp.FluxRegion(center=mp.Vector3(0,0,sz-dpml-sourcepad-monpad),
size=mp.Vector3(sx,sy,0))
refl = sim.add_flux(freq_centre, freq_width, nfreq, refl_flux_region)
tran_flux_region = mp.FluxRegion(center=mp.Vector3(0,0,dpml+monpad),
size=mp.Vector3(sx,sy,0))
tran = sim.add_flux(freq_centre, freq_width, nfreq, tran_flux_region)
sim.load_minus_flux_data(refl, refl_norm_data)
### NAIVE VISUALISATION ###
plt.figure(dpi=150)
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(0,0,+0.5*sz), size=mp.Vector3(sx,sy,0)), eps_parameters={'cmap':'Set2','interpolation':'none'})
plt.savefig(outdir+'/pyrxy.png')
plt.close()
plt.figure(dpi=150)
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(0,0,+0.5*sz), size=mp.Vector3(sx,0,sz)), eps_parameters={'cmap':'Set2','interpolation':'none'})
plt.savefig(outdir+'/pyrxz.png')
plt.close()
plt.figure(dpi=150)
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(0,0,+0.5*sz), size=mp.Vector3(0,sy,sz)), eps_parameters={'cmap':'Set2','interpolation':'none'})
plt.savefig(outdir+'/pyryz.png')
plt.close()
### RUN SIMULATION ###
# for now, checkpt is same
sim.run(until_after_sources=[max_time, mp.stop_when_fields_decayed(25, mp.Ex, checkpt, decay_factor)])
### SAVE DATA ###
refl_flux = mp.get_fluxes(refl)
tran_flux = mp.get_fluxes(tran)
flux_freq = mp.get_flux_freqs(refl)
wl=np.array([])
Rs=np.array([])
Ts=np.array([])
for i in range(nfreq):
wl = np.append(wl, 1/flux_freq[i])
Rs = np.append(Rs,-refl_flux[i]/tran_norm_flux[i])
Ts = np.append(Ts,tran_flux[i]/tran_norm_flux[i])
print('saving output csv files...')
print(outdir+"/wl_pyr.csv")
np.savetxt(outdir+"/wl_pyr.csv", wl)
np.savetxt(outdir+"/rs_pyr.csv", Rs)
np.savetxt(outdir+"/ts_pyr.csv", Ts)
print('files saved')
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
