#################################
# Domain, Resolution & Numerics
#

# time-scale with highly kinetic dynamics
stop_time = 0.2e-12            # [s]
# time-scale for converged ion energy
#   notes: - effective acc. time depends on laser pulse
#          - ions will start to leave the box
#stop_time = 1.0e-12           # [s]

# quick tests at ultra-low res. (CI)
#amr.n_cell = 384 512

# proper resolution for 10 n_c excl. acc. length
# (>=1x V100)
amr.n_cell = 2688 3712

# proper resolution for 30 n_c (dx<=3.33nm) incl. acc. length
# (>=6x V100)
#amr.n_cell = 7488 14720

# simulation box, no MR
#   note: increase z (space & cells) for converging ion energy
amr.max_level = 0
geometry.dims = 2
geometry.prob_lo = -7.5e-6 -5.e-6
geometry.prob_hi =  7.5e-6 25.e-6

# Boundary condition
boundary.field_lo = pml pml
boundary.field_hi = pml pml

# Order of particle shape factors
algo.particle_shape = 3

# improved plasma stability for 2D with very low initial target temperature
# when using Esirkepov current deposition with energy-conserving field gather
interpolation.galerkin_scheme = 0

# numerical tuning
warpx.cfl = 0.999
warpx.use_filter = 1          # bilinear current/charge filter


#################################
# Performance Tuning
#
# simple tuning:
#   the numprocs product must be equal to the number of MPI ranks and splits
#   the domain on the coarsest level equally into grids;
#   slicing in the 2nd dimension is preferred for ideal performance
warpx.numprocs = 1 2   # 2 MPI ranks
#warpx.numprocs = 1 4  # 4 MPI ranks

# detail tuning instead of warpx.numprocs:
#   It is important to have enough cells in a block & grid, otherwise
#   performance will suffer.
#   Use larger values for GPUs, try to fill a GPU well with memory and place
#   few large grids on each device (you can go as low as 1 large grid / device
#   if you do not need load balancing).
#   Slicing in the 2nd dimension is preferred for ideal performance
#amr.blocking_factor = 64
#amr.max_grid_size_x = 2688
#amr.max_grid_size_y = 128  # this is confusingly named and means z in 2D

# load balancing
#   The grid & block parameters above are needed for load balancing:
#   an average of ~10 grids per MPI rank (and device) are a good granularity
#   to allow efficient load-balancing as the simulation evolves
algo.load_balance_intervals = 100
algo.load_balance_costs_update = Heuristic

# particle bin-sorting on GPU (ideal defaults not investigated in 2D)
#   Try larger values than the defaults below and report back! :)
#warpx.sort_intervals = 4    # default on CPU: -1 (off); on GPU: 4
#warpx.sort_bin_size = 1 1 1


#################################
# Target Profile
#

#   definitions for target extent and pre-plasma
my_constants.L    = 0.05e-6            # [m] scale length (>0)
my_constants.Lcut = 2.0e-6             # [m] hard cutoff from surface
my_constants.r0 = 2.5e-6               # [m] radius or half-thickness
my_constants.eps_z = 0.05e-6           # [m] small offset in z to make zmin, zmax interval larger than 2*(r0 + Lcut)
my_constants.zmax = r0 + Lcut + eps_z  # [m] upper limit in z for particle creation

particles.species_names = electrons hydrogen

# particle species
hydrogen.species_type = hydrogen
hydrogen.injection_style = NUniformPerCell
hydrogen.num_particles_per_cell_each_dim = 2 2
hydrogen.momentum_distribution_type = at_rest
# minimum and maximum z position between which particles are initialized
# --> should be set for dense targets limit memory consumption during initialization
hydrogen.zmin = -zmax
hydrogen.zmax = zmax
hydrogen.profile = parse_density_function
hydrogen.addRealAttributes = orig_x orig_z
hydrogen.attribute.orig_x(x,y,z,ux,uy,uz,t) = "x"
hydrogen.attribute.orig_z(x,y,z,ux,uy,uz,t) = "z"

