diff --git a/src/kernels/particle_pusher_gr.hpp b/src/kernels/particle_pusher_gr.hpp index 2c0d8edd..c6fffb95 100644 --- a/src/kernels/particle_pusher_gr.hpp +++ b/src/kernels/particle_pusher_gr.hpp @@ -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); } /** @@ -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 Inline void interpolateFields(index_t p, vec_t& e, vec_t& b) const; @@ -528,80 +529,185 @@ namespace kernel::gr { template requires IsCompatibleWithPusherGR + template Inline void Pusher_kernel::interpolateFields(index_t p, vec_t& e0, vec_t& 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(N_GHOSTS) }; - const int j { i2(p) + static_cast(N_GHOSTS) }; - const auto dx1_ { static_cast(dx1(p)) }; - const auto dx2_ { static_cast(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(N_GHOSTS) }; + const int j { i2(p) + static_cast(N_GHOSTS) }; + const auto dx1_ { static_cast(dx1(p)) }; + const auto dx2_ { static_cast(dx2(p)) }; + + // direct interpolation - Arno + int indx = static_cast(dx1_ + HALF); + int indy = static_cast(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(indx + ONE) - (dx1_ + HALF); + real_t pondpx = ONE - pondmx; + real_t pondmy = static_cast(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(N_GHOSTS) }; + const int j { i2(p) + static_cast(N_GHOSTS) }; + const auto dx1_ { static_cast(dx1(p)) }; + const auto dx2_ { static_cast(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(i, dx1_, ip_min, S1p); + prtl_shape::order(j, dx2_, jp_min, S2p); + // dual shape function - staggered + prtl_shape::order(i, dx1_, id_min, S1d); + prtl_shape::order(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"); + } } } @@ -707,7 +813,7 @@ namespace kernel::gr { vec_t Dp_cntrv { ZERO }, Bp_cntrv { ZERO }, Dp_hat { ZERO }, Bp_hat { ZERO }; - interpolateFields(p, Dp_cntrv, Bp_cntrv); + interpolateFields(p, Dp_cntrv, Bp_cntrv); metric.template transform(xp, Dp_cntrv, Dp_hat); metric.template transform(xp, Bp_cntrv, Bp_hat);