# ----------------------------------------------------------------------------------------
# SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI
# ----------------------------------------------------------------------------------------
from math import pi, sqrt, exp
from numpy import s_
# SIMULATION PARAMTERS MATCHED TO 2019 ZVP QED paper simulation
# VERSION 2.0
# adjusted gaussian time profile
# adjusted electron momentum phase space for better calculation of acceleration
# VERSION 3.0
# plasma width now one wavelength long
# halved simulation size
# corrected relative size of plasma mirror and beam waist
# added Ex diagnostic
# electron vx diag removed log scaling
# removed px diag
# VERSION 4.0
# - shortened sim time since not interested in when filamentation times
# - matched beam frequency to ELI
# - fixed incorrect diagnostic due to sim size variation
# - decreased sim size again (less space before plasma - moved everything accordingly)
# - Ex->Ey in diagnostic
# - changed vx range to include negative direction
# VERSION 5.0
# - Shortened the temporal gaussian of the beam (to match Aimee ross)
# - reduced sim in x direction and accordingly moved the screen
# - changed range of KE bins
# - increased resolution of Ey so can see what is going on here, probs reduce/remove in future sims
# VERSION 6.0
# Matched to Aimee Ross input deck for Savin QED sims
# - set temp electrons to 300K (at this T collisions can be ignored but also not physical??)
# - made ions into Aluminium
# - DID NOT EDIT SPATIAL PROFILE OF BEAM - in aimee sims, beam only hits small part of surface not the whole thing but I am interested in edge effects see variation for this
# - changed beam profile again so peak amplitude on the fourth wave (take care: this pulse length is much shorter than feasible in experiment)
# VERSION 7.0
# - Added Ex diag to see pseudo capacitor
# - Consider shortening sim since don't care about when surface brakes up
# VERAION 7.beam
# Diags to split beam and plasma electrons (possible since beam electrons close to speed of light
# - resy = 128
# - added breight wheeler block
l0 = 2. * pi # laser wavelength [in code units]
t0 = l0 # optical cycle
Lsim = [4.*l0, 10.*l0] # length of the simulation
Tsim = 10.*t0 # duration of the simulation
resx = 128 # longitudinal cells in one laser wavelength (had to adjust so that patches works well in smilei - was 100)
resy = 128 # transverse cells in one laser wavelength (was 25 in QED 2019 sim)
rest = resx*11 # nb of timesteps in one optical cycle
Lcell = [l0/resx,l0/resy] # length of a cell
ne = 50. # electron density in simulation units (reference density = critical density?)
# Te = 115 / 511e3
Te = 0.0000000505 # electron temperature in simulation units COULD NOT FIND WHAT PARAMETER WAS USED FOR THIS IN THEIR SIMS
omega0 = 1.78*10**15 # reference frequency in SI units
time_peak_beam = 4*l0 # = 25.1 in S units and 14.11 in SI units
# We want 4x4 electron and 2x2 ion per cell
nppc_ele = 22*22 # number of ele per cell (had to change these slightly too
nppc_ion = 5*5 # number of ion per cell
cells_x = int(Lsim[0]/Lcell[0]) # cells in x direction
cells_y = int(Lsim[1]/Lcell[1]) # cells in y direction
ndump = 50 # number of timesteps between writing diagnostics
# density function for a flat plasma mirror
def density(x,y):
x_wall = 2*l0
if 2*l0 <= y <= 8*l0:
if x - x_wall <= 0:
return ne*exp((x - x_wall)/(0.2*l0))
elif x - x_wall <= l0:
return ne
return 0
def density_i(x,y):
return density(x,y)/13
Main(
geometry = "2Dcartesian",
interpolation_order = 2 ,
cell_length = Lcell,
grid_length = Lsim,
number_of_patches = [ 64, 32 ], # must be a power of 2, total number must be >= to number of MPI processes and preferably more than number of mp theads
timestep = t0/rest,
simulation_time = Tsim,
EM_boundary_conditions = [
['silver-muller'],
['silver-muller'],
],
reference_angular_frequency_SI = omega0,
print_every = ndump,
random_seed = smilei_mpi_rank
)
Species(
name = 'electron',
position_initialization = 'regular', # This should usually be 'random' or 'regular'
momentum_initialization = 'maxwell-juettner',
particles_per_cell = nppc_ele,
mass = 1.0,
charge = -1.0,
number_density = density,
temperature = [ Te ], # temperature in units of m_ec^2
radiation_model = "Monte-Carlo",
radiation_photon_species = "photon",
boundary_conditions = [
['remove'],
]
)
# Aluminium ions
Species(
name = 'ion',
position_initialization = 'regular',
momentum_initialization = 'cold',
particles_per_cell = nppc_ion,
mass = 1836.*27,
charge = 13.0,
number_density = density_i,
temperature = [ 0 ],
boundary_conditions = [
['remove'],
]
)
Species(
name = 'positron',
position_initialization = 'regular',
momentum_initialization = 'cold',
particles_per_cell = nppc_ele,
mass = 1.0,
charge = 1.0,
number_density = 0,
temperature = [ 0 ],
radiation_model = "Monte-Carlo",
radiation_photon_species = "photon",
boundary_conditions = [
['remove'],
]
)
Species(
name = 'photon',
position_initialization = 'regular',
momentum_initialization = 'cold',
particles_per_cell = nppc_ion,
mass = 0.,
charge = 0.,
number_density = 0,
temperature = [ 0 ],
boundary_conditions = [
['remove'],
],
multiphoton_Breit_Wheeler = ['electron','positron']
)
RadiationReaction(
)
MultiphotonBreitWheeler(
)
LaserGaussian2D(
a0 = 100., # normalized amplitude
omega = 1.,
focus = [2*l0, 5*l0], # coordinates of laser focus
waist = 6.*l0,
incidence_angle = 0.,
time_envelope = tgaussian(duration=2*time_peak_beam,fwhm=2*time_peak_beam/3,center=time_peak_beam)
)
# DIAGNOSTICS
# The purpose of this sim is to output the fast electrons and the HHG obtained so I require
# - electron density
# - electron energy
# - electron momentum
# - em fields
# - em energy
# When outputting the fields and electrons, will be useful to use a probe to get at surface of new object and a screen diagnostic for electrons
DiagFields(
every = 50,
fields = ["Ex", "Ey"],
subgrid = s_[::2, ::2]
)
def beam(particles):
weights = particles.weight
px = particles.px
weights[px>1]=0
return weights
def plasma(particles):
weights = particles.weight
px = particles.px
weights[px<=1]=0
return weights
DiagParticleBinning(
name = "beam_density",
deposited_quantity = beam,
every = ndump,
time_average = 1,
species = ["electron"],
axes = [
["x", 0., Lsim[0], 2000],
["y", 0., Lsim[1], 2000],
]
)
DiagParticleBinning(
name = "plasma_density",
deposited_quantity = plasma,
every = ndump,
time_average = 1,
species = ["electron"],
axes = [
["x", 0., Lsim[0], 2000],
["y", 0., Lsim[1], 2000],
]
)