electrons.species_type = electron
electrons.injection_style = NUniformPerCell
electrons.num_particles_per_cell_each_dim = 2 2
electrons.momentum_distribution_type = "gaussian"
electrons.ux_th = .01
electrons.uz_th = .01
# minimum and maximum z position between which particles are initialized
# --> should be set for dense targets limit memory consumption during initialization
electrons.zmin = -zmax
electrons.zmax = zmax

# ionization physics (field ionization/ADK)
#   [i1] none (fully pre-ionized):
electrons.profile = parse_density_function
#   [i2] field ionization (ADK):
#hydrogen.do_field_ionization = 1
#hydrogen.physical_element = H
#hydrogen.ionization_initial_level = 0
#hydrogen.ionization_product_species = electrons
#electrons.profile = constant
#electrons.density = 0.0

# collisional physics (binary MC model after Nanbu/Perez)
#collisions.collision_names = c_eH c_ee c_HH
#c_eH.species = electrons hydrogen
#c_ee.species = electrons electrons
#c_HH.species = hydrogen hydrogen
#c_eH.CoulombLog = 15.9
#c_ee.CoulombLog = 15.9
#c_HH.CoulombLog = 15.9

# number density: "fully ionized" electron density as reference
#   [material 1] cryogenic H2
my_constants.nc    = 1.742e27  # [m^-3]  1.11485e21 * 1.e6 / 0.8**2
my_constants.n0    = 30.0      # [n_c]
#   [material 2] liquid crystal
#my_constants.n0    = 192
#   [material 3] PMMA
#my_constants.n0    = 230
#   [material 4] Copper (ion density: 8.49e28/m^3; times ionization level)
#my_constants.n0    = 1400

# density profiles (target extent, pre-plasma and cutoffs defined above particle species list)

# [target 1] flat foil (thickness = 2*r0)
electrons.density_function(x,y,z) = "nc*n0*(
    if(abs(z)<=r0, 1.0, if(abs(z)<r0+Lcut, exp((-abs(z)+r0)/L), 0.0)) )"
hydrogen.density_function(x,y,z) = "nc*n0*(
    if(abs(z)<=r0, 1.0, if(abs(z)<r0+Lcut, exp((-abs(z)+r0)/L), 0.0)) )"

# [target 2] cylinder
#electrons.density_function(x,y,z) = "nc*n0*(
#    ((x*x+z*z)<=(r0*r0)) +
#    (sqrt(x*x+z*z)>r0)*(sqrt(x*x+z*z)<r0+Lcut)*exp( (-sqrt(x*x+z*z)+r0)/L ) )"
#hydrogen.density_function(x,y,z) = "nc*n0*(
#    ((x*x+z*z)<=(r0*r0)) +
#    (sqrt(x*x+z*z)>r0)*(sqrt(x*x+z*z)<r0+Lcut)*exp( (-sqrt(x*x+z*z)+r0)/L ) )"

# [target 3] sphere
#electrons.density_function(x,y,z) = "nc*n0*(
#    ((x*x+y*y+z*z)<=(r0*r0)) +
#    (sqrt(x*x+y*y+z*z)>r0)*(sqrt(x*x+y*y+z*z)<r0+Lcut)*exp( (-sqrt(x*x+y*y+z*z)+r0)/L ) )"
#hydrogen.density_function(x,y,z) = "nc*n0*(
#    ((x*x+y*y+z*z)<=(r0*r0)) +
#    (sqrt(x*x+y*y+z*z)>r0)*(sqrt(x*x+y*y+z*z)<r0+Lcut)*exp( (-sqrt(x*x+y*y+z*z)+r0)/L ) )"


