Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Laser with non-electron plasma #1002

Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion src/Hipace.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
5 changes: 3 additions & 2 deletions src/particles/deposition/ExplicitDeposition.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<depos_order>(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
Expand Down Expand Up @@ -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
Expand Down
8 changes: 7 additions & 1 deletion src/particles/deposition/PlasmaDepositCurrent.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<int> gpu_n_qsa_violation(n_qsa_violation);
Expand Down Expand Up @@ -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;
Expand All @@ -212,11 +217,12 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas
if constexpr (use_laser.value) {
doLaserGatherShapeN<depos_order_xy>(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
Expand Down
12 changes: 8 additions & 4 deletions src/particles/pusher/PlasmaParticleAdvance.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,8 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields,
const auto setPositionEnforceBC = EnforceBCandSetPos<PTileType>(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);
Expand Down Expand Up @@ -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++) {
Expand All @@ -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
Expand Down
6 changes: 3 additions & 3 deletions src/particles/pusher/PushPlasmaParticles.H
Original file line number Diff line number Diff line change
Expand Up @@ -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)
*/
Expand Down
Loading