smith_purcell_simulator.py

#!/usr/bin/env python
# coding: utf-8

# # About This Code

# ## Headers, Functions and Constants

# In[130]:


#Headers
import meep as mp
import numpy as np
import h5py
import yaml
import sys
from mpi4py import MPI
from scipy.io import loadmat
from scipy.io import savemat
import meep.materials as meepmat

#Constants
# -- Units noted.  
# -- Fundamental unit length of 1-um assumed for meep calculations (see meep manual)
tcon = 3.33564 #time conversion to meep units (fs/meep time)
q =  1.6e-19 # electron charge in coulomb
m =  9.10938356e-31 # electron mass in kg 

#Functions

#Functions for calculating index of refraction -- meep style
def calc_sig_d(n, k, fcen):
    eps = (n + 1j*k)**2
    eps_r = np.real(eps)
    return 2*np.pi*fcen*np.imag(eps)/eps_r

def calc_eps_r(n, k):
    eps = (n + 1j*k)**2
    return np.real(eps)


def run_simulation(save_prefix):

    #Open the settings yaml file
    with open(save_prefix + '_settings.yaml', 'r') as f:
        settings = yaml.safe_load(f)

    # ### Grating Settings
       
    #Grating Structure Properties
    grating_N = int(settings['grating_N']) #No. of grating elements
    grating_period = np.double(settings['grating_period']) #grating period in um
    grating_height = np.double(settings['grating_height']) #grating height in um
    grating_fill_factor = np.double(settings['grating_fill_factor']) #fill-factor from 0 to 1

    #If you want another grating under the underlayer...
    lower_grating_height = np.double(settings['lower_grating_height']) 

    #Grating Material Properties
    grating_material_eps = np.double(settings['grating_material_eps'])
    want_grating_material_by_name = settings['want_grating_material_by_name']
    grating_material_name = settings['grating_material_name']

    #Grating fill material properties
    # -- This is the material between each "tooth" of the grating
    grating_fill_material_eps = np.double(settings['grating_fill_material_eps'])
    want_grating_fill_material_by_name = settings['want_grating_fill_material_by_name']
    grating_fill_material_name = settings['grating_fill_material_name']

    #This underlayer is a layer of the grating type of whatever thickness. This is more for traditional
    #metalic gratings on-top of a planar film.  
    # -- Can set to 0 thickness if you do not want it
    underlayer_thick = np.double(settings['underlayer_thick'])
    want_underlayer_material_by_name = settings['want_underlayer_material_by_name']
    underlayer_material_eps = np.double(settings['underlayer_material_eps'])
    underlayer_material_name = settings['underlayer_material_name']


    #Substrate Material (semi-infinite film underneath grating):
    substrate_material_eps = np.double(settings['substrate_material_eps'])


    # ### Electron Settings

    e_height = np.double(settings['e_height']) #Height of electron above grating surface (um)
    W_e = np.double(settings['W_e']) # Kinetic energy of electron (eV)


    # ### FDTD-Specific Settings

    #Simulation Parameters
    # -- These parameters are important to running the simulation

    #Is this the reference simulation?
    reference = settings['reference']

    #PML Thicknesses
    dpml = np.double(settings['dpml']);

    #Monitor settings
    monitor_width = np.double(settings['monitor_width']) #Width in no. of gratings
    monitor_height = np.double(settings['monitor_height']) #height as a ratio of cell size

    # Sampling time for monitors in fs
    t_sample = np.double(settings['t_sample'])

    #EPS Averaging in Meep Calculation
    avgeps = np.double(settings['avgeps'])

    #Cell height (without PMLs)
    sy = np.double(settings['sy']) # size of y-dimension

    #Spatial Resolution
    res = np.double(settings['res']) # Resolution: points per meep-unit
    # e.g. if using 1 meep unit = 1 um,
    # 1/res um spacing

    # ### Frequency Output Range

    #Frequency range for calculations
    fcen  = np.double(settings['fcen']) #Central frequency (inverse wavelength in um)
    nfreq  = int(settings['nfreq']) #Number of frequency samples

    # ## Calculated Parameters
    # 
    # These are remaining parameters that are calculated based on user inputs.

    # ### Calculate Remaining Cell Sizes
    # 
    # We want to calculate the active cell width `sx` as it depends on the number of grating periods we include in our simulation and it is all well-defined.  We will also then calculate the total cell sizes `SX`, `SY`.  The z-sizes are all 0 for our 2D simulation here. 

    #Calculate the bandwidth for the monitors -- 1.9 fcen tends to work fine
    df = 1.9*fcen
    
    #Determine overall simulation cell size
    sx  = grating_period*grating_N*2.0 #Width of the cell in um

    #Full simulation cell sizes accounting for PMLs
    sX = 2.0*dpml+ sx
    sY = 2.0*dpml + sy
    sZ = 0.0

    # ### Electron Velocity
    v_e = 1e6*1e-15*tcon*np.sqrt(2.0*q*W_e/m)


    # ## Setup the FDTD Simulation

    # ### Define Materials

    # In[138]:


    if(want_grating_material_by_name):
        grating_material = eval('meepmat.' + grating_material_name)
    else:
        grating_material = mp.Medium(epsilon=grating_material_eps)
    
    if(want_grating_fill_material_by_name):
        grating_fill_material = eval('meepmat.' + grating_fill_material_name)
    else:
        grating_fill_material = mp.Medium(epsilon=grating_fill_material_eps)
    
    if(want_underlayer_material_by_name):
        underlayer_material = eval('meepmat.' + underlayer_material_name)
    else:
        underlayer_material = mp.Medium(epsilon=underlayer_material_eps)


    substrate_material = mp.Medium(epsilon=substrate_material_eps)


    # ### Define the Geometry

    # -- Cell --
    cell_size = mp.Vector3(sX, sY, sZ)

    # -- Geometry --
    geometry = []

    #Substrate Block
    geometry.append(mp.Block(center=mp.Vector3(0.0,
                                               float(-0.25*sY),
                                               0.0),
                             size=mp.Vector3(mp.inf,
                                             float(0.5*sY),
                                             mp.inf),
                             material=substrate_material))

    #Grating Fill Block
    geometry.append(mp.Block(center=mp.Vector3(0.0,
                                               float(0.5*grating_height),
                                               0.0),
                             size=mp.Vector3(mp.inf,
                                             float(grating_height),
                                             mp.inf),
                             material=grating_fill_material))

    #Underlayer Block
    geometry.append(mp.Block(center=mp.Vector3(0.0,
                                               float(-0.5*underlayer_thick),
                                               0.0),
                             size=mp.Vector3(mp.inf,
                                             float(underlayer_thick),
                                             mp.inf),
                             material=underlayer_material))


    if(not(reference)):

        #Loop through and create all desired objects...
        for i in range(grating_N):

            #Upper Grating
            geometry.append(mp.Block(center=mp.Vector3(float(grating_N*grating_period*-0.5 + grating_period*0.25 + i*grating_period),
                                                       float(grating_height*0.5),
                                                       0.0),
                                     size=mp.Vector3(float(grating_period*grating_fill_factor),
                                                     float(grating_height),
                                                     mp.inf),
                                     material=grating_material))
        
            #Lower Grating
            geometry.append(mp.Block(center=mp.Vector3(float(grating_N*grating_period*-0.5 + grating_period*0.25 + i*grating_period),
                                                       float(lower_grating_height*-0.5 - underlayer_thick),
                                                       0.0),
                                     size=mp.Vector3(float(grating_period*grating_fill_factor),
                                                     float(lower_grating_height),
                                                     mp.inf),
                                     material=grating_material))       


    # ### Boundary Conditions

    # -- Periodic Bdr. Cond -- 
    k_point = mp.Vector3(0, 0, 0)

    # ### Symmetries
    # 
    # This simulation does not have symmetry.  However, you should note that if your simulation has symmetry, then you can speed up the simulation significantly by including it. 
    # 
    # This symmetry is NOT symmetry of the objects, but of the fields.  If it is inversion symmetric, then you can include a negative pre-factor to account for it.  

    # ### Define the PML Layers

    # -- PML Layers --
    pml_layers = [mp.PML(thickness=np.double(dpml), direction=mp.X),
                  mp.PML(thickness=np.double(dpml), direction=mp.Y),
                  mp.PML(thickness=np.double(dpml), direction=mp.Z)]


    # ### Define the Field Monitors

    # In[142]:


    #Create a "volume" monitor (technically 2D as z = 0) for fields:
    monitor_xy = mp.Volume(center=mp.Vector3(0, 0),
                           size=mp.Vector3(monitor_width*grating_period, monitor_height*sy, 0))


    # ### Sources
    # 
    # Here will will have a delta source that is a point current source with zero size.  This reflects a moving charge in space.


    #Y-Polarized Gaussian source next to the dpml at -0.5 sz
    sources = [mp.Source(mp.ContinuousSource(frequency=1e-10),
                         component=mp.Ex,
                         center=mp.Vector3(sx*-0.5, grating_height+e_height),
                         size=mp.Vector3(0,0,0))]


    # ### Define our Simulation

    # In[144]:


    # -- Load Settings Into Simulation -- 
    sim = mp.Simulation(resolution=res,
                        cell_size=cell_size,
                        boundary_layers=pml_layers,
                        sources=sources,
                        k_point=k_point,
                        geometry=geometry,
                        dimensions=2,
                        filename_prefix=save_prefix)


    # ### Define and Add the Flux Regions
    # 
    # We define the flux regions and add them to our simulation.

    #Air Side
    air_fr = mp.FluxRegion(center=mp.Vector3(0, sy*0.05, 0), 
                           size=mp.Vector3(sx, 0, 0),
                           direction=mp.Y)
    air_flux = sim.add_flux(fcen, df, nfreq, air_fr)

    #Substrate Side
    substrate_fr = mp.FluxRegion(center=mp.Vector3(0, -1*sy*0.05, 0), 
                                 size=mp.Vector3(sx, 0, 0),
                                 direction=mp.Y)
    substrate_flux = sim.add_flux(fcen, df, nfreq, substrate_fr)


    # ## Run the Simulation
    # 
    # Now we have it all setup -- we just have to run the simulation.


    #Step function for moving the source with time given the electron velocity...
    # -- Moving delta current source is equivalent to moving charge in free-space
    def move_source(sim):
        sim.change_sources([mp.Source(mp.ContinuousSource(frequency=1e-10),
                                      component=mp.Ex,
                                      center=mp.Vector3(-0.5*sx + v_e*sim.meep_time(), 
                                                        grating_height + e_height,
                                                        0),
                                      size=mp.Vector3(0,0,0))])

    # -- Run Simulation -- 
    sim.run(move_source,
            mp.at_beginning(mp.output_epsilon),
            mp.in_volume(monitor_xy,
                         mp.to_appended("Ex_xy",
                                        mp.at_every(t_sample/tcon, mp.output_efield_x))),
            mp.in_volume(monitor_xy,
                         mp.to_appended("Ey_xy",
                                        mp.at_every(t_sample/tcon, mp.output_efield_y))),
            until = sx/v_e)


    # Get the flux outputs and save them


    air_flux_out = mp.get_fluxes(air_flux)
    substrate_flux_out = mp.get_fluxes(substrate_flux)
    flux_freqs = mp.get_flux_freqs(air_flux)

    savemat(save_prefix + '_output_flux.mat',
            {'flux_freqs':flux_freqs,
             'air_flux_out':air_flux_out,
             'substrate_flux_out':substrate_flux_out})

if __name__ == "__main__":
   run_simulation(sys.argv[1])