#################################
# Laser Pulse Profile
#
lasers.names        = laser1
laser1.position     = 0. 0. -4.0e-6     # point the laser plane (antenna)
laser1.direction    = 0. 0. 1.          # the plane's (antenna's) normal direction
laser1.polarization = 1. 0. 0.          # the main polarization vector
laser1.a0           = 16.0              # maximum amplitude of the laser field [V/m]
laser1.wavelength   = 0.8e-6            # central wavelength of the laser pulse [m]
laser1.profile      = Gaussian
laser1.profile_waist = 4.e-6            # beam waist (E(w_0)=E_0/e) [m]
laser1.profile_duration = 30.e-15       # pulse length (E(tau)=E_0/e; tau=tau_E=FWHM_I/1.17741) [s]
laser1.profile_t_peak = 50.e-15         # time until peak intensity reached at the laser plane [s]
laser1.profile_focal_distance = 4.0e-6  # focal distance from the antenna [m]

# e_max = a0 * 3.211e12 / lambda_0[mu]
#   a0 = 16, lambda_0 = 0.8mu -> e_max = 64.22 TV/m


#################################
# Diagnostics
#
diagnostics.diags_names = diag1 openPMDfw openPMDbw

diag1.intervals = 100
diag1.diag_type = Full
diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_hydrogen
# reduce resolution of output fields
diag1.coarsening_ratio = 4 4
diag1.electrons.variables = w ux uy uz
diag1.hydrogen.variables = w ux uz orig_x orig_z
# demonstration of a spatial and momentum filter
diag1.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)
diag1.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)

openPMDfw.intervals = 100
openPMDfw.diag_type = Full
openPMDfw.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_hydrogen
# reduce resolution of output fields
openPMDfw.coarsening_ratio = 4 4
openPMDfw.electrons.variables = w ux uy uz
openPMDfw.hydrogen.variables = w ux uy uz orig_x orig_z
openPMDfw.format = openpmd
openPMDfw.openpmd_backend = h5
openPMDfw.species = electrons hydrogen
# demonstration of a spatial and momentum filter
openPMDfw.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)
openPMDfw.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)

openPMDbw.intervals = 100
openPMDbw.diag_type = Full
openPMDbw.fields_to_plot = rho_hydrogen
# reduce resolution of output fields
openPMDbw.coarsening_ratio = 4 4
openPMDbw.electrons.variables = w ux uy uz
openPMDbw.hydrogen.variables = w ux uy uz orig_x orig_z
openPMDbw.format = openpmd
openPMDbw.openpmd_backend = h5
openPMDbw.species = electrons hydrogen
# demonstration of a momentum filter
openPMDbw.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz<0)
openPMDbw.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz<0)


#################################
# Reduced Diagnostics
#

# histograms with 2.0 degree acceptance angle in fw direction
# 2 deg * pi / 180 : 0.03490658503 rad
# half-angle +/-   : 0.017453292515 rad
warpx.reduced_diags_names                   = histuH histue histuzAll FieldProbe_Z FieldProbe_ScatPoint FieldProbe_ScatLine LBC PhaseSpaceIons PhaseSpaceElectrons

histuH.type                                 = ParticleHistogram
histuH.intervals                            = 100
histuH.path                                 = "./"
histuH.species                              = hydrogen
histuH.bin_number                           = 1000
histuH.bin_min                              =  0.0
histuH.bin_max                              = 35.0
histuH.histogram_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)"
histuH.filter_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)"

histue.type                                 = ParticleHistogram
histue.intervals                            = 100
histue.path                                 = "./"
histue.species                              = electrons
histue.bin_number                           = 1000
histue.bin_min                              = 0.0
histue.bin_max                              = 0.1
histue.histogram_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)"
histue.filter_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)"

# just a test entry to make sure that the histogram filter is purely optional:
# this one just records uz of all hydrogen ions, independent of their pointing
histuzAll.type                                 = ParticleHistogram
histuzAll.intervals                            = 100
histuzAll.path                                 = "./"
histuzAll.species                              = hydrogen
histuzAll.bin_number                           = 1000
histuzAll.bin_min                              = -35.0
histuzAll.bin_max                              =  35.0
histuzAll.histogram_function(t,x,y,z,ux,uy,uz) = "uz"

