Laser ionization

CMakeLists.txt

@@ -342,6 +342,12 @@ if(BUILD_TESTING)
                 WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
         )
 
+        add_test(NAME laser_ionization.1Rank
+                COMMAND bash ${HiPACE_SOURCE_DIR}/tests/laser_ionization.1Rank.sh
+                        $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR} ${HiPACE_COMPUTE}
+                WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+            )
+
         if (NOT HiPACE_COMPUTE STREQUAL CUDA)
 
             # These tests only run on CPU

docs/source/run/parameters.rst

@@ -446,6 +446,9 @@ When both are specified, the per-species value is used.
 * ``<plasma name>.can_ionize`` (`bool`) optional (default `0`)
     Whether this plasma can ionize. Can also be set to 1 by specifying ``<plasma name>.ionization_product``.
 
+* ``<plasma name>.can_laser_ionize`` (`bool`) optional (default `<plasma name>.can_ionize`)
+    Whether this plasma can be ionized by a laser.
+
 * ``<plasma name>.initial_ion_level`` (`int`) optional (default `-1`)
     The initial ionization state of the plasma. `0` for neutral gasses.
     If set, the plasma charge gets multiplied by this number. If the plasma species is not ionizable,

examples/laser_ionization/analysis_laser_ionization.py

@@ -0,0 +1,59 @@
+#! /usr/bin/env python3
+
+#
+# This file is part of HiPACE++.
+#
+# Authors: EyaDammak
+
+import argparse
+import numpy as np
+import math
+from openpmd_viewer import OpenPMDTimeSeries
+import statistics
+from scipy.constants import e as qe
+parser = argparse.ArgumentParser(
+    description='Script to analyze the equality of two simulations')
+parser.add_argument('--first',
+                    dest='first',
+                    required=True,
+                    help='Path to the directory containing output files')
+parser.add_argument('--second',
+                    dest='second',
+                    required=True,
+                    help='Path to the directory containing output files')
+args = parser.parse_args()
+
+ts_linear = OpenPMDTimeSeries(args.first)
+ts_circular = OpenPMDTimeSeries(args.second)
+
+a0_linear = 0.00885126
+a0_circular = 0.00787934
+
+nc = 1.75e27
+n0 = nc / 10000
+
+iteration = 0
+
+rho_elec_linear, _ = ts_linear.get_field(field='rho_elec', coord='z', iteration=iteration)
+rho_elec_mean_linear = np.mean(rho_elec_linear, axis=(1, 2))
+rho_average_linear = statistics.mean(rho_elec_mean_linear[0:10])
+fraction_linear = -rho_average_linear / qe / n0
+
+rho_elec_circular, _ = ts_circular.get_field(field='rho_elec', coord='z', iteration=iteration)
+rho_elec_mean_circular = np.mean(rho_elec_circular, axis=(1, 2))
+rho_average_circular = statistics.mean(rho_elec_mean_circular[0:10]) #average over a thickness in the ionized region
+fraction_circular = -rho_average_circular / qe / n0
+
+fraction_warpx_linear = 0.41014984 # result from WarpX simulation
+fraction_warpx_circular = 0.502250841 # result from WarpX simulation
+
+relative_diff_linear = np.abs( ( fraction_linear - fraction_warpx_linear ) / fraction_warpx_linear )
+relative_diff_circular = np.abs( ( fraction_circular - fraction_warpx_circular ) / fraction_warpx_circular )
+
+tolerance = 0.25
+print(f"fraction_warpx_linear = {fraction_warpx_linear}")
+print(f"fraction_hipace_linear = {fraction_linear}")
+print(f"fraction_warpx_circular = {fraction_warpx_circular}")
+print(f"fraction_hipace_circular = {fraction_circular}")
+
+assert ( (relative_diff_linear < tolerance) and (relative_diff_circular < tolerance) ), 'Test laser_ionization did not pass'

examples/laser_ionization/inputs_laser_ionization

