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.

Sign up for GitHub

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

Add derivative shape factors #831

Merged
merged 2 commits into from
Dec 2, 2022
Merged

Add derivative shape factors #831

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
26026c3
add derivative shape factor
AlexanderSinn Nov 29, 2022
98a5e0f
fix factor of -1
AlexanderSinn Nov 29, 2022
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
3 changes: 3 additions & 0 deletions docs/source/run/parameters.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,9 @@ General parameters
* ``hipace.depos_order_z`` (`int`) optional (default `0`)
Transverse particle shape order. Currently, only `0` is implemented.

* ``hipace.depos_derivative_type`` (`int`) optional (default `2`)
Type of derivative used in explicit deposition. `0`: analytic, `1`: nodal, `2`: centered

* ``hipace.outer_depos_loop`` (`bool`) optional (default `0`)
If the loop over depos_order is included in the loop over particles.

Expand Down
3 changes: 3 additions & 0 deletions src/Hipace.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -313,6 +313,9 @@ public:
/** Order of the field gather and current deposition shape factor in the longitudinal direction
*/
static int m_depos_order_z;
/** Type of derivative used in explicit deposition. 0: analytic, 1: nodal, 2: centered
*/
static int m_depos_derivative_type;
/* if the loop over depos_order is included in the loop over particles
*/
static bool m_outer_depos_loop;
Expand Down
4 changes: 4 additions & 0 deletions src/Hipace.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ amrex::Real Hipace::m_initial_time = 0.0;
int Hipace::m_verbose = 0;
int Hipace::m_depos_order_xy = 2;
int Hipace::m_depos_order_z = 0;
int Hipace::m_depos_derivative_type = 2;
bool Hipace::m_outer_depos_loop = false;
amrex::Real Hipace::m_predcorr_B_error_tolerance = 4e-2;
int Hipace::m_predcorr_max_iterations = 30;
Expand Down Expand Up @@ -128,6 +129,9 @@ Hipace::Hipace () :
"Multiple boxes per rank only implemented for one rank.");
queryWithParser(pph, "depos_order_xy", m_depos_order_xy);
queryWithParser(pph, "depos_order_z", m_depos_order_z);
queryWithParser(pph, "depos_derivative_type", m_depos_derivative_type);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_depos_order_xy != 0 || m_depos_derivative_type != 0,
"Analytic derivative with depos_order=0 would vanish");
queryWithParser(pph, "outer_depos_loop", m_outer_depos_loop);
queryWithParser(pph, "predcorr_B_error_tolerance", m_predcorr_B_error_tolerance);
queryWithParser(pph, "predcorr_max_iterations", m_predcorr_max_iterations);
Expand Down
139 changes: 60 additions & 79 deletions src/particles/deposition/ExplicitDeposition.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@

#include "Hipace.H"
#include "particles/particles_utils/ShapeFactors.H"
#include "particles/particles_utils/FieldGather.H"
#include "utils/Constants.H"
#include "utils/HipaceProfilerWrapper.H"
#include "utils/GPUUtil.H"
Expand Down Expand Up @@ -47,10 +48,10 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Multi
int const * const AMREX_RESTRICT a_ion_lev =
plasma.m_can_ionize ? soa.GetIntData(PlasmaIdx::ion_lev).data() : nullptr;

// Construct empty Array4 with one z slice so that Array3 constructor works for no laser
// Construct empty Array4 if there is no laser
const Array3<const amrex::Real> a_laser_arr = multi_laser.m_use_laser ?
multi_laser.getSlices(WhichLaserSlice::n00j00).const_array(pti) :
amrex::Array4<const amrex::Real>(nullptr, {0,0,0}, {0,0,1}, 0);
amrex::Array4<const amrex::Real>();

const amrex::Real x_pos_offset = GetPosOffset(0, gm, isl_fab.box());
const amrex::Real y_pos_offset = GetPosOffset(1, gm, isl_fab.box());
Expand All @@ -62,27 +63,37 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Multi

const PhysConst pc = get_phys_const();
const amrex::Real clight = pc.c;
const amrex::Real clight_inv = 1._rt/pc.c;
// The laser a0 is always normalized
const amrex::Real a_laser_fac = (pc.m_e/pc.q_e) * (pc.m_e/pc.q_e);
const amrex::Real charge_invvol_mu0 = plasma.m_charge * invvol * pc.mu0;
const amrex::Real charge_mass_ratio = plasma.m_charge / plasma.m_mass;

amrex::ParallelFor(
amrex::TypeList<amrex::CompileTimeOptions<0, 1, 2, 3>, // depos_order
amrex::CompileTimeOptions<false, true>, // can_ionize
amrex::CompileTimeOptions<false, true>>{}, // use_laser
{Hipace::m_depos_order_xy, plasma.m_can_ionize, multi_laser.m_use_laser},
amrex::TypeList<
amrex::CompileTimeOptions<0, 1, 2, 3>, // depos_order
amrex::CompileTimeOptions<0, 1, 2>, // derivative_type
amrex::CompileTimeOptions<false, true>, // can_ionize
amrex::CompileTimeOptions<false, true> // use_laser
>{}, {
Hipace::m_depos_order_xy,
Hipace::m_depos_derivative_type,
plasma.m_can_ionize,
multi_laser.m_use_laser
},
pti.numParticles(),
[=] AMREX_GPU_DEVICE (int ip, auto a_depos_order, auto can_ionize, auto use_laser) {
[=] AMREX_GPU_DEVICE (int ip, auto a_depos_order, auto a_derivative_type,
auto can_ionize, auto use_laser) noexcept {
constexpr int depos_order = a_depos_order.value;
constexpr int derivative_type = a_derivative_type.value;

const auto positions = pos_structs[ip];
if (positions.id() < 0) return;
const amrex::Real psi = psip[ip];
const amrex::Real psi_inv = 1._rt/psip[ip];
const amrex::Real xp = positions.pos(0);
const amrex::Real yp = positions.pos(1);
const amrex::Real vx = uxp[ip] / (psi * clight);
const amrex::Real vy = uyp[ip] / (psi * clight);
const amrex::Real vx = uxp[ip] * psi_inv * clight_inv;
const amrex::Real vy = uyp[ip] * psi_inv * clight_inv;

amrex::Real q_invvol_mu0 = charge_invvol_mu0;
amrex::Real q_mass_ratio = charge_mass_ratio;
Expand All @@ -96,75 +107,46 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Multi
const amrex::Real charge_density_mu0 = q_invvol_mu0 * wp[ip];

const amrex::Real xmid = (xp - x_pos_offset) * dx_inv;
amrex::Real sx_cell[depos_order + 1];
const int i_cell = compute_shape_factor<depos_order>(sx_cell, xmid);

// y direction
const amrex::Real ymid = (yp - y_pos_offset) * dy_inv;
amrex::Real sy_cell[depos_order + 1];
const int j_cell = compute_shape_factor<depos_order>(sy_cell, ymid);

amrex::Real Aabssqp = 0._rt;
// Rename variable for NVCC lambda capture to work
[[maybe_unused]] auto laser_arr = a_laser_arr;
if constexpr (use_laser.value) {
for (int iy=0; iy<=depos_order; iy++){
for (int ix=0; ix<=depos_order; ix++){
// Its important that Aabssqp is first fully gathered and not used
// directly per cell like AabssqDxp and AabssqDyp
const amrex::Real x00y00 = abssq(
laser_arr(i_cell+ix, j_cell+iy, 0),
laser_arr(i_cell+ix, j_cell+iy, 1) );
Aabssqp += sx_cell[ix]*sy_cell[iy]*x00y00;
}
}
// 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);
}

// calculate gamma/psi for plasma particles
const amrex::Real gamma_psi = 0.5_rt * (
(1._rt + 0.5_rt * Aabssqp) / (psi * psi) // TODO: fix units
(1._rt + 0.5_rt * Aabssqp) * psi_inv * psi_inv // TODO: fix units
+ vx * vx
+ vy * vy
+ 1._rt
);

for (int iy=0; iy <= depos_order+2; ++iy) {
// normal shape factor
amrex::Real shape_y = 0._rt;
// derivative shape factor
amrex::Real shape_dy = 0._rt;
if (iy != 0 && iy != depos_order + 2) {
shape_y = sy_cell[iy-1] * charge_density_mu0;
}
if (iy < depos_order + 1) {
shape_dy = sy_cell[iy];
}
if (iy > 1) {
shape_dy -= sy_cell[iy-2];
}
shape_dy *= dy_inv * 0.5_rt * clight * charge_density_mu0;

for (int ix=0; ix <= depos_order+2; ++ix) {
amrex::Real shape_x = 0._rt;
amrex::Real shape_dx = 0._rt;
if (ix != 0 && ix != depos_order + 2) {
shape_x = sx_cell[ix-1];
}
if (ix < depos_order + 1) {
shape_dx = sx_cell[ix];
}
if (ix > 1) {
shape_dx -= sx_cell[ix-2];
}
shape_dx *= dx_inv * 0.5_rt * clight;

if ((ix==0 || ix==depos_order + 2) && (iy==0 || iy==depos_order + 2)) {
// corners have a shape factor of zero
continue;
#ifdef AMREX_USE_GPU
#pragma unroll
AlexanderSinn marked this conversation as resolved.
Show resolved Hide resolved
#endif
for (int iy=0; iy <= depos_order+derivative_type; ++iy) {
#ifdef AMREX_USE_GPU
#pragma unroll
AlexanderSinn marked this conversation as resolved.
Show resolved Hide resolved
#endif
for (int ix=0; ix <= depos_order+derivative_type; ++ix) {

if constexpr (derivative_type == 2) {
if ((ix==0 || ix==depos_order + 2) && (iy==0 || iy==depos_order + 2)) {
// corners have a shape factor of zero
continue;
}
}

const int i = i_cell + ix - 1;
const int j = j_cell + iy - 1;
auto [shape_y, shape_dy, j] =
single_derivative_shape_factor<derivative_type, depos_order>(ymid, iy);
auto [shape_x, shape_dx, i] =
single_derivative_shape_factor<derivative_type, depos_order>(xmid, ix);

// get fields per cell instead of gathering them to avoid blurring
const amrex::Real Bz_v = arr(i,j,Bz);
Expand Down Expand Up @@ -192,38 +174,37 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Multi
AabssqDyp = (x00yp1-x00ym1) * 0.5_rt * dy_inv * laser_fac * clight;
}

amrex::Gpu::Atomic::Add(arr.ptr(i, j, Sy),
amrex::Gpu::Atomic::Add(arr.ptr(i, j, Sy), charge_density_mu0 * (
- shape_x * shape_y * (
- Bz_v * vx
+ ( Ez_v * vy
+ ExmBy_v * ( - vx * vy)
+ EypBx_v * (gamma_psi - vy * vy) ) / clight
- 0.25_rt * AabssqDyp * q_mass_ratio / psi
) * q_mass_ratio / psi
- shape_dx * shape_y * (
+ EypBx_v * (gamma_psi - vy * vy) ) * clight_inv
- 0.25_rt * AabssqDyp * q_mass_ratio * psi_inv
) * q_mass_ratio * psi_inv
+ ( - shape_dx * shape_y * dx_inv * (
- vx * vy
)
- shape_x * shape_dy * (
- shape_x * shape_dy * dy_inv * (
gamma_psi - vy * vy - 1._rt
)
);
)) * clight
));

amrex::Gpu::Atomic::Add(arr.ptr(i, j, Sx),
amrex::Gpu::Atomic::Add(arr.ptr(i, j, Sx), charge_density_mu0 * (
+ shape_x * shape_y * (
+ Bz_v * vy
+ ( Ez_v * vx
+ ExmBy_v * (gamma_psi - vx * vx)
+ EypBx_v * ( - vx * vy) ) / clight
- 0.25_rt * AabssqDxp * q_mass_ratio / psi
) * q_mass_ratio / psi
+ shape_dx * shape_y * (
+ EypBx_v * ( - vx * vy) ) * clight_inv
- 0.25_rt * AabssqDxp * q_mass_ratio * psi_inv
) * q_mass_ratio * psi_inv
+ ( + shape_dx * shape_y * dx_inv * (
gamma_psi - vx * vx - 1._rt
)
+ shape_x * shape_dy * (
+ shape_x * shape_dy * dy_inv * (
- vx * vy
)
);

)) * clight
));
AlexanderSinn marked this conversation as resolved.
Show resolved Hide resolved
}
}
});
Expand Down
16 changes: 8 additions & 8 deletions src/particles/deposition/PlasmaDepositCurrent.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas
const amrex::Real dy_inv = 1._rt/dx[1];

const PhysConst pc = get_phys_const();
const amrex::Real clight = pc.c;
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;
Expand Down Expand Up @@ -168,8 +169,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas
amrex::CompileTimeOptions<false, true>, // do_tiling
#endif
amrex::CompileTimeOptions<false, true> // use_laser
>{},
{
>{}, {
Hipace::m_depos_order_xy,
Hipace::m_outer_depos_loop,
plasma.m_can_ionize,
Expand Down Expand Up @@ -199,11 +199,11 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas

const auto positions = pos_structs[ip];
if (positions.id() < 0) return;
const amrex::Real psi = psip[ip];
const amrex::Real psi_inv = 1._rt/psip[ip];
const amrex::Real xp = positions.pos(0);
const amrex::Real yp = positions.pos(1);
const amrex::Real vx_c = uxp[ip] / psi;
const amrex::Real vy_c = uyp[ip] / psi;
const amrex::Real vx_c = uxp[ip] * psi_inv;
const amrex::Real vy_c = uyp[ip] * psi_inv;

// calculate charge of the plasma particles
amrex::Real q_invvol = charge_invvol;
Expand All @@ -226,7 +226,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas

// calculate gamma/psi for plasma particles
const amrex::Real gamma_psi = 0.5_rt * (
(1._rt + 0.5_rt * Aabssqp) / (psi * psi) // TODO: fix units
(1._rt + 0.5_rt * Aabssqp) * psi_inv * psi_inv // TODO: fix units
+ vx_c * vx_c * clightinv * clightinv
+ vy_c * vy_c * clightinv * clightinv
+ 1._rt
Expand Down Expand Up @@ -269,9 +269,9 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas
// wqx, wqy wqz are particle current in each direction
const amrex::Real wqx = charge_density * vx_c;
const amrex::Real wqy = charge_density * vy_c;
const amrex::Real wqz = charge_density * (gamma_psi-1._rt) / clightinv;
const amrex::Real wqz = charge_density * (gamma_psi-1._rt) * clight;
const amrex::Real wq = charge_density * gamma_psi;
const amrex::Real wchi = charge_density * q_mu0_mass_ratio / psi;
const amrex::Real wchi = charge_density * q_mu0_mass_ratio * psi_inv;

// Deposit current into isl_arr
if (jx != -1) { // deposit_jx_jy
Expand Down
Loading