Why is it difficult for my field to converge

[niubiworker]

After multiple corrections, the program can run, but I found that in the output txt document, it seems difficult for the field to converge to the expected value, which is to decay to e-3. I don't know why? Did my structure resonate?

test5-10-0.1.txt
code:

import numpy as np
import meep as mp
import matplotlib.pyplot as plt
from meep.materials import LiNbO3

a = 0.775
r,h,wl, = 0.315*a,0.3,1.55
period_num,y_num = 50,4
resolution,courant = 10,0.1
fcen,nfreq= 1/wl,10
df = 0.1*fcen
sx,sy,sz = 60,20*h,60
pml = 0.5*wl

def lno(a):
    b = 1+2.6734*a**2/(a**2-0.01764)+1.2290*a**2/(a**2-0.0591)+12.614*a**2/(a**2-474.60)
    return b
def lne(a):
    b = 1+2.9804*a**2/(a**2-0.02047)+0.5981*a**2/(a**2-0.0666)+8.9543*a**2/(a**2-416.08)
    return b

wg = 2*y_num*a*np.sqrt(3)
geometry = [mp.Block(center=mp.Vector3(0,0,0),
                     size=mp.Vector3(mp.inf,h,wg),
                     material=LiNbO3.epsilon(1/wl))]
geometry1 = [mp.Block(center=mp.Vector3(0,0,0),
                     size=mp.Vector3(mp.inf,h,wg),
                     material=mp.Medium(epsilon_diag=(lno(wl),lno(wl),lne(wl))))]
for i in np.linspace(-period_num*a,period_num*a,period_num*2+1):
    for j in np.linspace(-y_num*np.sqrt(3)*a,y_num*np.sqrt(3)*a,2*y_num+1):
        if j != 0:
            geometry1 += [mp.Cylinder(center=mp.Vector3(i,0,j,), 
                                      height=mp.inf, 
                                      material=mp.air, 
                                      radius=r,
                                      axis=mp.Vector3(y=1))]
            geometry1 += [mp.Cylinder(center=mp.Vector3(i+0.5*a,0,np.sqrt(3)/2*a+j,), 
                                      height=mp.inf, 
                                      material=mp.air,
                                      radius=r,
                                      axis=mp.Vector3(y=1))]
        else:
            geometry1 += [mp.Cylinder(center=mp.Vector3(i+0.5*a,0,np.sqrt(3)/2*a), 
                                      height=mp.inf, 
                                      material=mp.air,
                                      radius=r,
                                      axis=mp.Vector3(y=1))]  
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=df),
                     center=mp.Vector3(x=-0.5*sx+2*pml),
                     size=mp.Vector3(0,h,wg),
                     component=mp.Ex,)]
eig_sources = [mp.EigenModeSource(
    src=mp.GaussianSource(fcen,fwidth=df),
    center=mp.Vector3(-0.5*sx+2*pml,0,0),
    direction=mp.AUTOMATIC,
    size=mp.Vector3(0,h,wg),
    eig_band=1,
    eig_kpoint=mp.Vector3(0.5,0,0),
    eig_parity=mp.EVEN_Y,
    eig_lattice_center=mp.Vector3(-0.5*sx+2*pml-0.5*a,0,0),
    eig_lattice_size=mp.Vector3(a,h,2*wg),
    eig_match_freq=False
    )]
sim = mp.Simulation(cell_size=mp.Vector3(sx,sy,sz),
                    resolution=resolution,
                    boundary_layers=[mp.Absorber(pml)],
                    sources=eig_sources,
                    geometry=geometry,
                    Courant=courant,
                    symmetries=[mp.Mirror(mp.Y,phase=1),mp.Mirror(mp.Z,phase=-1)])
refl_fr = mp.FluxRegion(center=mp.Vector3(x=-0.5*sx+2*pml+a),
                        direction=mp.X,
                        size=mp.Vector3(0,h,2*wg))
refl = sim.add_flux(fcen,df,nfreq,refl_fr)
tran_fr = mp.FluxRegion(center=mp.Vector3(x=0.5*sx-pml),
                        direction=mp.X,
                        size=mp.Vector3(0,h,2*wg))
tran = sim.add_flux(fcen,df,nfreq,tran_fr)
'''sim.run(until = 1)
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(),size=mp.Vector3(sx,0,sz)))
plt.savefig("1.png")
plt.show()'''

pt = mp.Vector3(0.5*sx-pml-a,0,0)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ex,pt,1e-3))
straight_refl_data = sim.get_flux_data(refl)
straight_tran_flux = mp.get_fluxes(tran)

sim.reset_meep()

sim = mp.Simulation(cell_size=mp.Vector3(sx,sy,sz),
                    resolution=resolution,
                    boundary_layers=[mp.Absorber(pml)],
                    sources=eig_sources,
                    geometry=geometry1,
                    Courant=courant,
                    symmetries=[mp.Mirror(mp.Y,phase=1),mp.Mirror(mp.Z,phase=-1)])
ref1 = sim.add_flux(fcen,df,nfreq,refl_fr)
tran_fr = mp.FluxRegion(center=mp.Vector3(-0.5*sx+2*pml+period_num*a,0,0),
                        direction=mp.X,
                        size=mp.Vector3(0,2*h,2*wg))
tran = sim.add_flux(fcen,df,nfreq,tran_fr)
sim.load_minus_flux_data(ref1,straight_refl_data)
pt = mp.Vector3(-0.5*sx+2*pml+period_num*a-a,0,0)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, pt, 1e-3))
pc_ref1_flux = mp.get_fluxes(ref1)
pc_tran_flux = mp.get_fluxes(tran)
flux_freqs = mp.get_flux_freqs(ref1)
wl = []
Rs = []
Ts = []
for i in range(nfreq):
    wl = np.append(wl, 1 / flux_freqs[i])
    Rs = np.append(Rs, -pc_ref1_flux[i] / straight_tran_flux[i])
    Ts = np.append(Ts, pc_tran_flux[i] / straight_tran_flux[i])
plt.plot(wl, Rs, "bo-", label="reflectance")
plt.plot(wl, Ts, "ro-", label="transmittance")
plt.plot(wl, 1 - Rs - Ts, "go-", label="loss")
plt.xlim(min(wl),max(wl))
plt.xlabel("wavelength (μm)")
plt.legend(loc="upper right")
plt.savefig("test5-10-0.1-1.png")
'''def log1(base,x):
    return np.log(x)/np.log(base)
print(-10*log1(10,Ts))'''
plt.show()

[smartalecH]

Did my structure resonate?

From your plot, it certainly looks that way. I haven't looked at your script in detail, but it appears you are doing some sort of photonic crystal, which should have plenty of resonances.

This is a good usecase for harminv. Maybe check out this tutorial: https://meep.readthedocs.io/en/latest/Python_Tutorials/Resonant_Modes_and_Transmission_in_a_Waveguide_Cavity/#resonant-modes