@@ -0,0 +1,62 @@
+max_step = 0
+hipace.dt = 0
+hipace.verbose = 3
+
+amr.n_cell = 32 32 50
+
+
+my_constants.kp_inv = 10.e-6
+
+my_constants.nc = 1.75e27
+my_constants.n0 = nc/10000
+
+my_constants.a0 = 0.00885126
+
+my_constants.w0_um = 1.e8
+my_constants.tau_fs = 30/1.17741 
+my_constants.lambda0_nm = 800
+
+hipace.file_prefix = laser_ionization.1Rank
+hipace.do_tiling = 0
+hipace.deposit_rho_individual = 1
+
+geometry.prob_lo     = -10.*kp_inv   -10.*kp_inv   -15.*kp_inv #box shifted
+geometry.prob_hi     =  10.*kp_inv    10.*kp_inv    5.*kp_inv
+
+lasers.names = laser
+lasers.lambda0 = lambda0_nm * 1.e-9
+
+
+laser.a0 = a0
+laser.position_mean = 0. 0. 0
+laser.w0 = w0_um*1.e-6
+laser.tau = tau_fs*1.e-15
+laser.focal_distance = 50.e-6
+
+plasmas.names = elec ion
+
+elec.density(x,y,z) = n0
+elec.ppc = 0 0
+elec.u_mean = 0.0 0.0 0.0
+elec.element = electron
+elec.neutralize_background = false
+
+ion.density(x,y,z) = n0
+ion.ppc = 1 1
+ion.u_mean = 0.0 0.0 0.0
+ion.element = H
+ion.mass_Da = 1.008
+ion.initial_ion_level = 0
+ion.can_laser_ionize = 1
+ion.ionization_product = elec
+
+amr.max_level = 0
+
+diagnostic.output_period = 1
+
+hipace.depos_order_xy = 0
+
+boundary.field = Dirichlet
+boundary.particle = Periodic
+
+diagnostic.diag_type = xyz

src/Hipace.cpp

@@ -690,11 +690,14 @@ Hipace::SolveOneSlice (int islice, int step)
     // copy fields (and laser) to diagnostic array
     FillFieldDiagnostics(current_N_level, islice);
 
-    // plasma ionization
+    // plasma field ionization
     for (int lev=0; lev<current_N_level; ++lev) {
         m_multi_plasma.DoFieldIonization(lev, m_3D_geom[lev], m_fields);
     }
 
+    // plasma laser ionization
+    m_multi_plasma.DoLaserIonization(islice, m_multi_laser.GetLaserGeom(), m_multi_laser);
+
     // Push plasma particles
     for (int lev=0; lev<current_N_level; ++lev) {
         m_multi_plasma.AdvanceParticles(m_fields, m_3D_geom, false, lev);

src/laser/MultiLaser.H

@@ -171,6 +171,9 @@ public:
     /** Get the geometry of the Laser Box */
     const amrex::Geometry& GetLaserGeom () const { return m_laser_geom_3D; }
 
+    /** Whether the polarization is linear. If not, circular polarization is assumed */
+    bool LinearPolarization () const { return m_linear_polarization; }
+
     /** If the laser geometry includes this slice
      * \param[in] islice slice index
      */

src/particles/particles_utils/FieldGather.H

@@ -265,6 +265,80 @@ void doLaserGatherShapeN (const amrex::Real xp,
     }
 }
 
