Open 0xInfty opened 3 years ago
stevengj
commented
3 years ago Are you sure you're using the same material parameters in your Mie calculation? Units and dispersive media are tricky…
Have you checked that the results are converged with resolution? You might need a high resolution to resolve the skin depth of gold. (Metals in FDTD are often a challenge for this reason.)
oskooi
commented
2 years ago Sorry for the long overdue reply. There are several challenging aspects in setting up the Meep simulation to produce results which are accurate enough to agree with Mie scattering theory when lossy materials (such as Au) are involved. This is mainly due to the surface plasmon polaritons on the metal-dielectric interface which is a unique feature not found in the tutorial example which involved only lossless materials. To give just one example, aside from cranking up the grid resolution to resolve the skin depth of Au and increasing the runtime until the surface plasmons have fully decayed away, the size of the closed box of flux monitors must also be enlarged (compared to the original example in the tutorial) such that the monitors are not touching the Au sphere. Generally, the flux monitors need to be positioned sufficiently far from the sphere such that the evanescent tails of the surface plasmons have decayed away. Otherwise, discretization artifacts will be exacerbated which could affect the accuracy of the flux measurements. To demonstrate this effect, here are results for two different runs: with and without "padding" of the flux box around the sphere. A slight padding of the flux box has a large impact on the results particularly at the longer wavelengths (more padding leads to better agreement with the theory). More generally though, the uniform Yee grid of FDTD is not well suited for problems involving a large relative mismatch in the dimensions of the geometry and field oscillations.
rom meep.materials import Au
import meep as mp
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import PyMieScatt as ps
n_H2O = 1.333 # refractive index of water
r = 1.0 # radius of sphere
wvl_min = 2*np.pi*r/10
wvl_max = 2*np.pi*r/2
frq_min = 1/wvl_max
frq_max = 1/wvl_min
frq_cen = 0.5*(frq_min+frq_max)
dfrq = frq_max-frq_min
nfrq = 100
## needs to be large enough to resolve skin depth of Au sphere
resolution = 80
dpml = 0.5*wvl_max
dair = 0.5*wvl_max
pml_layers = [mp.PML(thickness=dpml)]
s = 2*(dpml+dair+r)
cell_size = mp.Vector3(s,s,s)
# is_integrated=True necessary for any planewave source extending into PML
sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
center=mp.Vector3(-0.5*s+dpml),
size=mp.Vector3(0,s,s),
component=mp.Ez)]
## symmetries sometimes cause field instabilities at certain resolutions
symmetries = [mp.Mirror(mp.Y),
mp.Mirror(mp.Z,phase=-1)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
sources=sources,
k_point=mp.Vector3(),
symmetries=symmetries,
default_material=mp.Medium(index=n_H2O))
# padding of closed flux box around Au sphere
dpad = 0.5*dair
box_x1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=-r-dpad),size=mp.Vector3(0,2*(r+dpad),2*(r+dpad))))
box_x2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=+r+dpad),size=mp.Vector3(0,2*(r+dpad),2*(r+dpad))))
box_y1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=-r-dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_y2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=+r+dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_z1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=-r-dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))
box_z2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=+r+dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))
sim.run(until_after_sources=10)
freqs = mp.get_flux_freqs(box_x1)
box_x1_data = sim.get_flux_data(box_x1)
box_x2_data = sim.get_flux_data(box_x2)
box_y1_data = sim.get_flux_data(box_y1)
box_y2_data = sim.get_flux_data(box_y2)
box_z1_data = sim.get_flux_data(box_z1)
box_z2_data = sim.get_flux_data(box_z2)
box_x1_flux0 = mp.get_fluxes(box_x1)
sim.reset_meep()
geometry = [mp.Sphere(material=Au,
center=mp.Vector3(),
radius=r)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
sources=sources,
k_point=mp.Vector3(),
symmetries=symmetries,
default_material=mp.Medium(index=n_H2O),
eps_averaging=False,
geometry=geometry)
box_x1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=-r-dpad),size=mp.Vector3(0,2*(r+dpad),2*(r+dpad))))
box_x2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=+r+dpad),size=mp.Vector3(0,2*(r+dpad),2*(r+dpad))))
box_y1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=-r-dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_y2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=+r+dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_z1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=-r-dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))
box_z2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=+r+dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))
sim.load_minus_flux_data(box_x1, box_x1_data)
sim.load_minus_flux_data(box_x2, box_x2_data)
sim.load_minus_flux_data(box_y1, box_y1_data)
sim.load_minus_flux_data(box_y2, box_y2_data)
sim.load_minus_flux_data(box_z1, box_z1_data)
sim.load_minus_flux_data(box_z2, box_z2_data)
sim.run(until_after_sources=100)
box_x1_flux = mp.get_fluxes(box_x1)
box_x2_flux = mp.get_fluxes(box_x2)
box_y1_flux = mp.get_fluxes(box_y1)
box_y2_flux = mp.get_fluxes(box_y2)
box_z1_flux = mp.get_fluxes(box_z1)
box_z2_flux = mp.get_fluxes(box_z2)
scatt_flux = np.asarray(box_x1_flux)-np.asarray(box_x2_flux)+np.asarray(box_y1_flux)-np.asarray(box_y2_flux)+np.asarray(box_z1_flux)-np.asarray(box_z2_flux)
intensity = np.asarray(box_x1_flux0)/(2*r)**2
scatt_cross_section = np.divide(scatt_flux,intensity)
scatt_eff_meep = scatt_cross_section*-1/(np.pi*r**2)
trace = lambda t: (t[0][0]+t[1][1]+t[2][2])/3
scatt_eff_theory = [ps.MieQ(np.sqrt(trace(Au.epsilon(f))),1000/f,2*r*1000,nMedium=n_H2O,asDict=True)['Qsca'] for f in freqs]
if mp.am_master():
plt.figure(dpi=150)
plt.loglog(2*np.pi*r*np.asarray(freqs),scatt_eff_meep,'bo-',label='Meep')
plt.loglog(2*np.pi*r*np.asarray(freqs),scatt_eff_theory,'ro-',label='theory')
plt.grid(True,which="both",ls="-")
plt.xlabel('(sphere circumference)/wavelength, 2πr/λ')
plt.ylabel('scattering efficiency, σ/πr$^{2}$')
plt.legend(loc='upper left')
plt.title('Mie Scattering of a Au Sphere in H$_2$O')
plt.tight_layout()
plt.savefig("mie_scattering_Au_in_H2O.png")
np.savez('scatt_eff_res{}.npz'.format(resolution),scatt_eff_meep=scatt_eff_meep)
WeiKuoLi
commented
7 months ago Hello Sir thank you for your nice and detailed explanation. However I want to point out that I think in your script: intensity = np.asarray(box_x1_flux0)/(2*r)**2
should be intensity = np.asarray(box_x1_flux0)/(2*(r+dpad))**2
Since the box is now bigger, so is the area of box_x1.
Cheers
Hello!
I've been racking my brains for two weeks with a problem I hope you'll be able to help me with. Sorry to borrow some of your time and thanks for reading.
I'm simulating an Au sphere with 103 nm diameter submerged in water. To know whether I've understood Meep or not, I'm comparing the results with Mie theory. However, the difference between the model's and the simulation's maximum scattering wavelength is huge. The simulation presents an alarming redshift of 21 nm because its maximum is at 601 nm aprox. and the theory's maximum is at 582 nm aprox. Moreover, the shape of the spectrum seems to be wrong too, because it's wider than it should be when compared to the theory.
Regarding my Meep environment, I've recently updated it and it says it has Meep version 1.17.2-alpha.706753. I'm running the code in parallel in my computer using
mpirun --use-hwthread-cpus -np 6with MPI version 1.3. That way I get to use most of the Intel Core i7-9700F CPU of my group's PC, wich also has 48 GB of RAM and 48 GB of Swap Space in its Ubuntu 20.04.2 LTS.The code is very similar to the one that appears on the tutorial. the only unknown function to you should be
import_mediumfrom myv_materialsmodule, but all it does is rescaling the materials' properties according tofrom_um_factor(actually you'll see it's mostly a copy of Meep's Materials library). Since I'm calculating scattering cross section, I have to run two simulations like this:What could the problem be here? Is it OK to use water as
default_materialin both simulations?One more time, thank you for taking the time to read this! I hope you can help me.