hamiltonian method

Author: EyaDammak

src/laser/MultiLaser.H

@@ -190,7 +190,6 @@ public:
      */
     bool UseLaser (const int islice) const { return m_use_laser && HasSlice(islice); }
 
-private:
 
     bool m_use_laser {false}; /**< whether a laser is used or not */
     /** Laser central wavelength defined by user.

src/particles/plasma/PlasmaParticleContainer.H

@@ -32,6 +32,7 @@ struct PlasmaIdx
         ux_half_step,       // momentum half a step behind the current slice for leapfrog pusher
         uy_half_step,       // at the same step for AB5 pusher
         psi_half_step,      // never effected by temp slice
+        phi_ion,
 #ifdef HIPACE_USE_AB5_PUSH
         Fx1, Fx2, Fx3, Fx4, Fx5,             // AB5 force terms
         Fy1, Fy2, Fy3, Fy4, Fy5,             //
@@ -239,6 +240,7 @@ public:
     bool m_can_laser_injection = false; /**< whether laser injection is enabled */
     amrex::Real m_uz_threshold = 1.;
     amrex::Real m_injection_weight_factor = 1.;
+    amrex::Real m_beta_p = 1.;
     std::string m_product_beam_name = ""; /**< name of Injection product beam */
     BeamParticleContainer* m_product_beam_pc = nullptr; /**< Injection product beam */
 

src/particles/plasma/PlasmaParticleContainer.cpp

@@ -77,6 +77,7 @@ PlasmaParticleContainer::ReadParameters ()
     queryWithParser(pp, "can_laser_injection", m_can_laser_injection);
     queryWithParser(pp, "uz_threshold", m_uz_threshold);
     queryWithParser(pp, "injection_weight_factor", m_injection_weight_factor);
+    queryWithParser(pp, "beta_p", m_beta_p);
 
     m_can_ionize = m_can_field_ionize || m_can_laser_ionize;
 
@@ -715,6 +716,7 @@ LaserIonization (const int islice,
                 arrdata_elec[PlasmaIdx::w      ][pidx] = arrdata_ion[PlasmaIdx::w     ][ip];
                 arrdata_elec[PlasmaIdx::ux     ][pidx] = ux * phys_const.c;
                 arrdata_elec[PlasmaIdx::uy     ][pidx] = uy * phys_const.c;
+                arrdata_elec[PlasmaIdx::phi_ion][pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz + 0.5_rt*amrex::abs(A*A)) - m_beta_p * uz - 1._rt;
                 arrdata_elec[PlasmaIdx::psi    ][pidx] = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz + 0.5_rt*amrex::abs(A*A))-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];
@@ -768,6 +770,10 @@ InjectionCondition (const int lev, const Fields& fields, const MultiLaser& laser
         const int ez_comp = Comps[WhichSlice::This]["Ez"];
         auto laser_geom = laser.GetLaserGeom();
 
+        amrex::Real a0 = laser.m_all_lasers[0].m_a0;
+        amrex::Real e_max = a0*a0/2*std::pow((1+a0*a0/2),-0.5_rt);
+        amrex::Real phi_min = e_max*e_max/2 - m_beta_p*std::sqrt((1+e_max*e_max/2)*(1+e_max*e_max/2)-1);
+
         // This kernel marks the plasma particles that has been injected in the wake
         amrex::ParallelFor(num_particles,
             [=] AMREX_GPU_DEVICE (int ip) {
@@ -804,7 +810,12 @@ InjectionCondition (const int lev, const Fields& fields, const MultiLaser& laser
 
                 amrex::Real condition = uz - uz_condition; // condition for injection
 
-                if (condition > 0 && Ezp < 0){
+                amrex::Real gam_p = std::sqrt(1._rt + ux*ux + uy*uy + uz*uz + 0.5_rt*amrex::abs(A*A));
+                amrex::Real gam_p_inv = 1/gam_p;
+                amrex::Real phi_ion = ptd_plasma.rdata(PlasmaIdx::phi_ion)[ip];
+                amrex::Real condition2 = phi_ion - phi_min + gam_p_inv - 1._rt;
+
+                if (condition2 > 0){
                     ptd_plasma.id(ip) = 3; // set the injected electron ID to 3
                 }
         });
@@ -925,7 +936,7 @@ PlasmaToBeam (const MultiLaser& laser, amrex::Vector<amrex::Geometry> const& gm,
                     amrex::Real psi = ptd_plasma.rdata(PlasmaIdx::psi)[ip];
                     ptd_beam.rdata(BeamIdx::uz)[pidx_beam] = (1+ux*ux+uy*uy - psi*psi + 0.5_rt*amrex::abs(A*A))/(2.*psi)*phys_const.c;
                     amrex::Real uz = ptd_beam.rdata(BeamIdx::uz)[pidx_beam] * clight_inv;
-                    const amrex::Real gam = std::sqrt(1. + ux*ux + uy*uy + uz*uz);
+                    const amrex::Real gam = std::sqrt(1. + ux*ux + uy*uy + uz*uz + 0.5_rt*amrex::abs(A*A));
                     ptd_beam.rdata(BeamIdx::w)[pidx_beam] = ptd_plasma.rdata(PlasmaIdx::w)[ip] * gam / (psi) * f;
                     // conservation of j_x and j_y
                     ptd_beam.idata(BeamIdx::nsubcycles)[pidx_beam] = 0;