Skip to content
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
41 changes: 7 additions & 34 deletions src/libpsc/psc_push_particles/inc_push.c
Original file line number Diff line number Diff line change
Expand Up @@ -5,39 +5,12 @@

#include "pushp.hxx"

template <typename real_t>
class PI
template <typename real_t, typename Real3 = Vec3<real_t>>
void find_idx_pos_1st_rel(Real3 final_pos, Real3 dxi, Int3& final_index,
Real3& final_pos_normalized, real_t shift)
{
public:
using Real3 = Vec3<real_t>;

PI(const Grid_t& grid) : dxi_{Real3{1., 1., 1.} / Real3(grid.domain.dx)} {}

// ----------------------------------------------------------------------
// find_idx_off_1st_rel

void find_idx_off_1st_rel(real_t xi[3], int lg[3], real_t og[3], real_t shift)
{
for (int d = 0; d < 3; d++) {
real_t pos = xi[d] * dxi_[d] + shift;
lg[d] = fint(pos);
og[d] = pos - lg[d];
}
for (int d = 0; d < 3; d++) {
final_pos_normalized[d] = final_pos[d] * dxi[d] + shift;
final_index[d] = fint(final_pos_normalized[d]);
}

// ----------------------------------------------------------------------
// find_idx_off_pos_1st_rel

void find_idx_off_pos_1st_rel(real_t xi[3], int lg[3], real_t og[3],
real_t pos[3], real_t shift)
{
for (int d = 0; d < 3; d++) {
pos[d] = xi[d] * dxi_[d] + shift;
lg[d] = fint(pos[d]);
og[d] = pos[d] - lg[d];
}
}

private:
Real3 dxi_;
};
}
28 changes: 12 additions & 16 deletions src/libpsc/psc_push_particles/push_particles_1vb.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ struct PushParticlesVb
static void push_mprts(Mparticles& mprts, MfieldsState& mflds)
{
const auto& grid = mprts.grid();
PI<real_t> pi(grid);
Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx);
real_t dq_kind[MAX_NR_KINDS];
auto& kinds = grid.kinds;
Expand All @@ -46,13 +45,8 @@ struct PushParticlesVb
flds.storage().view(_all, _all, _all, _s(JXI, JXI + 3)) = real_t(0);

for (auto prt : prts) {
Real3& x = prt.x();

real_t xm[3];
for (int d = 0; d < 3; d++) {
xm[d] = x[d] * dxi[d];
}
ip.set_coeffs(xm);
Real3 initial_pos_normalized = prt.x() * dxi;
ip.set_coeffs(initial_pos_normalized);

// FIELD INTERPOLATION
Real3 E = {ip.ex(EM), ip.ey(EM), ip.ez(EM)};
Expand All @@ -64,11 +58,12 @@ struct PushParticlesVb

// x^(n+0.5), p^(n+1.0) -> x^(n+1.5), p^(n+1.0)
auto v = advance.calc_v(prt.u());
advance.push_x(x, v);
advance.push_x(prt.x(), v);

int lf[3];
real_t of[3], xp[3];
pi.find_idx_off_pos_1st_rel(x, lf, of, xp, real_t(0.));
Int3 final_index;
Real3 final_pos_normalized;
find_idx_pos_1st_rel(prt.x(), dxi, final_index, final_pos_normalized,
real_t(0.));

// CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt
int lg[3];
Expand All @@ -81,7 +76,8 @@ struct PushParticlesVb
if (!Dim::InvarZ::value) {
lg[2] = ip.cz.g.l;
}
current.calc_j(J, xm, xp, lf, lg, prt.qni_wni(), v);
current.calc_j(J, initial_pos_normalized, final_pos_normalized,
final_index, lg, prt.qni_wni(), v);
}
}
}
Expand Down Expand Up @@ -112,14 +108,14 @@ struct PushParticlesVb
// field interpolation
real_t* x = prt.x;

real_t xm[3];
real_t initial_pos_normalized[3];
for (int d = 0; d < 3; d++) {
xm[d] = x[d] * dxi[d];
initial_pos_normalized[d] = x[d] * dxi[d];
}

// FIELD INTERPOLATION

ip.set_coeffs(xm);
ip.set_coeffs(initial_pos_normalized);
// FIXME, we're not using EM instead flds_em
real_t E[3] = {ip.ex(EM), ip.ey(EM), ip.ez(EM)};
real_t H[3] = {ip.hx(EM), ip.hy(EM), ip.hz(EM)};
Expand Down