+/**
+ * \brief Laser field gather for a single particle
+ *
+ * \tparam depos_order_xy Order of the transverse shape factor for the field gather
+ * \param[in] xp Particle position x
+ * \param[in] yp Particle position y
+ * \param[in] Ap a field at particle position
+ * \param[in] ADxp d/dx a field at particle position
+ * \param[in] ADzetap d/dzeta a field at particle position
+ * \param[in] laser_arr slice array for WhichSlice::This
+ * \param[in] dx_inv inverse cell spacing in x direction
+ * \param[in] dy_inv inverse cell spacing in y direction
+ * \param[in] dzeta_inv inverse cell spacing in zeta direction
+ * \param[in] x_pos_offset offset for converting positions to indexes
+ * \param[in] y_pos_offset offset for converting positions to indexes
+ */
+template <int depos_order_xy>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void doLaserGatherShapeN (const amrex::Real xp,
+                          const amrex::Real yp,
+                          amrex::GpuComplex<amrex::Real>& Ap,
+                          amrex::GpuComplex<amrex::Real>& ADxp,
+                          amrex::GpuComplex<amrex::Real>& ADzetap,
+                          Array3<amrex::Real const> const& laser_arr,
+                          const amrex::Real dx_inv,
+                          const amrex::Real dy_inv,
+                          const amrex::Real dzeta_inv,
+                          const amrex::Real x_pos_offset,
+                          const amrex::Real y_pos_offset)
+{
+    using namespace amrex::literals;
+
+    // x,y direction
+    const amrex::Real x = (xp-x_pos_offset)*dx_inv;
+    const amrex::Real y = (yp-y_pos_offset)*dy_inv;
+
+    // --- Compute shape factors
+    // x direction
+    // j_cell leftmost cell in x that the particle touches. sx_cell shape factor along x
+    amrex::Real sx_cell[depos_order_xy + 1];
+    const int i_cell = compute_shape_factor<depos_order_xy>(sx_cell, x);
+
+    // y direction
+    amrex::Real sy_cell[depos_order_xy + 1];
+    const int j_cell = compute_shape_factor<depos_order_xy>(sy_cell, y);
+
+    // Gather field Aabssq, AabssqDxp, and AabssqDxp on particle from field on grid slice_arr
+    // the derivative is calculated on the fly
+    for (int iy=0; iy<=depos_order_xy; iy++){
+        for (int ix=0; ix<=depos_order_xy; ix++){
+            using namespace WhichLaserSlice;
+            const amrex::Real s_cell = sx_cell[ix]*sy_cell[iy];
+            const int i = i_cell+ix;
+            const int j = j_cell+iy;
+            Ap += amrex::GpuComplex<amrex::Real> {
+                s_cell * laser_arr(i, j, n00j00_r),
+                s_cell * laser_arr(i, j, n00j00_i)
+            };
+            ADxp += amrex::GpuComplex<amrex::Real> {
+                s_cell * 0.5_rt * dx_inv *
+                    (laser_arr(i + 1, j, n00j00_r) - laser_arr(i - 1, j, n00j00_r)),
+                s_cell * 0.5_rt * dx_inv *
+                    (laser_arr(i + 1, j, n00j00_i) - laser_arr(i - 1, j, n00j00_i))
+            };
+            ADzetap += amrex::GpuComplex<amrex::Real> {
+                s_cell * dzeta_inv *
+                    (laser_arr(i, j, n00jp1_r) - laser_arr(i, j, n00j00_r)),
+                s_cell * dzeta_inv *
+                    (laser_arr(i, j, n00jp1_i) - laser_arr(i, j, n00j00_i))
+            };
+        }
+    }
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void doGatherEz (const amrex::Real xp,
                      const amrex::Real yp,

src/particles/plasma/MultiPlasma.H

@@ -98,6 +98,15 @@ public:
      */
     void DoFieldIonization (const int lev, const amrex::Geometry& geom, const Fields& fields);
 
+    /** Calculates Ionization Probability from laser and makes new Plasma Particles
+     *
+     * \param[in] islice zeta slice index
+     * \param[in] laser_geom Geometry of the laser
+     * \param[in] laser the laser class
+     */
+    void DoLaserIonization (const int islice, const amrex::Geometry& laser_geom,
+                            const MultiLaser& laser);
+
     /** \brief whether any plasma species uses a neutralizing background, e.g. no ion motion */
     bool AnySpeciesNeutralizeBackground () const;
 

src/particles/plasma/MultiPlasma.cpp

@@ -46,16 +46,16 @@ MultiPlasma::InitData (amrex::Vector<amrex::BoxArray> slice_ba,
         plasma.InitData(gm[0]);
 
         if(plasma.m_can_ionize) {
-            PlasmaParticleContainer* plasma_product = nullptr;
             for (int i=0; i<m_names.size(); ++i) {
                 if(m_names[i] == plasma.m_product_name) {
-                    plasma_product = &m_all_plasmas[i];
+                    plasma.m_product_pc = &m_all_plasmas[i];
                 }
             }
-            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(plasma_product != nullptr,
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(plasma.m_product_pc != nullptr,
                 "Must specify a valid product plasma for ionization using ionization_product");
-            plasma.InitIonizationModule(gm[0], plasma_product,
-            Hipace::m_background_density_SI); // geometry only for dz
+
+            plasma.InitIonizationModule(gm[0],
+                Hipace::m_background_density_SI); // geometry only for dz
         }
     }
 }
@@ -126,6 +126,15 @@ MultiPlasma::DoFieldIonization (
     }
 }
 
+void
+MultiPlasma::DoLaserIonization (
+    const int islice, const amrex::Geometry& laser_geom, const MultiLaser& laser)
+{
+    for (auto& plasma : m_all_plasmas) {
+        plasma.LaserIonization(islice, laser_geom, laser, Hipace::m_background_density_SI);
+    }
+}
+
 bool
 MultiPlasma::AnySpeciesNeutralizeBackground () const
 {

src/particles/plasma/PlasmaParticleContainer.H

@@ -84,10 +84,9 @@ public:
     /** Initialize ADK prefactors of ionizable plasmas
      *
      * \param[in] geom Geometry of the simulation, to get the cell size
-     * \param[in] product_pc The electron plasma PC that this Ion Plasma ionizes to
      * \param[in] background_density_SI background plasma density (only needed for normalized units)
      */
-    void InitIonizationModule (const amrex::Geometry& geom, PlasmaParticleContainer* product_pc,
+    void InitIonizationModule (const amrex::Geometry& geom,
                                const amrex::Real background_density_SI);
 
     /** Calculates Ionization Probability and generates new Plasma Particles
@@ -102,6 +101,18 @@ public:
                            const Fields& fields,
                            const amrex::Real background_density_SI);
 
+    /** Calculates Ionization Probability from a Laser and generates new Plasma Particles
+     *
+     * \param[in] islice zeta slice index
+     * \param[in] laser_geom Geometry of the laser
+     * \param[in] laser the laser class
+     * \param[in] background_density_SI background plasma density (only needed for normalized units)
+     */
+    void LaserIonization (const int islice,
+                          const amrex::Geometry& laser_geom,
+                          const MultiLaser& laser,
+                          const amrex::Real background_density_SI);
+
     /** Reorder particles to speed-up current deposition
      * \param[in] islice zeta slice index
      */
@@ -149,6 +160,8 @@ public:
      */
     void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom);
 
+    // init:
+
     amrex::Parser m_parser; /**< owns data for m_density_func */
     amrex::ParserExecutor<3> m_density_func; /**< Density function for the plasma */
     amrex::Real m_min_density {0.}; /**< minimal density at which particles are injected */
@@ -176,11 +189,19 @@ public:
     amrex::Real m_temperature_in_ev {0.}; /**< Temperature of the plasma in eV */
     /** whether to add a neutralizing background of immobile particles of opposite charge */
     bool m_neutralize_background = true;
+
+    // general:
+
     amrex::Real m_mass = 0; /**< mass of each particle of this species */
     amrex::Real m_charge = 0; /**< charge of each particle of this species, per Ion level */
-    int m_init_ion_lev = -1; /**< initial Ion level of each particle */
     int m_n_subcycles = 1; /**< number of subcycles in the plasma particle push */
-    bool m_can_ionize = false; /**< whether this plasma can ionize */
+
+    // ionization:
+
+    bool m_can_ionize = false; /**< whether this plasma can ionize from anything */
+    bool m_can_field_ionize = false; /**< whether this plasma can ionize from the field */
+    bool m_can_laser_ionize = false; /**< whether this plasma can ionize from a laser */
+    int m_init_ion_lev = -1; /**< initial Ion level of each particle */
     std::string m_product_name = ""; /**< name of Ionization product plasma */
     PlasmaParticleContainer* m_product_pc = nullptr; /**< Ionization product plasma */
     /** to calculate Ionization probability with ADK formula */
@@ -189,10 +210,18 @@ public:
     amrex::Gpu::DeviceVector<amrex::Real> m_adk_exp_prefactor;
     /** to calculate Ionization probability with ADK formula */
     amrex::Gpu::DeviceVector<amrex::Real> m_adk_power;
+    /** to calculate laser Ionization probability with ADK formula */
+    amrex::Gpu::DeviceVector<amrex::Real> m_laser_adk_prefactor;
+
+    // plasma sorting/reordering:
+
     /** After how many slices the particles are reordered. 0: off */
     int m_reorder_period = 0;
     /** 2D reordering index type. 0: cell, 1: node, 2: both */
     amrex::IntVect m_reorder_idx_type = {0, 0, 0};
+
+    // in-situ diagnostic:
+
     /** How often the insitu plasma diagnostics should be computed and written
      * Default is 0, meaning no output */
     int m_insitu_period {0};