diff options
29 files changed, 1831 insertions, 213 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index d411cc0a1..c32f1637f 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -336,6 +336,23 @@ Particle initialization * ``warpx.serialize_ics`` (`0 or 1`) Whether or not to use OpenMP threading for particle initialization. +* ``<species>.do_field_ionization`` (`0` or `1`) optional (default `0`) + Do field ionization for this species (using the ADK theory). + +* ``<species>.physical_element`` (`string`) + Only read if `do_field_ionization = 1`. Symbol of chemical element for + this species. Example: for Helium, use ``physical_element = He``. + +* ``<species>.ionization_product_species`` (`string`) + Only read if `do_field_ionization = 1`. Name of species in which ionized + electrons are stored. This species must be created as a regular species + in the input file (in particular, it must be in `particles.species_names`). + +* ``<species>.ionization_initial_level`` (`int`) optional (default `0`) + Only read if `do_field_ionization = 1`. Initial ionization level of the + species (must be smaller than the atomic number of chemical element given + in `physical_element`). + Laser initialization -------------------- diff --git a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py new file mode 100755 index 000000000..497a30097 --- /dev/null +++ b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py @@ -0,0 +1,53 @@ +#! /usr/bin/env python + +''' +Analysis script of a WarpX simulation of rigid injection in a boosted frame. + +A Gaussian electron beam starts from -5 microns, propagates rigidly up to +20 microns after which it expands due to emittance only (the focal position is +20 microns). The beam width is measured after ~50 microns, and compared with +the theory (with a 5% error allowed). + +The simulation runs in a boosted frame, and the analysis is done in the lab +frame, i.e., on the back-transformed diagnostics. +''' + +import sys, os, yt, glob +import numpy as np +import scipy.constants as scc +import read_raw_data +yt.funcs.mylog.setLevel(0) + +# filename = sys.argv[1] + +def get_particle_field(snapshot, species, field): + fn = snapshot + '/' + species + files = glob.glob(os.path.join(fn, field + '_*')) + files.sort() + all_data = np.array([]) + for f in files: + data = np.fromfile(f) + all_data = np.concatenate((all_data, data)) + return all_data + +# Read data from back-transformed diagnostics +snapshot = './lab_frame_data/snapshot00001' +header = './lab_frame_data/Header' +allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) +z = np.mean( get_particle_field(snapshot, 'beam', 'z') ) +w = np.std ( get_particle_field(snapshot, 'beam', 'x') ) + +# initial parameters +z0 = 20.e-6 +w0 = 1.e-6 +theta0 = np.arcsin(0.1) + +# Theoretical beam width after propagation if rigid ON +wth = np.sqrt( w0**2 + (z-z0)**2*theta0**2 ) +error = np.abs((w-wth)/wth) + +# Print error and assert small error +print("Beam position: " + str(z)) +print("Beam width : " + str(w)) +print("error: " + str(error)) +assert( error < 0.03 ) diff --git a/Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py b/Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py new file mode 100755 index 000000000..86214ad72 --- /dev/null +++ b/Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py @@ -0,0 +1,70 @@ +#! /usr/bin/env python + +''' +Analysis script of a WarpX simulation of rigid injection. + +A Gaussian electron beam starts from -5 microns, propagates rigidly up to +20 microns after which it expands due to emittance only (the focal position is +20 microns). The beam width is measured after ~50 microns, and compared with +the theory (with a 5% error allowed). + +As a help to the user, the script also compares beam width to the theory in +case rigid injection is OFF (i.e., the beam starts expanding from -5 microns), +in which case a warning is raised. +''' + +import sys +import yt +import numpy as np +import scipy.constants as scc +yt.funcs.mylog.setLevel(0) + +filename = sys.argv[1] + +# WarpX headers include more data when rigid injection is used, +# which gives an error with the last yt release. +# To avoid this issue, the three last lines of WarpXHeader are removed if +# needed. +def remove_rigid_lines(plotfile, nlines_if_rigid): + header_name = plotfile + '/WarpXHeader' + f = open(header_name, 'r') + file_lines = f.readlines() + nlines = len(file_lines) + f.close() + if nlines == nlines_if_rigid: + f = open(header_name, 'w') + f.writelines(file_lines[:-3]) + f.close() + +# Remove rigid injection header lines +remove_rigid_lines(filename, 18) +# Read beam parameters +ds = yt.load( filename ) +ad = ds.all_data() +# Beam longitudinal position +z = np.mean(ad['beam', 'particle_position_y'].v) +# Beam width +w = np.std(ad['beam', 'particle_position_x'].v) + +# initial parameters +z0 = 20.e-6 +z0_no_rigid = -5.e-6 +w0 = 1.e-6 +theta0 = np.arcsin(0.1) + +# Theoretical beam width after propagation if rigid OFF +# Inform the user if rigid injection simply off (just to be kind) +wth_no_rigid = np.sqrt( w0**2 + (z-z0_no_rigid)**2*theta0**2 ) +error_no_rigid = np.abs((w-wth_no_rigid)/wth_no_rigid) +if ( error_no_rigid < 0.05): + print("error no rigid: " + str(error_no_rigid)) + print("Looks like the beam defocuses as if rigid injection were OFF") + +# Theoretical beam width after propagation if rigid ON +wth = np.sqrt( w0**2 + (z-z0)**2*theta0**2 ) +error = np.abs((w-wth)/wth) +# Print error and assert small error +print("Beam position: " + str(z)) +print("Beam width : " + str(w)) +print("error: " + str(error)) +assert( error < 0.05 ) diff --git a/Examples/Modules/RigidInjection/inputs b/Examples/Modules/RigidInjection/inputs deleted file mode 100644 index 3c0d52ed0..000000000 --- a/Examples/Modules/RigidInjection/inputs +++ /dev/null @@ -1,88 +0,0 @@ -# Maximum number of time steps -max_step = 60 - -# number of grid points -amr.n_cell = 32 32 32 - -# Maximum allowable size of each subdomain in the problem domain; -# this is used to decompose the domain for parallel calculations. -amr.max_grid_size = 16 - -# Maximum level in hierarchy (for now must be 0, i.e., one level in total) -amr.max_level = 0 - -amr.plot_int = 1 # How often to write plotfiles. "<= 0" means no plotfiles. - -# Geometry -geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 0 # Is periodic? -geometry.prob_lo = -2. -2. -4. # physical domain -geometry.prob_hi = 2. 2. 4. - -# Verbosity -warpx.verbose = 1 - -# Algorithms - -# interpolation -interpolation.nox = 3 -interpolation.noy = 3 -interpolation.noz = 3 - -# CFL -warpx.cfl = 1.0 - -# Information about the particle species -particles.nspecies = 1 -particles.species_names = electrons -particles.rigid_injected_species = electrons - - - -# -# The electron species information -# - -electrons.charge = -q_e -electrons.mass = m_e -electrons.injection_style = "gaussian_beam" -electrons.x_rms = 0.1 -electrons.y_rms = 0.1 -electrons.z_rms = 0.1 -electrons.x_m = 0. -electrons.y_m = 0. -electrons.z_m = -1.5 -electrons.npart = 1000 -electrons.q_tot = -8.010883097437485e-07 - -electrons.profile = "constant" -electrons.density = 1 -electrons.momentum_distribution_type = "gaussian" -electrons.ux_m = 0. -electrons.uy_m = 0. -electrons.uz_m = 3. -electrons.ux_th = 0.01 -electrons.uy_th = 0.01 -electrons.uz_th = 0.01 - -electrons.xmin = -2 -electrons.xmax = 2 -electrons.ymin = -2 -electrons.ymax = 2 -electrons.zmin = -2 -electrons.zmax = 2 - -electrons.zinject_plane = 0. -electrons.projected = true -electrons.focused = false - -warpx.do_pml = 0 - -# Moving window -warpx.do_moving_window = 0 -warpx.moving_window_dir = z -warpx.moving_window_v = 1.0 # in units of the speed of light - -# Boosted frame -warpx.gamma_boost = 1.5 -warpx.boost_direction = z diff --git a/Examples/Modules/RigidInjection/inputs.BoostedFrame b/Examples/Modules/RigidInjection/inputs.BoostedFrame new file mode 100644 index 000000000..c7a60f14f --- /dev/null +++ b/Examples/Modules/RigidInjection/inputs.BoostedFrame @@ -0,0 +1,49 @@ +# stop_time = 1.5e-13 + +warpx.zmax_plasma_to_compute_max_step = 50.e-6 +warpx.gamma_boost = 5. +warpx.boost_direction = z +warpx.do_boosted_frame_diagnostic = 1 +warpx.num_snapshots_lab = 2 +warpx.dt_snapshots_lab = 1.8679589331096515e-13 + +amr.n_cell = 256 512 # 32 64 +amr.max_grid_size = 512 +amr.blocking_factor = 16 +amr.max_level = 0 +amr.plot_int = 100000 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 1 0 # Is periodic? +geometry.prob_lo = -50.e-6 -10.e-6 +geometry.prob_hi = 50.e-6 0.e-6 +warpx.cfl = .999 +warpx.do_pml = 1 +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 # in units of the speed of light +warpx.serialize_ics = 1 +particles.nspecies = 1 +particles.species_names = beam +particles.rigid_injected_species = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.x_rms = 1.e-6 +beam.y_rms = 1.e-6 +beam.z_rms = .5e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = -5.e-6 +beam.npart = 10000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 1000. +beam.ux_th = 100. +beam.uy_th = 100. +beam.uz_th = 0. +beam.zinject_plane = 20.e-6 +beam.rigid_advance = true +beam.projected = true +beam.focused = false diff --git a/Examples/Modules/RigidInjection/inputs.LabFrame b/Examples/Modules/RigidInjection/inputs.LabFrame new file mode 100644 index 000000000..730b46e01 --- /dev/null +++ b/Examples/Modules/RigidInjection/inputs.LabFrame @@ -0,0 +1,41 @@ +stop_time = 1.5e-13 +amr.n_cell = 32 64 +amr.max_grid_size = 256 +amr.blocking_factor = 16 +amr.max_level = 0 +amr.plot_int = 10000 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 1 0 # Is periodic? +geometry.prob_lo = -50.e-6 -10.e-6 +geometry.prob_hi = 50.e-6 0.e-6 +warpx.cfl = .999 +warpx.do_pml = 1 +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 # in units of the speed of light +warpx.serialize_ics = 1 +particles.nspecies = 1 +particles.species_names = beam +particles.rigid_injected_species = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.x_rms = 1.e-6 +beam.y_rms = 1.e-6 +beam.z_rms = .5e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = -5.e-6 +beam.npart = 2000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 1000. +beam.ux_th = 100. +beam.uy_th = 100. +beam.uz_th = 0. +beam.zinject_plane = 20.e-6 +beam.rigid_advance = true +beam.projected = true +beam.focused = false diff --git a/Examples/Modules/ionization/inputs.bf.rt b/Examples/Modules/ionization/inputs.bf.rt new file mode 100644 index 000000000..2fcb2fedc --- /dev/null +++ b/Examples/Modules/ionization/inputs.bf.rt @@ -0,0 +1,61 @@ +max_step = 420 +amr.n_cell = 16 800 +amr.max_grid_size = 64 +amr.blocking_factor = 16 +amr.plot_int = 10000 +geometry.coord_sys = 0 +geometry.is_periodic = 1 0 +geometry.prob_lo = -5.e-6 -40.e-6 +geometry.prob_hi = 5.e-6 0.e-6 +amr.max_level = 0 + +algo.maxwell_fdtd_solver = ckc +warpx.do_pml = 1 +warpx.cfl = .999 +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 +warpx.gamma_boost = 2. +warpx.boost_direction = z + +particles.nspecies = 2 +particles.species_names = electrons ions + +ions.mass = 2.3428415e-26 +ions.charge = q_e +ions.injection_style = nuniformpercell +ions.num_particles_per_cell_each_dim = 2 2 +ions.zmin = 0. +ions.zmax = 50.e-6 +ions.profile = constant +ions.density = 1. +ions.momentum_distribution_type = constant +ions.do_field_ionization = 1 +ions.ionization_initial_level = 2 +ions.ionization_product_species = electrons +ions.physical_element = N +ions.do_continuous_injection=1 + +electrons.mass = m_e +electrons.charge = -q_e +electrons.injection_style = nuniformpercell +electrons.num_particles_per_cell_each_dim = 2 2 +electrons.zmin = 0. +electrons.zmax = 50.e-6 +electrons.profile = constant +electrons.density = 2. +electrons.momentum_distribution_type = constant +electrons.do_continuous_injection=1 + +lasers.nlasers = 1 +lasers.names = laser1 +laser1.profile = Gaussian +laser1.position = 0. 0. -1.e-6 +laser1.direction = 0. 0. 1. +laser1.polarization = 1. 0. 0. +laser1.e_max = 7.224e12 # a0=1.8 +laser1.profile_waist = 1.e10 +laser1.profile_duration = 26.685e-15 +laser1.profile_t_peak = 60.e-15 +laser1.profile_focal_distance = 0 +laser1.wavelength = 0.8e-6 diff --git a/Examples/Modules/ionization/inputs.rt b/Examples/Modules/ionization/inputs.rt new file mode 100644 index 000000000..e3e4622a6 --- /dev/null +++ b/Examples/Modules/ionization/inputs.rt @@ -0,0 +1,54 @@ +max_step = 1600 +amr.n_cell = 16 800 +amr.max_grid_size = 64 +amr.blocking_factor = 16 +amr.plot_int = 10000 +geometry.coord_sys = 0 +geometry.is_periodic = 1 0 +geometry.prob_lo = -5.e-6 0.e-6 +geometry.prob_hi = 5.e-6 20.e-6 +amr.max_level = 0 + +algo.maxwell_fdtd_solver = ckc +warpx.do_pml = 1 +warpx.cfl = .999 + +particles.nspecies = 2 +particles.species_names = electrons ions + +ions.mass = 2.3428415e-26 +ions.charge = q_e +ions.injection_style = nuniformpercell +ions.num_particles_per_cell_each_dim = 2 1 +ions.zmin = 5.e-6 +ions.zmax = 15.e-6 +ions.profile = constant +ions.density = 1. +ions.momentum_distribution_type = constant +ions.do_field_ionization = 1 +ions.ionization_initial_level = 2 +ions.ionization_product_species = electrons +ions.physical_element = N + +electrons.mass = m_e +electrons.charge = -q_e +electrons.injection_style = nuniformpercell +electrons.num_particles_per_cell_each_dim = 1 1 +electrons.zmin = 1. +electrons.zmax = 15.e-6 +electrons.profile = constant +electrons.density = 2. +electrons.momentum_distribution_type = constant + +lasers.nlasers = 1 +lasers.names = laser1 +laser1.profile = Gaussian +laser1.position = 0. 0. 3.e-6 +laser1.direction = 0. 0. 1. +laser1.polarization = 1. 0. 0. +laser1.e_max = 7.224e12 # a0=1.8 +laser1.profile_waist = 1.e10 +laser1.profile_duration = 26.685e-15 +laser1.profile_t_peak = 60.e-15 +laser1.profile_focal_distance = 0 +laser1.wavelength = 0.8e-6 diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py new file mode 100755 index 000000000..b94541f90 --- /dev/null +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -0,0 +1,77 @@ +#! /usr/bin/env python + +""" +This script tests the result of the ionization module in WarpX. + +Input files inputs.rt and inputs.bf.rt are used to reproduce the test from +Chen, JCP, 2013, figure 2 (in the lab frame and in a boosted frame, +respectively): a plane-wave laser pulse propagates through a +uniform N2+ neutral plasma and further ionizes the Nitrogen atoms. This test +checks that, after the laser went through the plasma, ~32% of Nitrogen +ions are N5+, in agreement with theory from Chen's article. +""" + +import sys +import yt +import numpy as np +import scipy.constants as scc +yt.funcs.mylog.setLevel(0) + +# Open plotfile specified in command line, and get ion's ionization level. +filename = sys.argv[1] +ds = yt.load( filename ) +ad = ds.all_data() +ilev = ad['ions', 'particle_ionization_level'].v + +# Fraction of Nitrogen ions that are N5+. +N5_fraction = ilev[ilev == 5].size/float(ilev.size) + +print("Number of ions: " + str(ilev.size)) +print("Number of N5+ : " + str(ilev[ilev == 5].size)) +print("N5_fraction : " + str(N5_fraction)) + +do_plot = False +if do_plot: + import matplotlib.pyplot as plt + all_data_level_0 = ds.covering_grid(level=0,left_edge=ds.domain_left_edge, + dims=ds.domain_dimensions) + F = all_data_level_0['boxlib', 'Ex'].v.squeeze() + extent = [ ds.domain_left_edge[1], ds.domain_right_edge[1], + ds.domain_left_edge[0], ds.domain_right_edge[0] ] + ad = ds.all_data() + + # Plot ions with ionization levels + species = 'ions'; + xi = ad[species, 'particle_position_x'].v + zi = ad[species, 'particle_position_y'].v + ii = ad[species, 'particle_ionization_level'].v + plt.figure(figsize=(10,10)) + plt.subplot(211) + plt.imshow(np.abs(F), extent=extent, aspect='auto', + cmap='magma', origin='default') + plt.colorbar() + for lev in range(int(np.max(ii)+1)): + select = (ii == lev) + plt.scatter(zi[select],xi[select],s=.2, + label='ionization level: ' + str(lev)) + plt.legend() + plt.title("abs(Ex) (V/m) and ions") + plt.xlabel("z (m)") + plt.ylabel("x (m)") + plt.subplot(212) + plt.imshow(np.abs(F), extent=extent, aspect='auto', + cmap='magma', origin='default') + plt.colorbar() + + # Plot electrons + species = 'electrons'; + if species in [x[0] for x in ds.field_list]: + xe = ad[species, 'particle_position_x'].v + ze = ad[species, 'particle_position_y'].v + plt.scatter(ze,xe,s=.1,c='r',label='electrons') + plt.title("abs(Ex) (V/m) and electrons") + plt.xlabel("z (m)") + plt.ylabel("x (m)") + plt.savefig("image_ionization.pdf", bbox_inches='tight') + +assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 0fcb88e05..bc3239ebf 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -102,6 +102,39 @@ analysisRoutine = Examples/Tests/PML/analysis_pml_ckc.py #doVis = 0 #analysisRoutine = Examples/Tests/PML/analysis_pml_psatd.py +[RigidInjection_lab] +buildDir = . +inputFile = Examples/Modules/RigidInjection/inputs.LabFrame +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py + +[RigidInjection_boost_backtransformed] +buildDir = . +inputFile = Examples/Modules/RigidInjection/inputs.BoostedFrame +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +outputFile = "lab_frame_data" +doComparison = 0 +aux1File = Tools/read_raw_data.py +analysisRoutine = Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py + [nci_corrector] buildDir = . inputFile = Examples/Modules/nci_corrector/inputs2d @@ -117,6 +150,34 @@ compileTest = 0 doVis = 0 analysisRoutine = Examples/Modules/nci_corrector/ncicorr_analysis.py +[ionization_lab] +buildDir = . +inputFile = Examples/Modules/ionization/inputs.rt +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +analysisRoutine = Examples/Modules/ionization/ionization_analysis.py + +[ionization_boost] +buildDir = . +inputFile = Examples/Modules/ionization/inputs.bf.rt +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +analysisRoutine = Examples/Modules/ionization/ionization_analysis.py + [bilinear_filter] buildDir = . inputFile = Examples/Tests/SingleParticle/inputs @@ -294,23 +355,23 @@ particleTypes = electrons positrons analysisRoutine = Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py analysisOutputImage = langmuir_multi_2d_analysis.png -[Langmuir_multi_2d_psatd] -buildDir = . -inputFile = Examples/Tests/Langmuir/inputs.multi.2d.rt -dim = 2 -addToCompileString = USE_PSATD=TRUE -restartTest = 0 -useMPI = 1 -numprocs = 2 -useOMP = 1 -numthreads = 1 -compileTest = 0 -doVis = 0 -compareParticles = 1 -runtime_params = psatd.fftw_plan_measure=0 -particleTypes = electrons positrons -analysisRoutine = Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py -analysisOutputImage = langmuir_multi_2d_analysis.png +# [Langmuir_multi_2d_psatd] +# buildDir = . +# inputFile = Examples/Tests/Langmuir/inputs.multi.2d.rt +# dim = 2 +# addToCompileString = USE_PSATD=TRUE +# restartTest = 0 +# useMPI = 1 +# numprocs = 2 +# useOMP = 1 +# numthreads = 1 +# compileTest = 0 +# doVis = 0 +# compareParticles = 1 +# runtime_params = psatd.fftw_plan_measure=0 +# particleTypes = electrons positrons +# analysisRoutine = Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py +# analysisOutputImage = langmuir_multi_2d_analysis.png # [Langmuir_multi_2d_psatd_nodal] # buildDir = . diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 4ac5cb6a1..8bfa45a59 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -76,14 +76,15 @@ MultiParticleContainer::Checkpoint (const std::string& dir) const void MultiParticleContainer::WritePlotFile (const std::string& dir) const { - Vector<std::string> int_names; - Vector<int> int_flags; for (unsigned i = 0, n = species_names.size(); i < n; ++i) { auto& pc = allcontainers[i]; if (pc->plot_species) { Vector<std::string> real_names; + Vector<std::string> int_names; + Vector<int> int_flags; + real_names.push_back("weight"); real_names.push_back("momentum_x"); @@ -102,6 +103,15 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("theta"); #endif + if(pc->do_field_ionization){ + int_names.push_back("ionization_level"); + // int_flags specifies, for each integer attribs, whether it is + // dumped to plotfiles. So far, ionization_level is the only + // integer attribs, and it is automatically dumped to plotfiles + // when ionization is on. + int_flags.resize(1, 1); + } + // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index ff55724f8..75643d748 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -276,6 +276,9 @@ WarpX::EvolveEM (int numsteps) void WarpX::OneStep_nosub (Real cur_time) { + // Loop over species. For each ionizable species, create particles in + // product species. + mypc->doFieldIonization(); // Push particle from x^{n} to x^{n+1} // from p^{n-1/2} to p^{n+1/2} // Deposit current j^{n+1/2} @@ -339,6 +342,10 @@ WarpX::OneStep_sub1 (Real curtime) { // TODO: we could save some charge depositions + // Loop over species. For each ionizable species, create particles in + // product species. + mypc->doFieldIonization(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels"); const int fine_lev = 1; const int coarse_lev = 0; diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 1df05bc0f..e390dd61a 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -159,33 +159,29 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) auto const& Eyfab = Ey->array(mfi); auto const& Ezfab = Ez->array(mfi); if (do_nodal) { - amrex::ParallelFor(tbx, + amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bx_nodal(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); - }); - amrex::ParallelFor(tby, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_by_nodal(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); }); } else if (WarpX::maxwell_fdtd_solver_id == 0) { - amrex::ParallelFor(tbx, + amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); - }); - amrex::ParallelFor(tby, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin); @@ -198,23 +194,21 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - amrex::ParallelFor(tbx, + amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bx_ckc(j,k,l,Bxfab,Eyfab,Ezfab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tby, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_by_ckc(j,k,l,Byfab,Exfab,Ezfab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bz_ckc(j,k,l,Bzfab,Exfab,Eyfab, @@ -356,33 +350,29 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) auto const& jzfab = jz->array(mfi); if (do_nodal) { - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_nodal(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_nodal(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2); }); } else { - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin); @@ -393,17 +383,15 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { auto const& Ffab = F->array(mfi); if (WarpX::maxwell_fdtd_solver_id == 0) { - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_f_yee(j,k,l,Exfab,Ffab,dtsdx_c2); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_f_yee(j,k,l,Eyfab,Ffab,dtsdy_c2); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_f_yee(j,k,l,Ezfab,Ffab,dtsdz_c2); @@ -417,23 +405,21 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_f_ckc(j,k,l,Exfab,Ffab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_f_ckc(j,k,l,Eyfab,Ffab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_f_ckc(j,k,l,Ezfab,Ffab, @@ -616,7 +602,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // Rescale current in r-z mode since the inverse volume factor was not // included in the current deposition. - amrex::ParallelFor(tbr, + amrex::ParallelFor(tbr, tbt, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Wrap the current density deposited in the guard cells around @@ -630,8 +616,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // between on axis and off-axis factors const amrex::Real r = std::abs(rmin + (i - irmin + 0.5)*dr); Jr_arr(i,j,0) /= (2.*MathConst::pi*r); - }); - amrex::ParallelFor(tbt, + }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Wrap the current density deposited in the guard cells around @@ -649,8 +634,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu } else { Jt_arr(i,j,0) /= (2.*MathConst::pi*r); } - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Wrap the current density deposited in the guard cells around diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 786ebc622..c141ae4b0 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -454,9 +454,12 @@ LaserParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev = nullptr; + DepositCharge(pti, wp, ion_lev, rho, 0, 0, + np_current, thread_num, lev, lev); if (crho) { - DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 0, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -509,13 +512,14 @@ LaserParticleContainer::Evolve (int lev, // Current Deposition // // Deposit inside domains - DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + int* ion_lev = nullptr; + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, lev, lev, dt); bool has_buffer = cjx; if (has_buffer){ // Deposit in buffers - DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, np_current, np-np_current, thread_num, lev, lev-1, dt); } @@ -528,9 +532,12 @@ LaserParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev = nullptr; + DepositCharge(pti, wp, ion_lev, rho, 1, 0, + np_current, thread_num, lev, lev); if (crho) { - DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 1, np_current, + np-np_current, thread_num, lev, lev-1); } } diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 26d42370e..4253f5e76 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -6,6 +6,10 @@ /* \brief Charge Deposition for thread thread_num * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param rho_arr : Array4 of charge density, either full array or tile. * \param np_to_depose : Number of particles for which current is deposited. * \param dx : 3D cell size @@ -18,6 +22,7 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const yp, const amrex::Real * const zp, const amrex::Real * const wp, + const int * const ion_lev, const amrex::Array4<amrex::Real>& rho_arr, const long np_to_depose, const std::array<amrex::Real,3>& dx, @@ -25,6 +30,9 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Dim3 lo, const amrex::Real q) { + // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer + // (do_ionization=1) + const bool do_ionization = ion_lev; const amrex::Real dxi = 1.0/dx[0]; const amrex::Real dzi = 1.0/dx[2]; #if (AMREX_SPACEDIM == 2) @@ -43,7 +51,10 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { // --- Get particle quantities - const amrex::Real wq = q*wp[ip]*invvol; + amrex::Real wq = q*wp[ip]*invvol; + if (do_ionization){ + wq *= ion_lev[ip]; + } // --- Compute shape factors // x direction diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index f1e30e1be..e8e7fab0a 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -7,6 +7,10 @@ * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx_arr : Array4 of current density, either full array or tile. * \param jy_arr : Array4 of current density, either full array or tile. * \param jz_arr : Array4 of current density, either full array or tile. @@ -26,6 +30,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, + const int * const ion_lev, const amrex::Array4<amrex::Real>& jx_arr, const amrex::Array4<amrex::Real>& jy_arr, const amrex::Array4<amrex::Real>& jz_arr, @@ -36,6 +41,9 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real stagger_shift, const amrex::Real q) { + // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer + // (do_ionization=1) + const bool do_ionization = ion_lev; const amrex::Real dxi = 1.0/dx[0]; const amrex::Real dzi = 1.0/dx[2]; const amrex::Real dts2dx = 0.5*dt*dxi; @@ -61,7 +69,10 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + uyp[ip]*uyp[ip]*clightsq + uzp[ip]*uzp[ip]*clightsq); - const amrex::Real wq = q*wp[ip]; + amrex::Real wq = q*wp[ip]; + if (do_ionization){ + wq *= ion_lev[ip]; + } const amrex::Real vx = uxp[ip]*gaminv; const amrex::Real vy = uyp[ip]*gaminv; const amrex::Real vz = uzp[ip]*gaminv; @@ -161,6 +172,10 @@ void doDepositionShapeN(const amrex::Real * const xp, * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param Jx_arr : Array4 of current density, either full array or tile. * \param Jy_arr : Array4 of current density, either full array or tile. * \param Jz_arr : Array4 of current density, either full array or tile. @@ -179,6 +194,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, + const int * ion_lev, const amrex::Array4<amrex::Real>& Jx_arr, const amrex::Array4<amrex::Real>& Jy_arr, const amrex::Array4<amrex::Real>& Jz_arr, @@ -189,6 +205,9 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Dim3 lo, const amrex::Real q) { + // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer + // (do_ionization=1) + const bool do_ionization = ion_lev; const amrex::Real dxi = 1.0/dx[0]; const amrex::Real dtsdx0 = dt*dxi; const amrex::Real xmin = xyzmin[0]; @@ -224,7 +243,10 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, + uzp[ip]*uzp[ip]*clightsq); // wqx, wqy wqz are particle current in each direction - const amrex::Real wq = q*wp[ip]; + amrex::Real wq = q*wp[ip]; + if (do_ionization){ + wq *= ion_lev[ip]; + } const amrex::Real wqx = wq*invdtdx; #if (defined WARPX_DIM_3D) const amrex::Real wqy = wq*invdtdy; diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 7c9ede411..aaed96423 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -7,7 +7,6 @@ #include <string> #include <AMReX_Particles.H> - #include <WarpXParticleContainer.H> #include <PhysicalParticleContainer.H> #include <RigidInjectedParticleContainer.H> @@ -128,6 +127,8 @@ public: /// std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false); + void doFieldIonization (); + void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; @@ -207,6 +208,9 @@ private: void ReadParameters (); + void mapSpeciesProduct (); + int getSpeciesID (std::string product_str); + // Number of species dumped in BoostedFrameDiagnostics int nspecies_boosted_frame_diags = 0; // map_species_boosted_frame_diags[i] is the species ID in diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 4980bc62a..143f4b7f9 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -115,6 +115,9 @@ MultiParticleContainer::InitData () pc->InitData(); } pc_tmp->InitData(); + // For each species, get the ID of its product species. + // This is used for ionization and pair creation processes. + mapSpeciesProduct(); } @@ -430,7 +433,7 @@ MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const } int -MultiParticleContainer::doContinuousInjection() const +MultiParticleContainer::doContinuousInjection () const { int warpx_do_continuous_injection = 0; for (int i=0; i<nspecies+nlasers; i++){ @@ -441,3 +444,254 @@ MultiParticleContainer::doContinuousInjection() const } return warpx_do_continuous_injection; } + +/* \brief Get ID of product species of each species. + * The users specifies the name of the product species, + * this routine get its ID. + */ +void +MultiParticleContainer::mapSpeciesProduct () +{ + for (int i=0; i<nspecies; i++){ + auto& pc = allcontainers[i]; + // If species pc has ionization on, find species with name + // pc->ionization_product_name and store its ID into + // pc->ionization_product. + if (pc->do_field_ionization){ + int i_product = getSpeciesID(pc->ionization_product_name); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + i != i_product, + "ERROR: ionization product cannot be the same species"); + pc->ionization_product = i_product; + } + } +} + +/* \brief Given a species name, return its ID. + */ +int +MultiParticleContainer::getSpeciesID (std::string product_str) +{ + int i_product; + bool found = 0; + // Loop over species + for (int i=0; i<nspecies; i++){ + // If species name matches, store its ID + // into i_product + if (species_names[i] == product_str){ + found = 1; + i_product = i; + } + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + found != 0, + "ERROR: could not find product species ID for ionization. Wrong name?"); + return i_product; +} + +namespace +{ + // For particle i in mfi, if is_ionized[i]=1, copy particle + // particle i from container pc_source into pc_product + void createIonizedParticles ( + int lev, const MFIter& mfi, + std::unique_ptr< WarpXParticleContainer>& pc_source, + std::unique_ptr< WarpXParticleContainer>& pc_product, + amrex::Gpu::ManagedDeviceVector<int>& is_ionized) + { + BL_PROFILE("createIonizedParticles"); + + const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + // Get source particle data + auto& ptile_source = pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + const int np_source = ptile_source.GetArrayOfStructs().size(); + if (np_source == 0) return; + // --- source AoS particle data + WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); + // --- source SoA particle data + auto& soa_source = ptile_source.GetStructOfArrays(); + GpuArray<Real*,PIdx::nattribs> attribs_source; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + attribs_source[ia] = soa_source.GetRealData(ia).data(); + } + // --- source runtime attribs + GpuArray<Real*,3> runtime_uold_source; + // Prepare arrays for boosted frame diagnostics. + runtime_uold_source[0] = soa_source.GetRealData(PIdx::ux).data(); + runtime_uold_source[1] = soa_source.GetRealData(PIdx::uy).data(); + runtime_uold_source[2] = soa_source.GetRealData(PIdx::uz).data(); + + // Indices of product particle for each ionized source particle. + // i_product[i]-1 is the location in product tile of product particle + // from source particle i. + amrex::Gpu::ManagedDeviceVector<int> i_product; + i_product.resize(np_source); + // 0<i<np_source + // 0<i_product<np_ionized + // Strictly speaking, i_product should be an exclusive_scan of + // is_ionized. However, for indices where is_ionized is 1, the + // inclusive scan gives the same result with an offset of 1. + // The advantage of inclusive_scan is that the sum of is_ionized + // is in the last element, so no other reduction is required to get + // number of particles. + // Gpu::inclusive_scan runs on the current GPU stream, and synchronizes + // with the CPU, so that the next line (executed by the CPU) has the + // updated values of i_product + amrex::Gpu::inclusive_scan(is_ionized.begin(), is_ionized.end(), i_product.begin()); + int np_ionized = i_product[np_source-1]; + if (np_ionized == 0) return; + int* AMREX_RESTRICT p_i_product = i_product.dataPtr(); + + // Get product particle data + auto& ptile_product = pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + // old and new (i.e., including ionized particles) number of particles + // for product species + const int np_product_old = ptile_product.GetArrayOfStructs().size(); + const int np_product_new = np_product_old + np_ionized; + // Allocate extra space in product species for ionized particles. + ptile_product.resize(np_product_new); + // --- product AoS particle data + // First element is the first newly-created product particle + WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; + // --- product SoA particle data + auto& soa_product = ptile_product.GetStructOfArrays(); + GpuArray<Real*,PIdx::nattribs> attribs_product; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + // First element is the first newly-created product particle + attribs_product[ia] = soa_product.GetRealData(ia).data() + np_product_old; + } + // --- product runtime attribs + GpuArray<Real*,6> runtime_attribs_product; + bool do_boosted_product = WarpX::do_boosted_frame_diagnostic + && pc_product->DoBoostedFrameDiags(); + if (do_boosted_product) { + std::map<std::string, int> comps_product = pc_product->getParticleComps(); + runtime_attribs_product[0] = soa_product.GetRealData(comps_product[ "xold"]).data() + np_product_old; + runtime_attribs_product[1] = soa_product.GetRealData(comps_product[ "yold"]).data() + np_product_old; + runtime_attribs_product[2] = soa_product.GetRealData(comps_product[ "zold"]).data() + np_product_old; + runtime_attribs_product[3] = soa_product.GetRealData(comps_product["uxold"]).data() + np_product_old; + runtime_attribs_product[4] = soa_product.GetRealData(comps_product["uyold"]).data() + np_product_old; + runtime_attribs_product[5] = soa_product.GetRealData(comps_product["uzold"]).data() + np_product_old; + } + + int pid_product; +#pragma omp critical (doFieldIonization_nextid) + { + // ID of first newly-created product particle + pid_product = pc_product->NextID(); + // Update NextID to include particles created in this function + pc_product->setNextID(pid_product+np_ionized); + } + const int cpuid = ParallelDescriptor::MyProc(); + + // Loop over all source particles. If is_ionized, copy particle data + // to corresponding product particle. + amrex::For( + np_source, [=] AMREX_GPU_DEVICE (int is) noexcept + { + if(p_is_ionized[is]){ + // offset of 1 due to inclusive scan + int ip = p_i_product[is]-1; + // is: index of ionized particle in source species + // ip: index of corresponding new particle in product species + WarpXParticleContainer::ParticleType& p_product = particles_product[ip]; + WarpXParticleContainer::ParticleType& p_source = particles_source[is]; + // Copy particle from source to product: AoS + p_product.id() = pid_product + ip; + p_product.cpu() = cpuid; + p_product.pos(0) = p_source.pos(0); + p_product.pos(1) = p_source.pos(1); +#if (AMREX_SPACEDIM == 3) + p_product.pos(2) = p_source.pos(2); +#endif + // Copy particle from source to product: SoA + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + attribs_product[ia][ip] = attribs_source[ia][is]; + } + // Update xold etc. if boosted frame diagnostics required + // for product species. Fill runtime attribs with a copy of + // current properties (xold = x etc.). + if (do_boosted_product) { + runtime_attribs_product[0][ip] = p_source.pos(0); + runtime_attribs_product[1][ip] = p_source.pos(1); + runtime_attribs_product[2][ip] = p_source.pos(2); + runtime_attribs_product[3][ip] = runtime_uold_source[0][ip]; + runtime_attribs_product[4][ip] = runtime_uold_source[1][ip]; + runtime_attribs_product[5][ip] = runtime_uold_source[2][ip]; + } + } + } + ); + } +} + +void +MultiParticleContainer::doFieldIonization () +{ + BL_PROFILE("MPC::doFieldIonization"); + // Loop over all species. + // Ionized particles in pc_source create particles in pc_product + for (auto& pc_source : allcontainers){ + + // Skip if not ionizable + if (!pc_source->do_field_ionization){ continue; } + + // Get product species + auto& pc_product = allcontainers[pc_source->ionization_product]; + + for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ + + // When using runtime components, AMReX requires to touch all tiles + // in serial and create particles tiles with runtime components if + // they do not exist (or if they were defined by default, i.e., + // without runtime component). +#ifdef _OPENMP + // Touch all tiles of source species in serial if runtime attribs + for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + if ( (pc_source->NumRuntimeRealComps()>0) || (pc_source->NumRuntimeIntComps()>0) ) { + pc_source->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + } +#endif + // Touch all tiles of product species in serial + for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + pc_product->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + + // Enable tiling + MFItInfo info; + if (pc_source->do_tiling && Gpu::notInLaunchRegion()) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + pc_product->do_tiling, + "For ionization, either all or none of the " + "particle species must use tiling."); + info.EnableTiling(pc_source->tile_size); + } + +#ifdef _OPENMP + info.SetDynamic(true); +#pragma omp parallel +#endif + // Loop over all grids (if not tiling) or grids and tiles (if tiling) + for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + { + // Ionization mask: one element per source particles. + // 0 if not ionized, 1 if ionized. + amrex::Gpu::ManagedDeviceVector<int> is_ionized; + pc_source->buildIonizationMask(mfi, lev, is_ionized); + // Create particles in pc_product + createIonizedParticles(lev, mfi, pc_source, pc_product, is_ionized); + } + } // lev + } // pc_source +} diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b80619733..c7494fbdf 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -22,6 +22,8 @@ public: virtual void InitData () override; + void InitIonizationModule (); + #ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks) override; @@ -104,6 +106,9 @@ public: virtual void PostRestart () final {} void SplitParticles(int lev); + + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) override; // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 95909feba..79a6f88f1 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -6,6 +6,7 @@ #include <WarpX.H> #include <WarpXConst.H> #include <WarpXWrappers.h> +#include <IonizationEnergiesTable.H> #include <FieldGather.H> #include <WarpXAlgorithmSelection.H> @@ -37,6 +38,14 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); + pp.query("do_field_ionization", do_field_ionization); +#ifdef AMREX_USE_GPU + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + do_field_ionization == 0, + "Field ionization does not work on GPU so far, because the current " + "version of Redistribute in AMReX does not work with runtime parameters"); +#endif + pp.query("plot_species", plot_species); int do_user_plot_vars; do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); @@ -77,6 +86,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) void PhysicalParticleContainer::InitData() { + // Init ionization module here instead of in the PhysicalParticleContainer + // constructor because dt is required + if (do_field_ionization) {InitIonizationModule();} AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute } @@ -196,6 +208,12 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) + { + // need to create old values + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + } + // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); } @@ -278,6 +296,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) GetParticles(lev)[index]; tmp_particle_data.resize(finestLevel()+1); tmp_particle_data[lev][index]; + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { + DefineAndReturnParticleTile(lev, mfi.index(), mfi.LocalTileIndex()); + } } #endif @@ -401,6 +422,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int cpuid = ParallelDescriptor::MyProc(); auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { + DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); @@ -412,6 +438,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pa[ia] = soa.GetRealData(ia).data() + old_size; } + int* pi; + if (do_field_ionization) { + pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; + } + const GpuArray<Real,AMREX_SPACEDIM> overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -420,6 +451,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) std::size_t shared_mem_bytes = plasma_injector->sharedMemoryNeeded(); int lrrfac = rrfac; + bool loc_do_field_ionization = do_field_ionization; + int loc_ionization_initial_level = ionization_initial_level; + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). // The invalid ones are given negative ID and are deleted during the @@ -540,6 +574,10 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) u.z = gamma_boost * ( u.z -beta_boost*gamma_lab ); } + if (loc_do_field_ionization) { + pi[ip] = loc_ionization_initial_level; + } + u.x *= PhysConst::c; u.y *= PhysConst::c; u.z *= PhysConst::c; @@ -1120,9 +1158,18 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + // Deposit charge before particle push, in component 0 of MultiFab rho. + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 0, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 0, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -1240,13 +1287,20 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt); } } else { + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + // Deposit inside domains - DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, lev, lev, dt); if (has_buffer){ // Deposit in buffers - DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, np_current, np-np_current, thread_num, lev, lev-1, dt); } @@ -1261,9 +1315,18 @@ PhysicalParticleContainer::Evolve (int lev, } if (rho) { - DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + // Deposit charge after particle push, in component 1 of MultiFab rho. + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 1, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 1, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -1468,25 +1531,38 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, copy_attribs(pti, x, y, z); } + int* AMREX_RESTRICT ion_lev = nullptr; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } + // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( pti.numParticles(), + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } UpdateMomentumBoris( ux[i], uy[i], uz[i], gi[i], - Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); UpdatePosition( x[i], y[i], z[i], ux[i], uy[i], uz[i], dt ); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { - amrex::ParallelFor( pti.numParticles(), + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i], - Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ux[i], uy[i], uz[i], dt ); } ); } else { @@ -1895,3 +1971,136 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } } } + + +void PhysicalParticleContainer::InitIonizationModule () +{ + if (!do_field_ionization) return; + ParmParse pp(species_name); + if (charge != PhysConst::q_e){ + amrex::Warning( + "charge != q_e for ionizable species: overriding user value and setting charge = q_e."); + charge = PhysConst::q_e; + } + pp.query("ionization_initial_level", ionization_initial_level); + pp.get("ionization_product_species", ionization_product_name); + pp.get("physical_element", physical_element); + // Add runtime integer component for ionization level + AddIntComp("ionization_level"); + // Get atomic number and ionization energies from file + int ion_element_id = ion_map_ids[physical_element]; + ion_atomic_number = ion_atomic_numbers[ion_element_id]; + ionization_energies.resize(ion_atomic_number); + int offset = ion_energy_offsets[ion_element_id]; + for(int i=0; i<ion_atomic_number; i++){ + ionization_energies[i] = table_ionization_energies[i+offset]; + } + // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) + // For now, we assume l=0 and m=0. + // The approximate expressions are used, + // without Gamma function + Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + std::pow(PhysConst::alpha,4)/PhysConst::r_e; + Real UH = table_ionization_energies[0]; + Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; + + const Real dt = WarpX::GetInstance().getdt(0); + + adk_power.resize(ion_atomic_number); + adk_prefactor.resize(ion_atomic_number); + adk_exp_prefactor.resize(ion_atomic_number); + for (int i=0; i<ion_atomic_number; ++i){ + Real n_eff = (i+1) * std::sqrt(UH/ionization_energies[i]); + Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); + adk_power[i] = -(2*n_eff - 1); + Real Uion = ionization_energies[i]; + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); + adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; + } +} + +/* \brief create mask of ionized particles (1 if ionized, 0 otherwise) + * + * \param mfi: tile or grid + * \param lev: MR level + * \param ionization_mask: Array with as many elements as particles in mfi. + * This function initialized the array, and set each element to 1 or 0 + * depending on whether the particle is ionized or not. + */ +void +PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) +{ + BL_PROFILE("PPC::buildIonizationMask"); + // Get pointers to ionization data from pc_source + const Real * const AMREX_RESTRICT p_ionization_energies = ionization_energies.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_prefactor = adk_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_exp_prefactor = adk_exp_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_power = adk_power.dataPtr(); + + // Current tile info + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + // Get GPU-friendly arrays of particle data + auto& ptile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + // Only need attribs (i.e., SoA data) + auto& soa = ptile.GetStructOfArrays(); + const int np = ptile.GetArrayOfStructs().size(); + + // If no particle, nothing to do. + if (np == 0) return; + // Otherwise, resize ionization_mask, and get poiters to attribs arrays. + ionization_mask.resize(np); + int * const AMREX_RESTRICT p_ionization_mask = ionization_mask.data(); + const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); + const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); + const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); + const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); + const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); + const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); + const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); + const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); + int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); + + Real c = PhysConst::c; + Real c2_inv = 1./c/c; + int atomic_number = ion_atomic_number; + + // Loop over all particles in grid/tile. If ionized, set mask to 1 + // and increment ionization level. + ParallelFor( + np, + [=] AMREX_GPU_DEVICE (long i) { + // Get index of ionization_level + p_ionization_mask[i] = 0; + if ( ion_lev[i]<atomic_number ){ + Real random_draw = amrex::Random(); + // Compute electric field amplitude in the particle's frame of + // reference (particularly important when in boosted frame). + Real ga = std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i]) * c2_inv); + Real E = std::sqrt( + - ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * c2_inv + + ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) * ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) + + ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) * ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) + + ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) * ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) + ); + // Compute probability of ionization p + Real w_dtau = 1./ ga * p_adk_prefactor[ion_lev[i]] * + std::pow(E,p_adk_power[ion_lev[i]]) * + std::exp( p_adk_exp_prefactor[ion_lev[i]]/E ); + Real p = 1. - std::exp( - w_dtau ); + + if (random_draw < p){ + // increment particle's ionization level + ion_lev[i] += 1; + // update mask + p_ionization_mask[i] = 1; + } + } + } + ); +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index a21506ec3..c4bb5e410 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -85,6 +85,10 @@ public: RealVector& GetAttribs (int comp) { return GetStructOfArrays().GetRealData(comp); } + + IntVector& GetiAttribs (int comp) { + return GetStructOfArrays().GetIntData(comp); + } }; class MultiParticleContainer; @@ -176,6 +180,7 @@ public: virtual void DepositCharge(WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, amrex::MultiFab* rho, int icomp, const long offset, @@ -189,6 +194,7 @@ public: RealVector& uxp, RealVector& uyp, RealVector& uzp, + const int * const ion_lev, amrex::MultiFab* jx, amrex::MultiFab* jy, amrex::MultiFab* jz, @@ -258,6 +264,8 @@ public: static int NextID () { return ParticleType::NextID(); } + void setNextID(int next_id) { ParticleType::NextID(next_id); } + bool do_splitting = false; // split along axes (0) or diagonals (1) @@ -274,15 +282,22 @@ public: void AddIntComp (const std::string& name, bool comm=true) { - particle_comps[name] = NumIntComps(); + particle_icomps[name] = NumIntComps(); AddIntComp(comm); } int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) + {}; + + std::map<std::string, int> getParticleComps () { return particle_comps;} + protected: std::map<std::string, int> particle_comps; + std::map<std::string, int> particle_icomps; int species_id; @@ -299,6 +314,17 @@ protected: // support all features allowed by direct injection. int do_continuous_injection = 0; + int do_field_ionization = 0; + int ionization_product; + std::string ionization_product_name; + int ion_atomic_number; + int ionization_initial_level = 0; + amrex::Gpu::ManagedVector<amrex::Real> ionization_energies; + amrex::Gpu::ManagedVector<amrex::Real> adk_power; + amrex::Gpu::ManagedVector<amrex::Real> adk_prefactor; + amrex::Gpu::ManagedVector<amrex::Real> adk_exp_prefactor; + std::string physical_element; + int do_boosted_frame_diags = 1; amrex::Vector<amrex::FArrayBox> local_rho; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index db0f683c7..b34e71178 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -227,6 +227,10 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + } + particle_tile.push_back(p); } @@ -237,6 +241,10 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + } + for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_DIM_RZ @@ -385,6 +393,10 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, * \param pti : Particle iterator * \param wp : Array of particle weights * \param uxp uyp uzp : Array of particle + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx jy jz : Full array of current density * \param offset : Index of first particle for which current is deposited * \param np_to_depose: Number of particles for which current is deposited. @@ -398,6 +410,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, + const int * const ion_lev, MultiFab* jx, MultiFab* jy, MultiFab* jz, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, @@ -480,37 +493,40 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ - doEsirkepovDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<1>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doEsirkepovDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<2>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doEsirkepovDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<3>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } } else { if (WarpX::nox == 1){ - doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<1>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } else if (WarpX::nox == 2){ - doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<2>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } else if (WarpX::nox == 3){ - doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<3>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } } @@ -525,8 +541,27 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #endif } +/* \brief Charge Deposition for thread thread_num + * \param pti : Particle iterator + * \param wp : Array of particle weights + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. + * \param rho : Full array of charge density + * \param icomp : Component of rho into which charge is deposited. + 0: old value (before particle push). + 1: new value (after particle push). + * \param offset : Index of first particle for which charge is deposited + * \param np_to_depose: Number of particles for which charge is deposited. + Particles [offset,offset+np_tp_depose] deposit charge + * \param thread_num : Thread number (if tiling) + * \param lev : Level of box that contains particles + * \param depos_lev : Level on which particles deposit (if buffers are used) + */ void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) @@ -588,14 +623,14 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, BL_PROFILE_VAR_START(blp_ppc_chd); if (WarpX::nox == 1){ - doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); @@ -705,7 +740,15 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + + DepositCharge(pti, wp, ion_lev, rho.get(), 0, 0, np, + thread_num, lev, lev); } #ifdef _OPENMP } diff --git a/Source/Utils/IonizationEnergiesTable.H b/Source/Utils/IonizationEnergiesTable.H new file mode 100644 index 000000000..475236abd --- /dev/null +++ b/Source/Utils/IonizationEnergiesTable.H @@ -0,0 +1,136 @@ +#include <map> +#include <AMReX_AmrCore.H> + +#ifndef WARPX_IONIZATION_TABLE_H_ +#define WARPX_IONIZATION_TABLE_H_ + +std::map<std::string, int> ion_map_ids = { + {"H", 0}, + {"He", 1}, + {"Li", 2}, + {"Be", 3}, + {"B", 4}, + {"C", 5}, + {"N", 6}, + {"O", 7}, + {"F", 8}, + {"Ne", 9}, + {"Na", 10}, + {"Mg", 11}, + {"Al", 12}, + {"Si", 13}, + {"P", 14}, + {"S", 15}, + {"Cl", 16}, + {"Ar", 17}, + {"Kr", 18}, + {"Rb", 19}, + {"Xe", 20}, + {"Rn", 21} }; + +const int nelements = 22; + +const int ion_atomic_numbers[nelements] = { + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, + 11, 12, 13, 14, 15, 16, 17, 18, 36, 37, + 54, 86}; + +const int ion_energy_offsets[nelements] = { + 0, 1, 3, 6, 10, 15, 21, 28, 36, 45, + 55, 66, 78, 91, 105, 120, 136, 153, 171, 207, + 244, 298}; + +const int energies_tab_length = 384; + +const amrex::Real table_ionization_energies[energies_tab_length]{ + // H + 13.59843449, + // He + 24.58738880, 54.4177650, + // Li + 5.39171495, 75.6400964, 122.4543581, + // Be + 9.322699, 18.21115, 153.896203, 217.7185843, + // B + 8.298019, 25.15483, 37.93058, 259.3715, 340.226020, + // C + 11.2602880, 24.383154, 47.88778, 64.49352, 392.090515, 489.993194, + // N + 14.53413, 29.60125, 47.4453, 77.4735, 97.8901, 552.06732, 667.046116, + // O + 13.618055, 35.12112, 54.93554, 77.41350, 113.8990, 138.1189, 739.32682, + 871.40988, + // F + 17.42282, 34.97081, 62.70798, 87.175, 114.249, 157.16311, 185.1868, + 953.89804, 1103.11747, + // Ne + 21.564540, 40.96297, 63.4233, 97.1900, 126.247, 157.934, 207.271, + 239.0970, 1195.80783, 1362.19915, + // Na + 5.1390769, 47.28636, 71.6200, 98.936, 138.404, 172.23, 208.504, + 264.192, 299.856, 1465.13449, 1648.70218, + // Mg + 7.646236, 15.035271, 80.1436, 109.2654, 141.33, 186.76, 225.02, + 265.924, 327.99, 367.489, 1761.80487, 1962.66365, + // Al + 5.985769, 18.82855, 28.447642, 119.9924, 153.8252, 190.49, 241.76, + 284.64, 330.21, 398.65, 442.005, 2085.97700, 2304.14005, + // Si + 8.15168, 16.34585, 33.49300, 45.14179, 166.767, 205.279, 246.57, + 303.59, 351.28, 401.38, 476.273, 523.415, 2437.65813, 2673.17753, + // P + 10.486686, 19.76949, 30.20264, 51.44387, 65.02511, 220.430, 263.57, + 309.60, 372.31, 424.40, 479.44, 560.62, 611.741, 2816.90876, + 3069.8415, + // S + 10.36001, 23.33788, 34.86, 47.222, 72.5945, 88.0529, 280.954, + 328.794, 379.84, 447.7, 504.55, 564.41, 651.96, 706.994, + 3223.7807, 3494.1879, + // Cl + 12.967632, 23.81364, 39.80, 53.24, 67.68, 96.94, 114.2013, + 348.306, 400.851, 456.7, 530.0, 591.58, 656.30, 750.23, + 809.198, 3658.3437, 3946.2909, + // Ar + 15.7596117, 27.62967, 40.735, 59.58, 74.84, 91.290, 124.41, + 143.4567, 422.60, 479.76, 540.4, 619.0, 685.5, 755.13, + 855.5, 918.375, 4120.6656, 4426.2228, + // Kr + 13.9996053, 24.35984, 35.838, 50.85, 64.69, 78.49, 109.13, + 125.802, 233.0, 268, 308, 350, 391, 446, + 492, 540, 591, 640, 785, 831.6, 882.8, + 945, 999.0, 1042, 1155.0, 1205.23, 2928.9, 3072, + 3228, 3380, 3584, 3752.0, 3971, 4109.083, 17296.420, + 17936.209, + // Rb + 4.1771280, 27.28954, 39.247, 52.20, 68.44, 82.9, 98.67, + 132.79, 150.628, 277.12, 313.1, 356.0, 400, 443, + 502, 550, 601, 654, 706.0, 857, 905.3, + 958.9, 1024, 1080, 1125, 1242.5, 1294.57, 3133.3, + 3281, 3443, 3600, 3815, 3988, 4214, 4356.865, + 18305.884, 18965.516, + // Xe + 12.1298436, 20.975, 31.05, 42.20, 54.1, 66.703, 91.6, + 105.9778, 179.84, 202.0, 229.02, 255.0, 281, 314, + 343, 374, 404, 434, 549, 582, 616, + 650, 700, 736, 818, 857.0, 1493, 1571, + 1653, 1742, 1826, 1919, 2023, 2113, 2209, + 2300, 2556, 2637, 2726, 2811, 2975, 3068, + 3243, 3333.8, 7660, 7889, 8144, 8382, 8971, + 9243, 9581, 9810.37, 40271.724, 41299.71, + // Rn + 10.74850, 21.4, 29.4, 36.9, 52.9, 64.0, 88.0, + 102.0, 154.0, 173.9, 195.0, 218.0, 240, 264, + 293, 317, 342, 367, 488, 520, 550, + 580, 640, 680, 760, 800, 850, 920, + 980, 1050, 1110, 1180, 1250, 1310, 1390, + 1460, 1520, 1590, 1660, 1720, 2033, 2094, + 2158, 2227, 2293, 2357, 2467, 2535, 2606, + 2674, 2944, 3010, 3082, 3149, 3433, 3510, + 3699, 3777, 6169, 6318, 6476, 6646, 6807, + 6964, 7283, 7450, 7630, 7800, 8260, 8410, + 8570, 8710, 9610, 9780, 10120, 10290, 21770, + 22160, 22600, 22990, 26310, 26830, 27490, 27903.1, + 110842.0, 112843.7 }; + +#endif // #ifndef WARPX_IONIZATION_TABLE_H_ + diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index cd335dbcf..389b3670a 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -7,6 +7,7 @@ CEXE_headers += WarpXUtil.H CEXE_headers += WarpXAlgorithmSelection.H CEXE_sources += WarpXAlgorithmSelection.cpp CEXE_headers += NCIGodfreyTables.H +CEXE_headers += IonizationEnergiesTable.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index 4e226e55f..6901d4629 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -3,20 +3,24 @@ #include <AMReX_REAL.H> -// Physical constants -namespace PhysConst +// Math constants +namespace MathConst { - static constexpr amrex::Real c = 299792458.; - static constexpr amrex::Real ep0 = 8.854187817e-12; - static constexpr amrex::Real mu0 = 1.2566370614359173e-06; - static constexpr amrex::Real q_e = 1.6021764620000001e-19; - static constexpr amrex::Real m_e = 9.10938291e-31; - static constexpr amrex::Real m_p = 1.6726231000000001e-27; + static constexpr amrex::Real pi = 3.14159265358979323846; } -namespace MathConst +// Physical constants +namespace PhysConst { - static constexpr amrex::Real pi = 3.14159265358979323846; + static constexpr amrex::Real c = 299792458.; + static constexpr amrex::Real ep0 = 8.854187817e-12; + static constexpr amrex::Real mu0 = 1.2566370614359173e-06; + static constexpr amrex::Real q_e = 1.6021764620000001e-19; + static constexpr amrex::Real m_e = 9.10938291e-31; + static constexpr amrex::Real m_p = 1.6726231000000001e-27; + static constexpr amrex::Real hbar = 1.05457180010e-34; + static constexpr amrex::Real alpha = mu0/(4*MathConst::pi)*q_e*q_e*c/hbar; + static constexpr amrex::Real r_e = 1./(4*MathConst::pi*ep0) * q_e*q_e/(m_e*c*c); } #endif diff --git a/Source/Utils/atomic_data.txt b/Source/Utils/atomic_data.txt new file mode 100644 index 000000000..140f5a26a --- /dev/null +++ b/Source/Utils/atomic_data.txt @@ -0,0 +1,427 @@ +# Reference: +# Kramida, A., Ralchenko, Yu., Reader, J., and NIST ASD Team (2014). +# NIST Atomic Spectra Database (ver. 5.2), [Online]. +# Available: http://physics.nist.gov/asd [2017, March 3]. + +# License (from https://www.nist.gov/director/licensing): +# This data was developed by employees of the National Institute of Standards +# and Technology (NIST), an agency of the Federal Government. Pursuant to title +# 17 United States Code Section 105, works of NIST employees are not subject to +# copyright protection in the United States and are considered to be in the +# public domain. +# The data is provided by NIST as a public service and is expressly provided +# "AS IS." NIST MAKES NO WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR STATUTORY, +# INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTY OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE, NON-INFRINGEMENT AND DATA ACCURACY. NIST +# does not warrant or make any representations regarding the use of the data or +# the results thereof, including but not limited to the correctness, accuracy, +# reliability or usefulness of the data. NIST SHALL NOT BE LIABLE AND YOU +# HEREBY RELEASE NIST FROM LIABILITY FOR ANY INDIRECT, CONSEQUENTIAL, SPECIAL, +# OR INCIDENTAL DAMAGES (INCLUDING DAMAGES FOR LOSS OF BUSINESS PROFITS, +# BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, AND THE LIKE), WHETHER +# ARISING IN TORT, CONTRACT, OR OTHERWISE, ARISING FROM OR RELATING TO THE DATA +# (OR THE USE OF OR INABILITY TO USE THIS DATA), EVEN IF NIST HAS BEEN ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGES. +# To the extent that NIST may hold copyright in countries other than the United +# States, you are hereby granted the non-exclusive irrevocable and unconditional +# right to print, publish, prepare derivative works and distribute the NIST +# data, in any medium, or authorize others to do so on your behalf, on a +# royalty-free basis throughout the world. +# You may improve, modify, and create derivative works of the data or any +# portion of the data, and you may copy and distribute such modifications or +# works. Modified works should carry a notice stating that you changed the data +# and should note the date and nature of any such change. Please explicitly +# acknowledge the National Institute of Standards and Technology as the source +# of the data: Data citation recommendations are provided below. +# Permission to use this data is contingent upon your acceptance of the terms +# of this agreement and upon your providing appropriate acknowledgments of +# NIST's creation of the data. + +----------------------------------------------------------------------------- +At. num | Sp. Name | Ion Charge | Ionization Energy (eV) | +--------|---------------|------------|--------------------------------------| + 1 | H I | 0 | (13.59843449) | + 2 | He I | 0 | 24.58738880 | + 2 | He II | +1 | (54.4177650) | + 3 | Li I | 0 | 5.39171495 | + 3 | Li II | +1 | [75.6400964] | + 3 | Li III | +2 | (122.4543581) | + 4 | Be I | 0 | 9.322699 | + 4 | Be II | +1 | 18.21115 | + 4 | Be III | +2 | [153.896203] | + 4 | Be IV | +3 | (217.7185843) | + 5 | B I | 0 | 8.298019 | + 5 | B II | +1 | 25.15483 | + 5 | B III | +2 | 37.93058 | + 5 | B IV | +3 | 259.3715 | + 5 | B V | +4 | (340.226020) | + 6 | C I | 0 | 11.2602880 | + 6 | C II | +1 | 24.383154 | + 6 | C III | +2 | 47.88778 | + 6 | C IV | +3 | 64.49352 | + 6 | C V | +4 | [392.090515] | + 6 | C VI | +5 | (489.993194) | + 7 | N I | 0 | 14.53413 | + 7 | N II | +1 | [29.60125] | + 7 | N III | +2 | [47.4453] | + 7 | N IV | +3 | 77.4735 | + 7 | N V | +4 | 97.8901 | + 7 | N VI | +5 | [552.06732] | + 7 | N VII | +6 | (667.046116) | + 8 | O I | 0 | 13.618055 | + 8 | O II | +1 | 35.12112 | + 8 | O III | +2 | [54.93554] | + 8 | O IV | +3 | 77.41350 | + 8 | O V | +4 | 113.8990 | + 8 | O VI | +5 | [138.1189] | + 8 | O VII | +6 | [739.32682] | + 8 | O VIII | +7 | (871.40988) | + 9 | F I | 0 | 17.42282 | + 9 | F II | +1 | 34.97081 | + 9 | F III | +2 | [62.70798] | + 9 | F IV | +3 | [87.175] | + 9 | F V | +4 | [114.249] | + 9 | F VI | +5 | 157.16311 | + 9 | F VII | +6 | [185.1868] | + 9 | F VIII | +7 | [953.89804] | + 9 | F IX | +8 | (1103.11747) | + 10 | Ne I | 0 | 21.564540 | + 10 | Ne II | +1 | 40.96297 | + 10 | Ne III | +2 | 63.4233 | + 10 | Ne IV | +3 | [97.1900] | + 10 | Ne V | +4 | 126.247 | + 10 | Ne VI | +5 | 157.934 | + 10 | Ne VII | +6 | 207.271 | + 10 | Ne VIII | +7 | [239.0970] | + 10 | Ne IX | +8 | [1195.80783] | + 10 | Ne X | +9 | (1362.19915) | + 11 | Na I | 0 | 5.1390769 | + 11 | Na II | +1 | 47.28636 | + 11 | Na III | +2 | [71.6200] | + 11 | Na IV | +3 | [98.936] | + 11 | Na V | +4 | 138.404 | + 11 | Na VI | +5 | [172.23] | + 11 | Na VII | +6 | [208.504] | + 11 | Na VIII | +7 | [264.192] | + 11 | Na IX | +8 | (299.856) | + 11 | Na X | +9 | (1465.13449) | + 11 | Na XI | +10 | (1648.70218) | + 12 | Mg I | 0 | 7.646236 | + 12 | Mg II | +1 | 15.035271 | + 12 | Mg III | +2 | 80.1436 | + 12 | Mg IV | +3 | [109.2654] | + 12 | Mg V | +4 | [141.33] | + 12 | Mg VI | +5 | [186.76] | + 12 | Mg VII | +6 | [225.02] | + 12 | Mg VIII | +7 | [265.924] | + 12 | Mg IX | +8 | [327.99] | + 12 | Mg X | +9 | (367.489) | + 12 | Mg XI | +10 | (1761.80487) | + 12 | Mg XII | +11 | (1962.66365) | + 13 | Al I | 0 | 5.985769 | + 13 | Al II | +1 | 18.82855 | + 13 | Al III | +2 | 28.447642 | + 13 | Al IV | +3 | 119.9924 | + 13 | Al V | +4 | [153.8252] | + 13 | Al VI | +5 | [190.49] | + 13 | Al VII | +6 | [241.76] | + 13 | Al VIII | +7 | [284.64] | + 13 | Al IX | +8 | [330.21] | + 13 | Al X | +9 | [398.65] | + 13 | Al XI | +10 | (442.005) | + 13 | Al XII | +11 | (2085.97700) | + 13 | Al XIII | +12 | (2304.14005) | + 14 | Si I | 0 | 8.15168 | + 14 | Si II | +1 | 16.34585 | + 14 | Si III | +2 | 33.49300 | + 14 | Si IV | +3 | 45.14179 | + 14 | Si V | +4 | 166.767 | + 14 | Si VI | +5 | [205.279] | + 14 | Si VII | +6 | [246.57] | + 14 | Si VIII | +7 | [303.59] | + 14 | Si IX | +8 | [351.28] | + 14 | Si X | +9 | [401.38] | + 14 | Si XI | +10 | [476.273] | + 14 | Si XII | +11 | (523.415) | + 14 | Si XIII | +12 | (2437.65813) | + 14 | Si XIV | +13 | (2673.17753) | + 15 | P I | 0 | 10.486686 | + 15 | P II | +1 | [19.76949] | + 15 | P III | +2 | 30.20264 | + 15 | P IV | +3 | 51.44387 | + 15 | P V | +4 | 65.02511 | + 15 | P VI | +5 | [220.430] | + 15 | P VII | +6 | [263.57] | + 15 | P VIII | +7 | [309.60] | + 15 | P IX | +8 | [372.31] | + 15 | P X | +9 | [424.40] | + 15 | P XI | +10 | [479.44] | + 15 | P XII | +11 | [560.62] | + 15 | P XIII | +12 | (611.741) | + 15 | P XIV | +13 | (2816.90876) | + 15 | P XV | +14 | (3069.8415) | + 16 | S I | 0 | 10.36001 | + 16 | S II | +1 | 23.33788 | + 16 | S III | +2 | [34.86] | + 16 | S IV | +3 | [47.222] | + 16 | S V | +4 | 72.5945 | + 16 | S VI | +5 | 88.0529 | + 16 | S VII | +6 | 280.954 | + 16 | S VIII | +7 | 328.794 | + 16 | S IX | +8 | [379.84] | + 16 | S X | +9 | [447.7] | + 16 | S XI | +10 | [504.55] | + 16 | S XII | +11 | [564.41] | + 16 | S XIII | +12 | [651.96] | + 16 | S XIV | +13 | (706.994) | + 16 | S XV | +14 | (3223.7807) | + 16 | S XVI | +15 | (3494.1879) | + 17 | Cl I | 0 | 12.967632 | + 17 | Cl II | +1 | 23.81364 | + 17 | Cl III | +2 | [39.80] | + 17 | Cl IV | +3 | [53.24] | + 17 | Cl V | +4 | [67.68] | + 17 | Cl VI | +5 | [96.94] | + 17 | Cl VII | +6 | 114.2013 | + 17 | Cl VIII | +7 | 348.306 | + 17 | Cl IX | +8 | 400.851 | + 17 | Cl X | +9 | [456.7] | + 17 | Cl XI | +10 | [530.0] | + 17 | Cl XII | +11 | [591.58] | + 17 | Cl XIII | +12 | [656.30] | + 17 | Cl XIV | +13 | [750.23] | + 17 | Cl XV | +14 | (809.198) | + 17 | Cl XVI | +15 | (3658.3437) | + 17 | Cl XVII | +16 | (3946.2909) | + 18 | Ar I | 0 | 15.7596117 | + 18 | Ar II | +1 | 27.62967 | + 18 | Ar III | +2 | [40.735] | + 18 | Ar IV | +3 | [59.58] | + 18 | Ar V | +4 | [74.84] | + 18 | Ar VI | +5 | 91.290 | + 18 | Ar VII | +6 | [124.41] | + 18 | Ar VIII | +7 | [143.4567] | + 18 | Ar IX | +8 | [422.60] | + 18 | Ar X | +9 | [479.76] | + 18 | Ar XI | +10 | [540.4] | + 18 | Ar XII | +11 | [619.0] | + 18 | Ar XIII | +12 | [685.5] | + 18 | Ar XIV | +13 | [755.13] | + 18 | Ar XV | +14 | [855.5] | + 18 | Ar XVI | +15 | (918.375) | + 18 | Ar XVII | +16 | (4120.6656) | + 18 | Ar XVIII | +17 | (4426.2228) | + 36 | Kr I | 0 | 13.9996053 | + 36 | Kr II | +1 | 24.35984 | + 36 | Kr III | +2 | 35.838 | + 36 | Kr IV | +3 | 50.85 | + 36 | Kr V | +4 | [64.69] | + 36 | Kr VI | +5 | [78.49] | + 36 | Kr VII | +6 | (109.13) | + 36 | Kr VIII | +7 | 125.802 | + 36 | Kr IX | +8 | [233.0] | + 36 | Kr X | +9 | (268) | + 36 | Kr XI | +10 | (308) | + 36 | Kr XII | +11 | (350) | + 36 | Kr XIII | +12 | (391) | + 36 | Kr XIV | +13 | (446) | + 36 | Kr XV | +14 | (492) | + 36 | Kr XVI | +15 | (540) | + 36 | Kr XVII | +16 | (591) | + 36 | Kr XVIII | +17 | (640) | + 36 | Kr XIX | +18 | [785] | + 36 | Kr XX | +19 | [831.6] | + 36 | Kr XXI | +20 | [882.8] | + 36 | Kr XXII | +21 | [945] | + 36 | Kr XXIII | +22 | [999.0] | + 36 | Kr XXIV | +23 | [1042] | + 36 | Kr XXV | +24 | [1155.0] | + 36 | Kr XXVI | +25 | [1205.23] | + 36 | Kr XXVII | +26 | [2928.9] | + 36 | Kr XXVIII | +27 | [3072] | + 36 | Kr XXIX | +28 | [3228] | + 36 | Kr XXX | +29 | [3380] | + 36 | Kr XXXI | +30 | [3584] | + 36 | Kr XXXII | +31 | [3752.0] | + 36 | Kr XXXIII | +32 | [3971] | + 36 | Kr XXXIV | +33 | (4109.083) | + 36 | Kr XXXV | +34 | (17296.420) | + 36 | Kr XXXVI | +35 | (17936.209) | + 37 | Rb I | 0 | 4.1771280 | + 37 | Rb II | +1 | 27.28954 | + 37 | Rb III | +2 | [39.247] | + 37 | Rb IV | +3 | [52.20] | + 37 | Rb V | +4 | (68.44) | + 37 | Rb VI | +5 | (82.9) | + 37 | Rb VII | +6 | [98.67] | + 37 | Rb VIII | +7 | [132.79] | + 37 | Rb IX | +8 | 150.628 | + 37 | Rb X | +9 | 277.12 | + 37 | Rb XI | +10 | (313.1) | + 37 | Rb XII | +11 | (356.0) | + 37 | Rb XIII | +12 | (400) | + 37 | Rb XIV | +13 | (443) | + 37 | Rb XV | +14 | (502) | + 37 | Rb XVI | +15 | (550) | + 37 | Rb XVII | +16 | (601) | + 37 | Rb XVIII | +17 | (654) | + 37 | Rb XIX | +18 | (706.0) | + 37 | Rb XX | +19 | [857] | + 37 | Rb XXI | +20 | [905.3] | + 37 | Rb XXII | +21 | [958.9] | + 37 | Rb XXIII | +22 | [1024] | + 37 | Rb XXIV | +23 | [1080] | + 37 | Rb XXV | +24 | [1125] | + 37 | Rb XXVI | +25 | [1242.5] | + 37 | Rb XXVII | +26 | [1294.57] | + 37 | Rb XXVIII | +27 | [3133.3] | + 37 | Rb XXIX | +28 | [3281] | + 37 | Rb XXX | +29 | [3443] | + 37 | Rb XXXI | +30 | [3600] | + 37 | Rb XXXII | +31 | [3815] | + 37 | Rb XXXIII | +32 | [3988] | + 37 | Rb XXXIV | +33 | [4214] | + 37 | Rb XXXV | +34 | (4356.865) | + 37 | Rb XXXVI | +35 | (18305.884) | + 37 | Rb XXXVII | +36 | (18965.516) | + 54 | Xe I | 0 | 12.1298436 | + 54 | Xe II | +1 | [20.975] | + 54 | Xe III | +2 | [31.05] | + 54 | Xe IV | +3 | 42.20 | + 54 | Xe V | +4 | [54.1] | + 54 | Xe VI | +5 | 66.703 | + 54 | Xe VII | +6 | [91.6] | + 54 | Xe VIII | +7 | 105.9778 | + 54 | Xe IX | +8 | 179.84 | + 54 | Xe X | +9 | (202.0) | + 54 | Xe XI | +10 | [229.02] | + 54 | Xe XII | +11 | (255.0) | + 54 | Xe XIII | +12 | (281) | + 54 | Xe XIV | +13 | (314) | + 54 | Xe XV | +14 | (343) | + 54 | Xe XVI | +15 | (374) | + 54 | Xe XVII | +16 | (404) | + 54 | Xe XVIII | +17 | (434) | + 54 | Xe XIX | +18 | (549) | + 54 | Xe XX | +19 | (582) | + 54 | Xe XXI | +20 | (616) | + 54 | Xe XXII | +21 | (650) | + 54 | Xe XXIII | +22 | (700) | + 54 | Xe XXIV | +23 | (736) | + 54 | Xe XXV | +24 | (818) | + 54 | Xe XXVI | +25 | [857.0] | + 54 | Xe XXVII | +26 | (1493) | + 54 | Xe XXVIII | +27 | (1571) | + 54 | Xe XXIX | +28 | (1653) | + 54 | Xe XXX | +29 | (1742) | + 54 | Xe XXXI | +30 | (1826) | + 54 | Xe XXXII | +31 | (1919) | + 54 | Xe XXXIII | +32 | (2023) | + 54 | Xe XXXIV | +33 | (2113) | + 54 | Xe XXXV | +34 | (2209) | + 54 | Xe XXXVI | +35 | (2300) | + 54 | Xe XXXVII | +36 | (2556) | + 54 | Xe XXXVIII | +37 | (2637) | + 54 | Xe XXXIX | +38 | (2726) | + 54 | Xe XL | +39 | (2811) | + 54 | Xe XLI | +40 | (2975) | + 54 | Xe XLII | +41 | (3068) | + 54 | Xe XLIII | +42 | (3243) | + 54 | Xe XLIV | +43 | [3333.8] | + 54 | Xe XLV | +44 | (7660) | + 54 | Xe XLVI | +45 | (7889) | + 54 | Xe XLVII | +46 | (8144) | + 54 | Xe XLVIII | +47 | (8382) | + 54 | Xe XLIX | +48 | (8971) | + 54 | Xe L | +49 | (9243) | + 54 | Xe LI | +50 | (9581) | + 54 | Xe LII | +51 | (9810.37) | + 54 | Xe LIII | +52 | (40271.724) | + 54 | Xe LIV | +53 | (41299.71) | + 86 | Rn I | 0 | 10.74850 | + 86 | Rn II | +1 | [21.4] | + 86 | Rn III | +2 | [29.4] | + 86 | Rn IV | +3 | (36.9) | + 86 | Rn V | +4 | (52.9) | + 86 | Rn VI | +5 | (64.0) | + 86 | Rn VII | +6 | (88.0) | + 86 | Rn VIII | +7 | (102.0) | + 86 | Rn IX | +8 | (154.0) | + 86 | Rn X | +9 | (173.9) | + 86 | Rn XI | +10 | (195.0) | + 86 | Rn XII | +11 | (218.0) | + 86 | Rn XIII | +12 | (240) | + 86 | Rn XIV | +13 | (264) | + 86 | Rn XV | +14 | (293) | + 86 | Rn XVI | +15 | (317) | + 86 | Rn XVII | +16 | (342) | + 86 | Rn XVIII | +17 | (367) | + 86 | Rn XIX | +18 | (488) | + 86 | Rn XX | +19 | (520) | + 86 | Rn XXI | +20 | (550) | + 86 | Rn XXII | +21 | (580) | + 86 | Rn XXIII | +22 | (640) | + 86 | Rn XXIV | +23 | (680) | + 86 | Rn XXV | +24 | (760) | + 86 | Rn XXVI | +25 | (800) | + 86 | Rn XXVII | +26 | (850) | + 86 | Rn XXVIII | +27 | (920) | + 86 | Rn XXIX | +28 | (980) | + 86 | Rn XXX | +29 | (1050) | + 86 | Rn XXXI | +30 | (1110) | + 86 | Rn XXXII | +31 | (1180) | + 86 | Rn XXXIII | +32 | (1250) | + 86 | Rn XXXIV | +33 | (1310) | + 86 | Rn XXXV | +34 | (1390) | + 86 | Rn XXXVI | +35 | (1460) | + 86 | Rn XXXVII | +36 | (1520) | + 86 | Rn XXXVIII | +37 | (1590) | + 86 | Rn XXXIX | +38 | (1660) | + 86 | Rn XL | +39 | (1720) | + 86 | Rn XLI | +40 | (2033) | + 86 | Rn XLII | +41 | (2094) | + 86 | Rn XLIII | +42 | (2158) | + 86 | Rn XLIV | +43 | (2227) | + 86 | Rn XLV | +44 | (2293) | + 86 | Rn XLVI | +45 | (2357) | + 86 | Rn XLVII | +46 | (2467) | + 86 | Rn XLVIII | +47 | (2535) | + 86 | Rn XLIX | +48 | (2606) | + 86 | Rn L | +49 | (2674) | + 86 | Rn LI | +50 | (2944) | + 86 | Rn LII | +51 | (3010) | + 86 | Rn LIII | +52 | (3082) | + 86 | Rn LIV | +53 | (3149) | + 86 | Rn LV | +54 | (3433) | + 86 | Rn LVI | +55 | (3510) | + 86 | Rn LVII | +56 | (3699) | + 86 | Rn LVIII | +57 | [3777] | + 86 | Rn LIX | +58 | (6169) | + 86 | Rn LX | +59 | (6318) | + 86 | Rn LXI | +60 | (6476) | + 86 | Rn LXII | +61 | (6646) | + 86 | Rn LXIII | +62 | (6807) | + 86 | Rn LXIV | +63 | (6964) | + 86 | Rn LXV | +64 | (7283) | + 86 | Rn LXVI | +65 | (7450) | + 86 | Rn LXVII | +66 | (7630) | + 86 | Rn LXVIII | +67 | (7800) | + 86 | Rn LXIX | +68 | (8260) | + 86 | Rn LXX | +69 | (8410) | + 86 | Rn LXXI | +70 | (8570) | + 86 | Rn LXXII | +71 | (8710) | + 86 | Rn LXXIII | +72 | (9610) | + 86 | Rn LXXIV | +73 | (9780) | + 86 | Rn LXXV | +74 | (10120) | + 86 | Rn LXXVI | +75 | (10290) | + 86 | Rn LXXVII | +76 | (21770) | + 86 | Rn LXXVIII | +77 | (22160) | + 86 | Rn LXXIX | +78 | (22600) | + 86 | Rn LXXX | +79 | (22990) | + 86 | Rn LXXXI | +80 | (26310) | + 86 | Rn LXXXII | +81 | (26830) | + 86 | Rn LXXXIII | +82 | (27490) | + 86 | Rn LXXXIV | +83 | (27903.1) | + 86 | Rn LXXXV | +84 | (110842.0) | + 86 | Rn LXXXVI | +85 | (112843.7) | +----------------------------------------------------------------------------- diff --git a/Source/Utils/write_atomic_data_cpp.py b/Source/Utils/write_atomic_data_cpp.py new file mode 100644 index 000000000..b085d50eb --- /dev/null +++ b/Source/Utils/write_atomic_data_cpp.py @@ -0,0 +1,73 @@ +''' +This python script reads ionization tables in atomic_data.txt (generated from +the NIST website) and extracts ionization levels into C++ file +IonizationEnergiesTable.H, which contains tables + metadata. +''' + +import re, os +import numpy as np +from scipy.constants import e + +filename = os.path.join( '.', 'atomic_data.txt' ) +with open(filename) as f: + text_data = f.read() + +# Read full table from file and get names, atomic numbers and offsets +# position in table of ionization energies for all species +regex_command = '\n\s+(\d+)\s+\|\s+([A-Z]+[a-z]*)\s+\w+\s+\|\s+\+*(\d+)\s+\|\s+\(*\[*(\d+\.*\d*)' +list_of_tuples = re.findall( regex_command, text_data ) +ion_atom_numbers = [int(i) for i in list(dict.fromkeys( [x[0] for x in list_of_tuples] ))] +ion_names = list(dict.fromkeys( [x[1] for x in list_of_tuples] )) +ion_offsets = np.concatenate(([0], np.cumsum(np.array(ion_atom_numbers)[:-1])), axis=0) + +# Head of CPP file +cpp_string = '' +cpp_string += '#include <map>\n' +cpp_string += '#include <AMReX_AmrCore.H>\n\n' +cpp_string += '#ifndef WARPX_IONIZATION_TABLE_H_\n' +cpp_string += '#define WARPX_IONIZATION_TABLE_H_\n\n' + +# Map each element to ID in table +cpp_string += 'std::map<std::string, int> ion_map_ids = {' +for count, name in enumerate(ion_names): + cpp_string += '\n {"' + name + '", ' + str(count) + '},' +cpp_string = cpp_string[:-1] +cpp_string += ' };\n\n' + +# Atomic number of each species +cpp_string += 'const int nelements = ' + str(len(ion_names)) + ';\n\n' +cpp_string += 'const int ion_atomic_numbers[nelements] = {' +for count, atom_num in enumerate(ion_atom_numbers): + if count%10==0: cpp_string += '\n ' + cpp_string += str(atom_num) + ', ' +cpp_string = cpp_string[:-2] +cpp_string += '};\n\n' + +# Offset of each element in table of ionization energies +cpp_string += 'const int ion_energy_offsets[nelements] = {' +for count, offset in enumerate(ion_offsets): + if count%10==0: cpp_string += '\n ' + cpp_string += str(offset) + ', ' +cpp_string = cpp_string[:-2] +cpp_string += '};\n\n' + +# Table of ionization energies +cpp_string += 'const int energies_tab_length = ' + str(len(list_of_tuples)) + ';\n\n' +cpp_string += 'const amrex::Real table_ionization_energies[energies_tab_length]{' +for element in ion_names: + cpp_string += '\n // ' + element + regex_command = \ + '\n\s+(\d+)\s+\|\s+%s\s+\w+\s+\|\s+\+*(\d+)\s+\|\s+\(*\[*(\d+\.*\d*)' \ + %element + list_of_tuples = re.findall( regex_command, text_data ) + for count, energy in enumerate([x[2] for x in list_of_tuples]): + if count%7==0: cpp_string += '\n ' + cpp_string += energy + ', ' +cpp_string = cpp_string[:-2] +cpp_string += ' };\n\n' + +# Write the string to file +cpp_string += '#endif // #ifndef WARPX_IONIZATION_TABLE_H_\n' +f= open("IonizationEnergiesTable.H","w") +f.write(cpp_string) +f.close() diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1b653fd7f..7b90225dc 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -38,7 +38,7 @@ Vector<int> WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::use_picsar_deposition = 1; +long WarpX::use_picsar_deposition = 0; long WarpX::current_deposition_algo; long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; diff --git a/Tools/read_raw_data.py b/Tools/read_raw_data.py index bca14b582..363c3a007 100644 --- a/Tools/read_raw_data.py +++ b/Tools/read_raw_data.py @@ -312,7 +312,7 @@ def _read_buffer(snapshot, header_fn, _component_names): for i in range(header.ncomp): comp_data = arr[i*size:(i+1)*size].reshape(shape, order='F') data = all_data[_component_names[i]] - data[[slice(l,h+1) for l, h in zip(lo, hi)]] = comp_data + data[tuple([slice(l,h+1) for l, h in zip(lo, hi)])] = comp_data all_data[_component_names[i]] = data return all_data |