In trying to educate myself about Meep and FDTD, I have been running various toy simulations, and most recently have tried to replicate metal reflectivity results as presented at
http://juluribk.com/2011/04/27/plasmonic-materials-in-meep/ ; the project file linked from that page is also the source of the LD model parameters I have used. I am using the github source of Dec 8th commits, on Arch Linux. According to the value of resolution I choose, the simulation may succeed or fail with a message such as: "meep: Lorentzian pole at too high a frequency 1e-20 for stability with dt = 0.00125313: reduce the Courant factor, increase the resolution, or use a different dielectric model" but the pattern of success/failure bewilders me. For example: resolution (set in run_parameters.py) 50 - succeeds 100 - fails 150 - succeeds 200 - succeeds 300 - fails 350 - succeeds 400 - fails 401 - succeeds 399 - fails 398 - succeeds When the simulation succeeds, the results closely match those in the source reference and are almost superimposable apart from the resolution=50 run. I think the simulation files are pretty simple - attached. Am I doing something fundamentally wrong - or is there a rational explanation, or bug? Many thanks, Ian
metal_mirror.sh
Description: application/shellscript
import meep as mp import metals from run_parameters import * #k_point=mp.Vector3(0,0,0) #sym = mp.Mirror(direction=mp.Y, phase=1) cell=mp.Vector3(x_size, y_size) sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=mp.Vector3(0, source_y), size=mp.Vector3(x_size, 0))] pml_layers=[mp.PML(pmld)] sim=mp.Simulation(cell_size=cell, boundary_layers=pml_layers, #geometry=geometry, sources=sources, #symmetries=[sym], resolution=resolution) refl_fr = mp.FluxRegion(center=mp.Vector3(0, fr_y), size=mp.Vector3(x_size/2-pmld, 0)) #pt=mp.Vector3(0, 0) refl = sim.add_flux(fcen, df, nfreq, refl_fr) sim.run(mp.at_beginning(mp.output_epsilon), mp.to_appended("ez", mp.at_every(0.1, mp.output_efield_z)), until=sim_duration) sim.save_flux('refl-flux', refl) sim.display_fluxes(refl)
import meep as mp import metals from run_parameters import * cell=mp.Vector3(x_size, y_size) sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=mp.Vector3(0, source_y), size=mp.Vector3(x_size, 0))] geometry=[mp.Block(mp.Vector3(x_size, metal_thickness, 1e20), center=mp.Vector3(0, metal_y), material=metals.silver)] pml_layers=[mp.PML(pmld)] sim=mp.Simulation(cell_size=cell, boundary_layers=pml_layers, geometry=geometry, sources=sources, resolution=resolution) refl_fr = mp.FluxRegion(center=mp.Vector3(0, fr_y), size=mp.Vector3(x_size/2-pmld, 0)) #pt=mp.Vector3(0, 0) refl = sim.add_flux(fcen, df, nfreq, refl_fr) sim.load_minus_flux('refl-flux', refl) sim.run(mp.at_beginning(mp.output_epsilon), mp.to_appended("ez", mp.at_every(0.1, mp.output_efield_z)), until=sim_duration) sim.display_fluxes(refl)
import meep as mp susc_aluminium = [ mp.LorentzianSusceptibility(frequency=8.06e-21, gamma=0.0379, sigma=1.173e42), mp.LorentzianSusceptibility(frequency=0.13066, gamma=0.2686, sigma=1940.97) ] susc_silver = [ mp.LorentzianSusceptibility(frequency=1e-20, gamma=0.038715, sigma=4.4625e41), mp.LorentzianSusceptibility(frequency=0.65815, gamma=3.1343, sigma=7.9247), mp.LorentzianSusceptibility(frequency=3.6142, gamma=0.36456, sigma=0.50133), mp.LorentzianSusceptibility(frequency=6.6017, gamma=0.052426, sigma=0.013329), mp.LorentzianSusceptibility(frequency=7.3259, gamma=0.7388, sigma=0.82655), mp.LorentzianSusceptibility(frequency=16.365, gamma=1.9511, sigma=1.1133) ] silver = mp.Medium(epsilon=1.0, E_susceptibilities=susc_silver)
x_size=.6 y_size=.5 pmld=.10 metal_thickness=.1 sim_duration=11 fcen = 3.333 df = 3.0 source_y=-y_size/2+1.5*pmld metal_y=y_size/2-metal_thickness/2-1.5*pmld fr_y=((y_size/2-1.5*pmld-metal_thickness)+(1.5*pmld-y_size/2))/2 nfreq=250 resolution=398 if __name__=="__main__": print("Simulation space from Y = ", -y_size/2, " to Y = ", y_size/2) print("Source position = ", source_y) print("Metal position = ", metal_y, " from ", metal_y-metal_thickness/2, " to ", metal_y+metal_thickness/2) print("Flux measured at = ", fr_y)
_______________________________________________ meep-discuss mailing list meep-discuss@ab-initio.mit.edu http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss