Skip to content

Commit

Permalink
add doc fix indent, assert AB5
Browse files Browse the repository at this point in the history
  • Loading branch information
@SeverinDiederichs
SeverinDiederichs committed Aug 16, 2023
1 parent a2616f8 commit 8028fd7
Show file tree
Hide file tree
Showing 3 changed files with 158 additions and 147 deletions.
4 changes: 4 additions & 0 deletions docs/source/run/parameters.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -376,6 +376,10 @@ When both are specified, the per-species value is used.
A good starting point is a period of 4 to reorder plasma particles on every fourth zeta-slice.
To disable reordering set this to 0.

* ``<plasma name> or plasmas.n_subcycles`` (`int`) optional (default `1`)
Number of sub-cycles within the plasma pusher. Currently only implemented for the leapfrog pusher. Must be larger or equal to 1. Sub-cycling is needed if plasma particles move
significantly in the transverse direction during a single longitudinal cell. If they move too many cells such that they do not sample certain small transverse structures in the wakefields, sub-cycling is needed and fixes the issue.

* ``<plasma name> or plasmas.reorder_idx_type`` (2 `int`) optional (default `0 0` or `1 1`)
Change if plasma particles are binned to cells (0), nodes (1) or both (2)
for both x and y direction as part of the reordering.
Expand Down
7 changes: 7 additions & 0 deletions src/particles/plasma/PlasmaParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,13 @@ PlasmaParticleContainer::ReadParameters ()
}

queryWithParserAlt(pp, "n_subcycles", m_n_subcycles, pp_alt);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_n_subcycles >= 1,
"n_subcycles must be larger or equal to 1 sub-cycle (default is 1)");
#ifdef HIPACE_USE_AB5_PUSH
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_n_subcycles == 1,
"Plasma subcycling only implemeted for leapfrog pusher!"
"Please set plasmas.n_subcycles = 1");
#endif
queryWithParser(pp, "mass_Da", mass_Da);
if(mass_Da != 0) {
m_mass = phys_const.m_p * mass_Da / 1.007276466621;
Expand Down
294 changes: 147 additions & 147 deletions src/particles/pusher/PlasmaParticleAdvance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,162 +107,162 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields,
amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev)[ip];

for (int i = 0; i < n_subcycles; i++) {
// define field at particle position reals
amrex::Real ExmByp = 0._rt, EypBxp = 0._rt, Ezp = 0._rt;
amrex::Real Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt;
amrex::Real Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt;

doGatherShapeN<depos_order.value>(xp, yp, ExmByp, EypBxp, Ezp, Bxp, Byp,
Bzp, slice_arr, psi_comp, ez_comp, bx_comp, by_comp,
bz_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset);

if (use_laser.value) {
doLaserGatherShapeN<depos_order.value>(xp, yp, Aabssqp, AabssqDxp,
AabssqDyp, a_arr, dx_inv, dy_inv, x_pos_offset, y_pos_offset);
}

amrex::Real q_mass_clight_ratio = charge_mass_clight_ratio;
if (can_ionize) {
q_mass_clight_ratio *= ptd.idata(PlasmaIdx::ion_lev)[ip];
}
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;
// define field at particle position reals
amrex::Real ExmByp = 0._rt, EypBxp = 0._rt, Ezp = 0._rt;
amrex::Real Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt;
amrex::Real Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt;

doGatherShapeN<depos_order.value>(xp, yp, ExmByp, EypBxp, Ezp, Bxp, Byp,
Bzp, slice_arr, psi_comp, ez_comp, bx_comp, by_comp,
bz_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset);

if (use_laser.value) {
doLaserGatherShapeN<depos_order.value>(xp, yp, Aabssqp, AabssqDxp,
AabssqDyp, a_arr, dx_inv, dy_inv, x_pos_offset, y_pos_offset);
}

amrex::Real q_mass_clight_ratio = charge_mass_clight_ratio;
if (can_ionize) {
q_mass_clight_ratio *= ptd.idata(PlasmaIdx::ion_lev)[ip];
}
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;

#ifndef HIPACE_USE_AB5_PUSH

constexpr int nsub = 4;
const amrex::Real sdz = dz/nsub;

amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip];
amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip];
amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip];

// full push in momentum
// from t-1/2 to t+1/2
// using the fields at t
for (int isub=0; isub<nsub; ++isub) {

const amrex::Real psi_inv = 1._rt/psi;

auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush(
ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

const DualNumber ux_dual{ux, dz_ux};
const DualNumber uy_dual{uy, dz_uy};
const DualNumber psi_inv_dual{psi_inv, -psi_inv*psi_inv*dz_psi};

auto [dz_ux_dual, dz_uy_dual, dz_psi_dual] = PlasmaMomentumPush(
ux_dual, uy_dual, psi_inv_dual, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

ux += sdz*dz_ux + 0.5_rt*sdz*sdz*dz_ux_dual.epsilon;
uy += sdz*dz_uy + 0.5_rt*sdz*sdz*dz_uy_dual.epsilon;
psi += sdz*dz_psi + 0.5_rt*sdz*sdz*dz_psi_dual.epsilon;

}

// full push in position
// from t to t+1
// using the momentum at t+1/2
xp += dz*clight_inv*(ux * (1._rt/psi));
yp += dz*clight_inv*(uy * (1._rt/psi));

if (setPositionEnforceBC(ptd, ip, xp, yp)) return;

if (!temp_slice) {
// update values of the last non temp slice
// the next push always starts from these
ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux;
ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy;
ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi;
ptd.rdata(PlasmaIdx::x_prev)[ip] = xp;
ptd.rdata(PlasmaIdx::y_prev)[ip] = yp;
}

// half push in momentum
// from t+1/2 to t+1
// still using the fields at t as an approximation
// the result is used for current deposition etc. but not in the pusher
for (int isub=0; isub<(nsub/2); ++isub) {

const amrex::Real psi_inv = 1._rt/psi;

auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush(
ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

const DualNumber ux_dual{ux, dz_ux};
const DualNumber uy_dual{uy, dz_uy};
const DualNumber psi_inv_dual{psi_inv, -psi_inv*psi_inv*dz_psi};

auto [dz_ux_dual, dz_uy_dual, dz_psi_dual] = PlasmaMomentumPush(
ux_dual, uy_dual, psi_inv_dual, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

ux += sdz*dz_ux + 0.5_rt*sdz*sdz*dz_ux_dual.epsilon;
uy += sdz*dz_uy + 0.5_rt*sdz*sdz*dz_uy_dual.epsilon;
psi += sdz*dz_psi + 0.5_rt*sdz*sdz*dz_psi_dual.epsilon;

}
ptd.rdata(PlasmaIdx::ux)[ip] = ux;
ptd.rdata(PlasmaIdx::uy)[ip] = uy;
ptd.rdata(PlasmaIdx::psi)[ip] = psi;
} // for (int isub=0; isub<nsub; ++isub) {
constexpr int nsub = 4;
const amrex::Real sdz = dz/nsub;

amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip];
amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip];
amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip];

// full push in momentum
// from t-1/2 to t+1/2
// using the fields at t
for (int isub=0; isub<nsub; ++isub) {

const amrex::Real psi_inv = 1._rt/psi;

auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush(
ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

const DualNumber ux_dual{ux, dz_ux};
const DualNumber uy_dual{uy, dz_uy};
const DualNumber psi_inv_dual{psi_inv, -psi_inv*psi_inv*dz_psi};

auto [dz_ux_dual, dz_uy_dual, dz_psi_dual] = PlasmaMomentumPush(
ux_dual, uy_dual, psi_inv_dual, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

ux += sdz*dz_ux + 0.5_rt*sdz*sdz*dz_ux_dual.epsilon;
uy += sdz*dz_uy + 0.5_rt*sdz*sdz*dz_uy_dual.epsilon;
psi += sdz*dz_psi + 0.5_rt*sdz*sdz*dz_psi_dual.epsilon;

}

// full push in position
// from t to t+1
// using the momentum at t+1/2
xp += dz*clight_inv*(ux * (1._rt/psi));
yp += dz*clight_inv*(uy * (1._rt/psi));

if (setPositionEnforceBC(ptd, ip, xp, yp)) return;

if (!temp_slice) {
// update values of the last non temp slice
// the next push always starts from these
ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux;
ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy;
ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi;
ptd.rdata(PlasmaIdx::x_prev)[ip] = xp;
ptd.rdata(PlasmaIdx::y_prev)[ip] = yp;
}

// half push in momentum
// from t+1/2 to t+1
// still using the fields at t as an approximation
// the result is used for current deposition etc. but not in the pusher
for (int isub=0; isub<(nsub/2); ++isub) {

const amrex::Real psi_inv = 1._rt/psi;

auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush(
ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

const DualNumber ux_dual{ux, dz_ux};
const DualNumber uy_dual{uy, dz_uy};
const DualNumber psi_inv_dual{psi_inv, -psi_inv*psi_inv*dz_psi};

auto [dz_ux_dual, dz_uy_dual, dz_psi_dual] = PlasmaMomentumPush(
ux_dual, uy_dual, psi_inv_dual, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

ux += sdz*dz_ux + 0.5_rt*sdz*sdz*dz_ux_dual.epsilon;
uy += sdz*dz_uy + 0.5_rt*sdz*sdz*dz_uy_dual.epsilon;
psi += sdz*dz_psi + 0.5_rt*sdz*sdz*dz_psi_dual.epsilon;

}
ptd.rdata(PlasmaIdx::ux)[ip] = ux;
ptd.rdata(PlasmaIdx::uy)[ip] = uy;
ptd.rdata(PlasmaIdx::psi)[ip] = psi;
#else
amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip];
amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip];
amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip];
const amrex::Real psi_inv = 1._rt/psi;

auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush(
ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

ptd.rdata(PlasmaIdx::Fx1)[ip] = clight_inv*(ux * psi_inv);
ptd.rdata(PlasmaIdx::Fy1)[ip] = clight_inv*(uy * psi_inv);
ptd.rdata(PlasmaIdx::Fux1)[ip] = dz_ux;
ptd.rdata(PlasmaIdx::Fuy1)[ip] = dz_uy;
ptd.rdata(PlasmaIdx::Fpsi1)[ip] = dz_psi;

const amrex::Real ab5_coeffs[5] = {
( 1901._rt / 720._rt ) * dz, // a1 times dz
( -1387._rt / 360._rt ) * dz, // a2 times dz
( 109._rt / 30._rt ) * dz, // a3 times dz
( -637._rt / 360._rt ) * dz, // a4 times dz
( 251._rt / 720._rt ) * dz // a5 times dz
};
amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip];
amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip];
amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip];
const amrex::Real psi_inv = 1._rt/psi;

auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush(
ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp,
Aabssqp, AabssqDxp, AabssqDyp, clight_inv, q_mass_clight_ratio);

ptd.rdata(PlasmaIdx::Fx1)[ip] = clight_inv*(ux * psi_inv);
ptd.rdata(PlasmaIdx::Fy1)[ip] = clight_inv*(uy * psi_inv);
ptd.rdata(PlasmaIdx::Fux1)[ip] = dz_ux;
ptd.rdata(PlasmaIdx::Fuy1)[ip] = dz_uy;
ptd.rdata(PlasmaIdx::Fpsi1)[ip] = dz_psi;

const amrex::Real ab5_coeffs[5] = {
( 1901._rt / 720._rt ) * dz, // a1 times dz
( -1387._rt / 360._rt ) * dz, // a2 times dz
( 109._rt / 30._rt ) * dz, // a3 times dz
( -637._rt / 360._rt ) * dz, // a4 times dz
( 251._rt / 720._rt ) * dz // a5 times dz
};

#ifdef AMREX_USE_GPU
#pragma unroll
#endif
for (int iab=0; iab<5; ++iab) {
xp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fx1 + iab)[ip];
yp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fy1 + iab)[ip];
ux += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fux1 + iab)[ip];
uy += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fuy1 + iab)[ip];
psi += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fpsi1 + iab)[ip];
}

if (setPositionEnforceBC(ptd, ip, xp, yp)) return;

if (!temp_slice) {
// update values of the last non temp slice
// the next push always starts from these
ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux;
ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy;
ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi;
ptd.rdata(PlasmaIdx::x_prev)[ip] = xp;
ptd.rdata(PlasmaIdx::y_prev)[ip] = yp;
}

ptd.rdata(PlasmaIdx::ux)[ip] = ux;
ptd.rdata(PlasmaIdx::uy)[ip] = uy;
ptd.rdata(PlasmaIdx::psi)[ip] = psi;
for (int iab=0; iab<5; ++iab) {
xp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fx1 + iab)[ip];
yp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fy1 + iab)[ip];
ux += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fux1 + iab)[ip];
uy += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fuy1 + iab)[ip];
psi += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fpsi1 + iab)[ip];
}

if (setPositionEnforceBC(ptd, ip, xp, yp)) return;

if (!temp_slice) {
// update values of the last non temp slice
// the next push always starts from these
ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux;
ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy;
ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi;
ptd.rdata(PlasmaIdx::x_prev)[ip] = xp;
ptd.rdata(PlasmaIdx::y_prev)[ip] = yp;
}

ptd.rdata(PlasmaIdx::ux)[ip] = ux;
ptd.rdata(PlasmaIdx::uy)[ip] = uy;
ptd.rdata(PlasmaIdx::psi)[ip] = psi;
#endif
} // end loop over subcycles
});
}

Expand Down

0 comments on commit 8028fd7

Please sign in to comment.