FieldProbe_Z.type = FieldProbe
FieldProbe_Z.intervals = 100
FieldProbe_Z.integrate = 0
FieldProbe_Z.probe_geometry = Line
FieldProbe_Z.x_probe = 0.0
FieldProbe_Z.z_probe = -5.0e-6
FieldProbe_Z.x1_probe = 0.0
FieldProbe_Z.z1_probe = 25.0e-6
FieldProbe_Z.resolution = 3712

FieldProbe_ScatPoint.type = FieldProbe
FieldProbe_ScatPoint.intervals = 1
FieldProbe_ScatPoint.integrate = 0
FieldProbe_ScatPoint.probe_geometry = Point
FieldProbe_ScatPoint.x_probe = 0.0
FieldProbe_ScatPoint.z_probe = 15e-6

FieldProbe_ScatLine.type = FieldProbe
FieldProbe_ScatLine.intervals = 100
FieldProbe_ScatLine.integrate = 1
FieldProbe_ScatLine.probe_geometry = Line
FieldProbe_ScatLine.x_probe = -2.5e-6
FieldProbe_ScatLine.z_probe = 15e-6
FieldProbe_ScatLine.x1_probe = 2.5e-6
FieldProbe_ScatLine.z1_probe = 15e-6
FieldProbe_ScatLine.resolution = 201

# check computational load per box
LBC.type = LoadBalanceCosts
LBC.intervals = 100

PhaseSpaceIons.type                                 = ParticleHistogram2D
PhaseSpaceIons.intervals                            = 100
PhaseSpaceIons.species                              = hydrogen
PhaseSpaceIons.bin_number_abs                       = 1000
PhaseSpaceIons.bin_number_ord                       = 1000
PhaseSpaceIons.bin_min_abs                          = -5.e-6
PhaseSpaceIons.bin_max_abs                          = 25.e-6
PhaseSpaceIons.bin_min_ord                          = -20.e-6
PhaseSpaceIons.bin_max_ord                          = 20.e-6
PhaseSpaceIons.histogram_function_abs(t,x,y,z,ux,uy,uz,w) = "z"
PhaseSpaceIons.histogram_function_ord(t,x,y,z,ux,uy,uz,w) = "uz"
PhaseSpaceIons.value_function(t,x,y,z,ux,uy,uz,w) = "w"
# PhaseSpaceIons.filter_function(t,x,y,z,ux,uy,uz,w) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)"

PhaseSpaceElectrons.type                                 = ParticleHistogram2D
PhaseSpaceElectrons.intervals                            = 100
PhaseSpaceElectrons.species                              = electrons
PhaseSpaceElectrons.bin_number_abs                       = 1000
PhaseSpaceElectrons.bin_number_ord                       = 1000
PhaseSpaceElectrons.bin_min_abs                          = -5.e-6
PhaseSpaceElectrons.bin_max_abs                          = 25.e-6
PhaseSpaceElectrons.bin_min_ord                          = -20
PhaseSpaceElectrons.bin_max_ord                          = 20
PhaseSpaceElectrons.histogram_function_abs(t,x,y,z,ux,uy,uz,w) = "z"
PhaseSpaceElectrons.histogram_function_ord(t,x,y,z,ux,uy,uz,w) = "uz"
PhaseSpaceElectrons.value_function(t,x,y,z,ux,uy,uz,w) = "w"
PhaseSpaceElectrons.filter_function(t,x,y,z,ux,uy,uz,w) = "sqrt(x*x+y*y) < 1e-6"

#################################
# Physical Background
#
# This example is modeled after a target similar to the hydrogen jet here:
#   [1] https://doi.org/10.1038/s41598-017-10589-3
#   [2] https://arxiv.org/abs/1903.06428
#
authors = "Axel Huebl <axelhuebl@lbl.gov>"