diff --git a/src/Hipace.cpp b/src/Hipace.cpp index b2344c3a3b..6376401637 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -335,7 +335,6 @@ Hipace::Evolve () if (m_multi_laser.m_use_laser) { AMREX_ALWAYS_ASSERT(!m_adaptive_time_step.m_do_adaptive_time_step); - AMREX_ALWAYS_ASSERT(m_multi_plasma.GetNPlasmas() <= 1); } m_physical_time = step == 0 ? m_initial_time : m_multi_buffer.get_time(); diff --git a/src/particles/deposition/ExplicitDeposition.cpp b/src/particles/deposition/ExplicitDeposition.cpp index 1ea4209fbe..e15e2c06e0 100644 --- a/src/particles/deposition/ExplicitDeposition.cpp +++ b/src/particles/deposition/ExplicitDeposition.cpp @@ -126,16 +126,18 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Multi amrex::Real Aabssqp = 0._rt; // Rename variable for NVCC lambda capture to work [[maybe_unused]] auto laser_arr = a_laser_arr; + [[maybe_unused]] auto laser_fac = a_laser_fac; if constexpr (use_laser.value) { // Its important that Aabssqp is first fully gathered and not used // directly per cell like AabssqDxp and AabssqDyp doLaserGatherShapeN(xp, yp, Aabssqp, laser_arr, dx_inv, dy_inv, x_pos_offset, y_pos_offset); + Aabssqp *= laser_fac * q_mass_ratio * q_mass_ratio; } // calculate gamma/psi for plasma particles const amrex::Real gamma_psi = 0.5_rt * ( - (1._rt + 0.5_rt * Aabssqp) * psi_inv * psi_inv // TODO: fix units + (1._rt + 0.5_rt * Aabssqp) * psi_inv * psi_inv + vx * vx + vy * vy + 1._rt @@ -171,7 +173,6 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Multi amrex::Real AabssqDxp = 0._rt; amrex::Real AabssqDyp = 0._rt; // Rename variables for NVCC lambda capture to work - [[maybe_unused]] auto laser_fac = a_laser_fac; [[maybe_unused]] auto clight = a_clight; if constexpr (use_laser.value) { // avoid going outside of domain diff --git a/src/particles/deposition/PlasmaDepositCurrent.cpp b/src/particles/deposition/PlasmaDepositCurrent.cpp index 6dd1f27c16..a5cd99e178 100644 --- a/src/particles/deposition/PlasmaDepositCurrent.cpp +++ b/src/particles/deposition/PlasmaDepositCurrent.cpp @@ -90,6 +90,8 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas const amrex::Real clightinv = 1.0_rt/pc.c; const amrex::Real charge_invvol = charge * invvol; const amrex::Real charge_mu0_mass_ratio = charge * pc.mu0 / mass; + const amrex::Real laser_norm = (charge/pc.q_e) * (pc.m_e/mass) + * (charge/pc.q_e) * (pc.m_e/mass); int n_qsa_violation = 0; amrex::Gpu::DeviceScalar gpu_n_qsa_violation(n_qsa_violation); @@ -199,9 +201,12 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas // calculate charge of the plasma particles amrex::Real q_invvol = charge_invvol * ptd.rdata(PlasmaIdx::w)[ip]; amrex::Real q_mu0_mass_ratio = charge_mu0_mass_ratio; + [[maybe_unused]] amrex::Real laser_norm_ion = laser_norm; if constexpr (can_ionize.value) { q_invvol *= ptd.idata(PlasmaIdx::ion_lev)[ip]; q_mu0_mass_ratio *= ptd.idata(PlasmaIdx::ion_lev)[ip]; + laser_norm_ion *= + ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; } const amrex::Real xmid = (xp - x_pos_offset) * dx_inv; @@ -212,11 +217,12 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas if constexpr (use_laser.value) { doLaserGatherShapeN(xp, yp, Aabssqp, laser_arr, dx_inv, dy_inv, x_pos_offset, y_pos_offset); + Aabssqp *= laser_norm_ion; } // calculate gamma/psi for plasma particles const amrex::Real gamma_psi = 0.5_rt * ( - (1._rt + 0.5_rt * Aabssqp) * psi_inv * psi_inv // TODO: fix units + (1._rt + 0.5_rt * Aabssqp) * psi_inv * psi_inv + vx_c * vx_c * clightinv * clightinv + vy_c * vy_c * clightinv * clightinv + 1._rt diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 7baa63fc93..a3b7501ec0 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -77,7 +77,8 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, const auto setPositionEnforceBC = EnforceBCandSetPos(gm[0]); const amrex::Real dz = gm[0].CellSize(2) / n_subcycles; - const amrex::Real me_clight_mass_ratio = phys_const.c * phys_const.m_e/plasma.m_mass; + const amrex::Real laser_norm = (plasma.m_charge/phys_const.q_e) * (phys_const.m_e/plasma.m_mass) + * (plasma.m_charge/phys_const.q_e) * (phys_const.m_e/plasma.m_mass); const amrex::Real clight = phys_const.c; const amrex::Real clight_inv = 1._rt/phys_const.c; const amrex::Real charge_mass_clight_ratio = plasma.m_charge/(plasma.m_mass * phys_const.c); @@ -114,8 +115,11 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, amrex::Real Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt; amrex::Real q_mass_clight_ratio = charge_mass_clight_ratio; + amrex::Real laser_norm_ion = laser_norm; if (can_ionize) { q_mass_clight_ratio *= ptd.idata(PlasmaIdx::ion_lev)[ip]; + laser_norm_ion *= + ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; } for (int i = 0; i < n_subcycles; i++) { @@ -141,9 +145,9 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, Bxp *= clight; Byp *= clight; - Aabssqp *= 0.5_rt; // TODO: fix units of Aabssqp - AabssqDxp *= 0.25_rt * me_clight_mass_ratio; - AabssqDyp *= 0.25_rt * me_clight_mass_ratio; + Aabssqp *= 0.5_rt * laser_norm_ion; + AabssqDxp *= 0.25_rt * clight * laser_norm_ion; + AabssqDyp *= 0.25_rt * clight * laser_norm_ion; } #ifndef HIPACE_USE_AB5_PUSH diff --git a/src/particles/pusher/PushPlasmaParticles.H b/src/particles/pusher/PushPlasmaParticles.H index 1dd05ccabb..f035c04c0e 100644 --- a/src/particles/pusher/PushPlasmaParticles.H +++ b/src/particles/pusher/PushPlasmaParticles.H @@ -29,9 +29,9 @@ struct PlasmaMomentumDerivative { * \param[in] Bx_clight Bx * clight * \param[in] By_clight By * clight * \param[in] Bz Bz - * \param[in] Aabssq_norm Aabssqp * 0.5 TODO: fix units of Aabssqp - * \param[in] AabssqDx_norm AabssqDx * 0.25 * clight * m_e / mass - * \param[in] AabssqDy_norm AabssqDy * 0.25 * clight * m_e / mass + * \param[in] Aabssq_norm Aabssqp * 0.5 * (charge / q_e)^2 * (m_e / mass)^2 + * \param[in] AabssqDx_norm AabssqDx * 0.25 * clight * (charge / q_e)^2 * (m_e / mass)^2 + * \param[in] AabssqDy_norm AabssqDy * 0.25 * clight * (charge / q_e)^2 * (m_e / mass)^2 * \param[in] clight_inv 1 / clight * \param[in] charge_mass_clight_ratio charge / (mass * clight) */