Skip to content

Direct field interpolation with arbitrary shape function for GR #189

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

Draft
LudwigBoess wants to merge 3 commits into
base: 1.4.0rc
Choose a base branch
from
Draft

Direct field interpolation with arbitrary shape function for GR #189

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
945b696
Direct field interpolation with arbitrary shape function for GR
LudwigBoess Mar 18, 2026
601ba8d
applied formatting
LudwigBoess Mar 18, 2026
59030d5
removed 1D and 3D versions
LudwigBoess Mar 23, 2026
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
256 changes: 181 additions & 75 deletions src/kernels/particle_pusher_gr.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,8 +159,8 @@ namespace kernel::gr {
(boundaries[0].first == PrtlBC::HORIZON);
is_absorb_i1max = (boundaries[0].second == PrtlBC::ABSORB) ||
(boundaries[0].second == PrtlBC::HORIZON);
is_axis_i2min = (boundaries[1].first == PrtlBC::AXIS);
is_axis_i2max = (boundaries[1].second == PrtlBC::AXIS);
is_axis_i2min = (boundaries[1].first == PrtlBC::AXIS);
is_axis_i2max = (boundaries[1].second == PrtlBC::AXIS);
}

/**
Expand Down Expand Up @@ -280,11 +280,12 @@ namespace kernel::gr {
// Helper functions

/**
* @brief First order Yee mesh field interpolation to particle position.
* @brief Direct field interpolation to particle position with arbitrary shape function.
* @param p index of the particle.
* @param e interpolated e-field vector of size 3 [return].
* @param b interpolated b-field vector of size 3 [return].
*/
template <unsigned short O>
Inline void interpolateFields(index_t p,
vec_t<Dim::_3D>& e,
vec_t<Dim::_3D>& b) const;
Expand Down Expand Up @@ -528,80 +529,185 @@ namespace kernel::gr {

template <class M>
requires IsCompatibleWithPusherGR<M>
template <unsigned short O>
Inline void Pusher_kernel<M>::interpolateFields(index_t p,
vec_t<Dim::_3D>& e0,
vec_t<Dim::_3D>& b0) const {
if constexpr (D == Dim::_1D) {
raise::KernelError(HERE, "interpolateFields: 1D implementation called");
} else if constexpr (D == Dim::_2D) {
// first order interpolation
// Using fields em::e = D(t = n), and em::b0 = B(t = n)

const int i { i1(p) + static_cast<int>(N_GHOSTS) };
const int j { i2(p) + static_cast<int>(N_GHOSTS) };
const auto dx1_ { static_cast<real_t>(dx1(p)) };
const auto dx2_ { static_cast<real_t>(dx2(p)) };

// first order
real_t c000, c100, c010, c110, c00, c10;

// Ex1
// interpolate to nodes
c000 = HALF * (DB(i, j, em::dx1) + DB(i - 1, j, em::dx1));
c100 = HALF * (DB(i, j, em::dx1) + DB(i + 1, j, em::dx1));
c010 = HALF * (DB(i, j + 1, em::dx1) + DB(i - 1, j + 1, em::dx1));
c110 = HALF * (DB(i, j + 1, em::dx1) + DB(i + 1, j + 1, em::dx1));
// interpolate from nodes to the particle position
c00 = c000 * (ONE - dx1_) + c100 * dx1_;
c10 = c010 * (ONE - dx1_) + c110 * dx1_;
e0[0] = c00 * (ONE - dx2_) + c10 * dx2_;
// Ex2
c000 = HALF * (DB(i, j, em::dx2) + DB(i, j - 1, em::dx2));
c100 = HALF * (DB(i + 1, j, em::dx2) + DB(i + 1, j - 1, em::dx2));
c010 = HALF * (DB(i, j, em::dx2) + DB(i, j + 1, em::dx2));
c110 = HALF * (DB(i + 1, j, em::dx2) + DB(i + 1, j + 1, em::dx2));
c00 = c000 * (ONE - dx1_) + c100 * dx1_;
c10 = c010 * (ONE - dx1_) + c110 * dx1_;
e0[1] = c00 * (ONE - dx2_) + c10 * dx2_;
// Ex3
c000 = DB(i, j, em::dx3);
c100 = DB(i + 1, j, em::dx3);
c010 = DB(i, j + 1, em::dx3);
c110 = DB(i + 1, j + 1, em::dx3);
c00 = c000 * (ONE - dx1_) + c100 * dx1_;
c10 = c010 * (ONE - dx1_) + c110 * dx1_;
e0[2] = c00 * (ONE - dx2_) + c10 * dx2_;

// Bx1
c000 = HALF * (DB0(i, j, em::bx1) + DB0(i, j - 1, em::bx1));
c100 = HALF * (DB0(i + 1, j, em::bx1) + DB0(i + 1, j - 1, em::bx1));
c010 = HALF * (DB0(i, j, em::bx1) + DB0(i, j + 1, em::bx1));
c110 = HALF * (DB0(i + 1, j, em::bx1) + DB0(i + 1, j + 1, em::bx1));
c00 = c000 * (ONE - dx1_) + c100 * dx1_;
c10 = c010 * (ONE - dx1_) + c110 * dx1_;
b0[0] = c00 * (ONE - dx2_) + c10 * dx2_;
// Bx2
c000 = HALF * (DB0(i - 1, j, em::bx2) + DB0(i, j, em::bx2));
c100 = HALF * (DB0(i, j, em::bx2) + DB0(i + 1, j, em::bx2));
c010 = HALF * (DB0(i - 1, j + 1, em::bx2) + DB0(i, j + 1, em::bx2));
c110 = HALF * (DB0(i, j + 1, em::bx2) + DB0(i + 1, j + 1, em::bx2));
c00 = c000 * (ONE - dx1_) + c100 * dx1_;
c10 = c010 * (ONE - dx1_) + c110 * dx1_;
b0[1] = c00 * (ONE - dx2_) + c10 * dx2_;
// Bx3
c000 = INV_4 * (DB0(i - 1, j - 1, em::bx3) + DB0(i - 1, j, em::bx3) +
DB0(i, j - 1, em::bx3) + DB0(i, j, em::bx3));
c100 = INV_4 * (DB0(i, j - 1, em::bx3) + DB0(i, j, em::bx3) +
DB0(i + 1, j - 1, em::bx3) + DB0(i + 1, j, em::bx3));
c010 = INV_4 * (DB0(i - 1, j, em::bx3) + DB0(i - 1, j + 1, em::bx3) +
DB0(i, j, em::bx3) + DB0(i, j + 1, em::bx3));
c110 = INV_4 * (DB0(i, j, em::bx3) + DB0(i, j + 1, em::bx3) +
DB0(i + 1, j, em::bx3) + DB0(i + 1, j + 1, em::bx3));
c00 = c000 * (ONE - dx1_) + c100 * dx1_;
c10 = c010 * (ONE - dx1_) + c110 * dx1_;
b0[2] = c00 * (ONE - dx2_) + c10 * dx2_;
} else if constexpr (D == Dim::_3D) {
raise::KernelNotImplementedError(HERE);

// Zig-zag interpolation
if constexpr (O == 0u) {

if constexpr (D == Dim::_1D) {
raise::KernelError(HERE, "1D not applicable");
} else if constexpr (D == Dim::_2D) {
const int i { i1(p) + static_cast<int>(N_GHOSTS) };
const int j { i2(p) + static_cast<int>(N_GHOSTS) };
const auto dx1_ { static_cast<real_t>(dx1(p)) };
const auto dx2_ { static_cast<real_t>(dx2(p)) };

// direct interpolation - Arno
int indx = static_cast<int>(dx1_ + HALF);
int indy = static_cast<int>(dx2_ + HALF);

// first order
real_t c000, c100, c010, c110, c00, c10;

real_t ponpmx = ONE - dx1_;
real_t ponppx = dx1_;
real_t ponpmy = ONE - dx2_;
real_t ponppy = dx2_;

real_t pondmx = static_cast<real_t>(indx + ONE) - (dx1_ + HALF);
real_t pondpx = ONE - pondmx;
real_t pondmy = static_cast<real_t>(indy + ONE) - (dx2_ + HALF);
real_t pondpy = ONE - pondmy;

// Ex1
// Interpolate --- (dual, primal)
c000 = DB(i - 1 + indx, j, em::dx1);
c100 = DB(i + indx, j, em::dx1);
c010 = DB(i - 1 + indx, j + 1, em::dx1);
c110 = DB(i + indx, j + 1, em::dx1);
c00 = c000 * pondmx + c100 * pondpx;
c10 = c010 * pondmx + c110 * pondpx;
e0[0] = c00 * ponpmy + c10 * ponppy;
// Ex2
// Interpolate -- (primal, dual)
c000 = DB(i, j - 1 + indy, em::dx2);
c100 = DB(i + 1, j - 1 + indy, em::dx2);
c010 = DB(i, j + indy, em::dx2);
c110 = DB(i + 1, j + indy, em::dx2);
c00 = c000 * ponpmx + c100 * ponppx;
c10 = c010 * ponpmx + c110 * ponppx;
e0[1] = c00 * pondmy + c10 * pondpy;
// Ex3
// Interpolate -- (primal, primal)
c000 = DB(i, j, em::dx3);
c100 = DB(i + 1, j, em::dx3);
c010 = DB(i, j + 1, em::dx3);
c110 = DB(i + 1, j + 1, em::dx3);
c00 = c000 * ponpmx + c100 * ponppx;
c10 = c010 * ponpmx + c110 * ponppx;
e0[2] = c00 * ponpmy + c10 * ponppy;

// Bx1
// Interpolate -- (primal, dual)
c000 = DB0(i, j - 1 + indy, em::bx1);
c100 = DB0(i + 1, j - 1 + indy, em::bx1);
c010 = DB0(i, j + indy, em::bx1);
c110 = DB0(i + 1, j + indy, em::bx1);
c00 = c000 * ponpmx + c100 * ponppx;
c10 = c010 * ponpmx + c110 * ponppx;
b0[0] = c00 * pondmy + c10 * pondpy;
// Bx2
// Interpolate -- (dual, primal)
c000 = DB0(i - 1 + indx, j, em::bx2);
c100 = DB0(i + indx, j, em::bx2);
c010 = DB0(i - 1 + indx, j + 1, em::bx2);
c110 = DB0(i + indx, j + 1, em::bx2);
c00 = c000 * pondmx + c100 * pondpx;
c10 = c010 * pondmx + c110 * pondpx;
b0[1] = c00 * ponpmy + c10 * ponppy;
// Bx3
// Interpolate -- (dual, dual)
c000 = DB0(i - 1 + indx, j - 1 + indy, em::bx3);
c100 = DB0(i + indx, j - 1 + indy, em::bx3);
c010 = DB0(i - 1 + indx, j + indy, em::bx3);
c110 = DB0(i + indx, j + indy, em::bx3);
c00 = c000 * pondmx + c100 * pondpx;
c10 = c010 * pondmx + c110 * pondpx;
b0[2] = c00 * pondmy + c10 * pondpy;
} else if constexpr (D == Dim::_3D) {
raise::KernelError(HERE, "3D not applicable");
}
} else if constexpr (O >= 1u) {

if constexpr (D == Dim::_1D) {
raise::KernelError(HERE, "1D not applicable");
} else if constexpr (D == Dim::_2D) {

const int i { i1(p) + static_cast<int>(N_GHOSTS) };
const int j { i2(p) + static_cast<int>(N_GHOSTS) };
const auto dx1_ { static_cast<real_t>(dx1(p)) };
const auto dx2_ { static_cast<real_t>(dx2(p)) };

// primal and dual shape function
real_t S1p[O + 1], S1d[O + 1];
real_t S2p[O + 1], S2d[O + 1];
// minimum contributing cells
int ip_min, id_min;
int jp_min, jd_min;

// primal shape function - not staggered
prtl_shape::order<false, O>(i, dx1_, ip_min, S1p);
prtl_shape::order<false, O>(j, dx2_, jp_min, S2p);
// dual shape function - staggered
prtl_shape::order<true, O>(i, dx1_, id_min, S1d);
prtl_shape::order<true, O>(j, dx2_, jd_min, S2d);

// Ex1 -- dual, primal
e0[0] = ZERO;
for (int idx2 = 0; idx2 < O + 1; idx2++) {
real_t c0 = ZERO;
for (int idx1 = 0; idx1 < O + 1; idx1++) {
c0 += S1d[idx1] * DB(id_min + idx1, jp_min + idx2, em::dx1);
}
e0[0] += c0 * S2p[idx2];
}

// Ex2 -- primal, dual
e0[1] = ZERO;
for (int idx2 = 0; idx2 < O + 1; idx2++) {
real_t c0 = ZERO;
for (int idx1 = 0; idx1 < O + 1; idx1++) {
c0 += S1p[idx1] * DB(ip_min + idx1, jd_min + idx2, em::dx2);
}
e0[1] += c0 * S2d[idx2];
}

// Ex3 -- primal, primal
e0[2] = ZERO;
for (int idx2 = 0; idx2 < O + 1; idx2++) {
real_t c0 = ZERO;
for (int idx1 = 0; idx1 < O + 1; idx1++) {
c0 += S1p[idx1] * DB(ip_min + idx1, jp_min + idx2, em::dx3);
}
e0[2] += c0 * S2p[idx2];
}

// Bx1 -- primal, dual
b0[0] = ZERO;
for (int idx2 = 0; idx2 < O + 1; idx2++) {
real_t c0 = ZERO;
for (int idx1 = 0; idx1 < O + 1; idx1++) {
c0 += S1p[idx1] * DB0(ip_min + idx1, jd_min + idx2, em::bx1);
}
b0[0] += c0 * S2d[idx2];
}

// Bx2 -- dual, primal
b0[1] = ZERO;
for (int idx2 = 0; idx2 < O + 1; idx2++) {
real_t c0 = ZERO;
for (int idx1 = 0; idx1 < O + 1; idx1++) {
c0 += S1d[idx1] * DB0(id_min + idx1, jp_min + idx2, em::bx2);
}
b0[1] += c0 * S2p[idx2];
}

// Bx3 -- dual, dual
b0[2] = ZERO;
for (int idx2 = 0; idx2 < O + 1; idx2++) {
real_t c0 = ZERO;
for (int idx1 = 0; idx1 < O + 1; idx1++) {
c0 += S1d[idx1] * DB0(id_min + idx1, jd_min + idx2, em::bx3);
}
b0[2] += c0 * S2d[idx2];
}

} else if constexpr (D == Dim::_3D) {
raise::KernelError(HERE, "3D not applicable");
}
}
}

Expand Down Expand Up @@ -707,7 +813,7 @@ namespace kernel::gr {

vec_t<Dim::_3D> Dp_cntrv { ZERO }, Bp_cntrv { ZERO }, Dp_hat { ZERO },
Bp_hat { ZERO };
interpolateFields(p, Dp_cntrv, Bp_cntrv);
interpolateFields<SHAPE_ORDER>(p, Dp_cntrv, Bp_cntrv);
metric.template transform<Idx::U, Idx::T>(xp, Dp_cntrv, Dp_hat);
metric.template transform<Idx::U, Idx::T>(xp, Bp_cntrv, Bp_hat);

Expand Down
Loading