Laser-Ion Acceleration with a Planar Target
This example shows how to model laser-ion acceleration with planar targets of solid density [4, 5, 6]. The acceleration mechanism in this scenario depends on target parameters
Although laser-ion acceleration requires full 3D modeling for adequate description of the acceleration dynamics, especially the acceleration field lengths and decay times, this example models a 2D example. 2D modeling can often hint at a qualitative overview of the dynamics, but mostly saves computational costs since the plasma frequency (and Debye length) of the plasma determines the resolution need in laser-solid interaction modeling.
Note
TODO: The Python (PICMI) input file needs to be created.
Note
The resolution of this 2D case is extremely low by default. You will need a computing cluster for adequate resolution of the target density, see comments in the input file.
Warning
It is strongly advised to set the parameters <species>.zmin / zmax / xmin / ... when working with highly dense targets that are limited in one or multiple dimensions.
The particle creation routine will first create particles everywhere between these limits (defaulting to box size if unset), setting particles to invalid only afterwards based on the density profile.
Not setting these parameters can quickly lead to memory overflows.
Run
This example can be run either as:
Python script: (TODO) or
WarpX executable using an input file:
warpx.2d inputs_2d
For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.
Note
TODO: This input file should be created following the inputs_2d file.
#################################
# 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>"
Analyze
Note
This section is TODO.
Visualize
Note
This section is TODO.