Laser ionization with electron momentum

src/particles/plasma/PlasmaParticleContainer.H

@@ -212,6 +212,8 @@ public:
     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;
+    /** to calculate laser ionization momentum width for linear polarization */
+     amrex::Gpu::DeviceVector<amrex::Real> m_laser_dp_prefactor;
 
     // plasma sorting/reordering:
 

src/particles/plasma/PlasmaParticleContainer.cpp

@@ -426,7 +426,7 @@ IonizationModule (const int lev,
                 const long pidx = pid + old_size;
 
                 // Copy ion data to new electron
-                amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // sets the ionized electron ID to 2 (valid/invalid) for the new electron
+                amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // sets the ionized electron ID to 2 (valid/invalid) for the ionized electrons
                 amrex::ParticleCPUWrapper{idcpu_elec[pidx]} = lev; // current level
                 arrdata_elec[PlasmaIdx::x      ][pidx] = arrdata_ion[PlasmaIdx::x     ][ip];
                 arrdata_elec[PlasmaIdx::y      ][pidx] = arrdata_ion[PlasmaIdx::y     ][ip];
@@ -531,6 +531,7 @@ LaserIonization (const int islice,
         amrex::Real* AMREX_RESTRICT adk_exp_prefactor = m_adk_exp_prefactor.data();
         amrex::Real* AMREX_RESTRICT adk_power = m_adk_power.data();
         amrex::Real* AMREX_RESTRICT laser_adk_prefactor = m_laser_adk_prefactor.data();
+        amrex::Real* AMREX_RESTRICT laser_dp_prefactor = m_laser_dp_prefactor.data();
 
         long num_ions = ptile_ion.numParticles();
 
@@ -623,30 +624,84 @@ LaserIonization (const int islice,
         amrex::Gpu::DeviceScalar<uint32_t> ip_elec(0);
         uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr();
 
-        // This kernel adds the new ionized electrons to the Plasma Particle Container
-        amrex::ParallelFor(num_ions,
-            [=] AMREX_GPU_DEVICE (long ip) {
+        // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time
+        // and calculates the momentum of the ionized electron based on equations (B8) and (B9)
+        // from the Massimo, 2020 article. It computes the energy of the emitted electron
+        // and assigns the resulting properties (momentum, position, etc.) to the new electrons
+        // created in the plasma container.
+        amrex::AnyCTO(
+            amrex::TypeList<
+                amrex::CompileTimeOptions<0, 1, 2, 3>
+            >{}, {
+                Hipace::m_depos_order_xy
+            },
+            [&] (auto cto_func) {
+                amrex::ParallelForRNG(num_ions, cto_func);
+            },
+            [=] AMREX_GPU_DEVICE (long ip, const amrex::RandomEngine& engine,
+                                  auto depos_order_xy) {
 
             if(p_ion_mask[ip] != 0) {
+
+                // Avoid temp slice
+                const amrex::Real xp = x_prev[ip];
+                const amrex::Real yp = y_prev[ip];
+
+                if (amrex::ConstParticleIDWrapper(idcpup[ip]) < 0 ||
+                    !laser_bounds.contains(xp, yp)) return;
+
+                Complex A = 0;
+                Complex A_dx = 0;
+                Complex A_dzeta = 0;
+
+                doLaserGatherShapeN<depos_order_xy>(xp, yp, A, A_dx, A_dzeta, laser_arr,
+                    dx_inv, dy_inv, dzeta_inv, x_pos_offset, y_pos_offset);
+
+                const Complex Et = I * A * omega0 + A_dzeta * phys_const.c; // transverse component
+                const Complex El = - A_dx * phys_const.c; // longitudinal component
+
+                amrex::Real Ep = std::sqrt( amrex::abs(Et*Et) + amrex::abs(El*El) );
+                Ep *= phys_const.m_e * phys_const.c / phys_const.q_e;
+                Ep *= E0;
+
+                amrex::Real ux;
+                amrex::Real uy;
+                amrex::Real uz;
+                const int ion_lev_loc = ion_lev[ip];
+
+                if (linear_polarization) {
+                    amrex::Real width_p;
+                    amrex::Real p_pol;
+                    width_p = std::sqrt(laser_dp_prefactor[ion_lev_loc] * Ep) * std::sqrt(amrex::abs(A*A)); // equation (4) art. Massimo
+                    p_pol = amrex::RandomNormal(0.0, width_p, engine);
+                    ux = p_pol;
+                    uy = 0._rt;
+                    uz = (amrex::abs(A * A) / 4. + p_pol * p_pol / 2.);
+                } else {
+                    amrex::Real angle;
+                    angle = amrex::Random(engine) * 2 * MathConst::pi;
+                    ux = std::sqrt(amrex::abs(A*A)) / std::sqrt(2) * std::cos(angle);
+                    uy = std::sqrt(amrex::abs(A*A)) / std::sqrt(2) * std::sin(angle);
+                    uz = amrex::abs(A*A) / 2.;
+                }
+
                 const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec`
                 const long pidx = pid + old_size;
-
                 // Copy ion data to new electron
-                amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // sets the ionized electron ID to 2 (valid/invalid) for the new electron
+                amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // sets the ionized electron ID to 2 (valid/invalid) for the ionized electrons
                 amrex::ParticleCPUWrapper{idcpu_elec[pidx]} =
                     amrex::ParticleCPUWrapper{idcpu_ion[pidx]}; // current level
                 arrdata_elec[PlasmaIdx::x      ][pidx] = arrdata_ion[PlasmaIdx::x     ][ip];
                 arrdata_elec[PlasmaIdx::y      ][pidx] = arrdata_ion[PlasmaIdx::y     ][ip];
-
                 arrdata_elec[PlasmaIdx::w      ][pidx] = arrdata_ion[PlasmaIdx::w     ][ip];
-                arrdata_elec[PlasmaIdx::ux     ][pidx] = 0._rt;
-                arrdata_elec[PlasmaIdx::uy     ][pidx] = 0._rt;
-                arrdata_elec[PlasmaIdx::psi    ][pidx] = 1._rt;
+                arrdata_elec[PlasmaIdx::ux     ][pidx] = ux * phys_const.c;
+                arrdata_elec[PlasmaIdx::uy     ][pidx] = uy * phys_const.c;
+                arrdata_elec[PlasmaIdx::psi    ][pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz)-uz; //psi = gamma - uz
                 arrdata_elec[PlasmaIdx::x_prev ][pidx] = arrdata_ion[PlasmaIdx::x_prev][ip];
                 arrdata_elec[PlasmaIdx::y_prev ][pidx] = arrdata_ion[PlasmaIdx::y_prev][ip];
-                arrdata_elec[PlasmaIdx::ux_half_step ][pidx] = 0._rt;
-                arrdata_elec[PlasmaIdx::uy_half_step ][pidx] = 0._rt;
-                arrdata_elec[PlasmaIdx::psi_half_step][pidx] = 1._rt;
+                arrdata_elec[PlasmaIdx::ux_half_step ][pidx] = ux * phys_const.c;
+                arrdata_elec[PlasmaIdx::uy_half_step ][pidx] = uy * phys_const.c;
+                arrdata_elec[PlasmaIdx::psi_half_step][pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz)-uz;
 #ifdef HIPACE_USE_AB5_PUSH
 #ifdef AMREX_USE_GPU
 #pragma unroll

src/particles/plasma/PlasmaParticleContainerInit.cpp

@@ -411,7 +411,7 @@ InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_
     // 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
+    // without Gamma function.
     const PhysConst phys_const = make_constants_SI();
     const amrex::Real alpha = 0.0072973525693_rt;
     const amrex::Real a3 = alpha * alpha * alpha;
@@ -431,11 +431,13 @@ InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_
     m_adk_prefactor.resize(ion_atomic_number);
     m_adk_exp_prefactor.resize(ion_atomic_number);
     m_laser_adk_prefactor.resize(ion_atomic_number);
+    m_laser_dp_prefactor.resize(ion_atomic_number);
 
     amrex::Gpu::PinnedVector<amrex::Real> h_adk_power(ion_atomic_number);
     amrex::Gpu::PinnedVector<amrex::Real> h_adk_prefactor(ion_atomic_number);
     amrex::Gpu::PinnedVector<amrex::Real> h_adk_exp_prefactor(ion_atomic_number);
     amrex::Gpu::PinnedVector<amrex::Real> h_laser_adk_prefactor(ion_atomic_number);
+    amrex::Gpu::PinnedVector<amrex::Real> h_laser_dp_prefactor(ion_atomic_number);
 
     for (int i=0; i<ion_atomic_number; ++i)
     {
@@ -448,6 +450,7 @@ InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_
             * std::pow(2*std::pow((Uion/UH),3./2.)*Ea,2*n_eff - 1);
         h_adk_exp_prefactor[i] = -2./3. * std::pow( Uion/UH,3./2.) * Ea;
         h_laser_adk_prefactor[i] = (3./MathConst::pi) * std::pow(Uion/UH, -3./2.) / Ea;
+    h_laser_dp_prefactor[i] = 3./2. * std::pow(Uion/UH, -3./2.) / Ea;
     }
 
     amrex::Gpu::copy(amrex::Gpu::hostToDevice,
@@ -458,4 +461,6 @@ InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_
         h_adk_exp_prefactor.begin(), h_adk_exp_prefactor.end(), m_adk_exp_prefactor.begin());
     amrex::Gpu::copy(amrex::Gpu::hostToDevice,
         h_laser_adk_prefactor.begin(), h_laser_adk_prefactor.end(), m_laser_adk_prefactor.begin());
+    amrex::Gpu::copy(amrex::Gpu::hostToDevice,
+         h_laser_dp_prefactor.begin(), h_laser_dp_prefactor.end(), m_laser_dp_prefactor.begin());
 }

tests/checksum/benchmarks_json/laser_ionization.1Rank.json

@@ -1,24 +1,24 @@
 {
   "lev=0": {
-    "Bx": 2.6828383569303,
-    "By": 2.7679232216354,
-    "Bz": 1.8909285687383,
+    "Bx": 3.6289326913652,
+    "By": 3.8848759208256,
+    "Bz": 1.6623540799087,
     "ExmBy": 0.0,
     "EypBx": 0.0,
-    "Ez": 2448530507.7321,
+    "Ez": 2705084383.7236,
     "Psi": 0.0,
-    "Sx": 174799612003.59,
-    "Sy": 165714923911.0,
+    "Sx": 176801264949.4,
+    "Sy": 163740521603.12,
     "aabs": 0.19201022247138,
-    "chi": 105640532321080.0,
-    "jx": 340889782577.97,
+    "chi": 104705826169570.0,
+    "jx": 455020279202.25,
     "jx_beam": 0.0,
-    "jy": 344102223375.08,
+    "jy": 431292339864.21,
     "jy_beam": 0.0,
     "jz_beam": 0.0,
     "laserEnvelope": 30.678508193173,
-    "rho_elec": 477708750.68105,
-    "rho_ion": 477712998.48921,
+    "rho_elec": 473481992.40286,
+    "rho_ion": 473479247.80813,
     "rhomjz": 0.0
   }
 }