From 6c0eb5e193478b98156bf9e7d4189ed092c8e328 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 10:58:22 -0500 Subject: [PATCH 01/16] inc_curr_1vb_split: name ip --- src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index 393ed033bf..05683cdf55 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -63,7 +63,8 @@ struct Current1vbSplit { const int dim = 0; int im = fint(xm[dim]); - if (!Dim::InvarX::value && (fint(xp[dim]) != im)) { + int ip = fint(xp[dim]); + if (!Dim::InvarX::value && ip != im) { real_t x1[3]; calc_j2_split_along_dim(dim, im, x1, xm, xp); calc_j2_one_cell(curr_cache, qni_wni, xm, x1); @@ -78,7 +79,8 @@ struct Current1vbSplit { const int dim = 1; int im = fint(xm[dim]); - if (!Dim::InvarY::value && (fint(xp[dim]) != im)) { + int ip = fint(xp[dim]); + if (!Dim::InvarY::value && ip != im) { real_t x1[3]; calc_j2_split_along_dim(dim, im, x1, xm, xp); calc_j2_split_dim_x(curr_cache, qni_wni, xm, x1); @@ -93,7 +95,8 @@ struct Current1vbSplit { const int dim = 2; int im = fint(xm[dim]); - if (!Dim::InvarZ::value && (fint(xp[dim]) != im)) { + int ip = fint(xp[dim]); + if (!Dim::InvarZ::value && ip != im) { real_t x1[3]; calc_j2_split_along_dim(dim, im, x1, xm, xp); calc_j2_split_dim_y(curr_cache, qni_wni, xm, x1); From dd7da4074788927b7044249e9a8d55346e3b039b Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:02:27 -0500 Subject: [PATCH 02/16] inc_curr_1vb_split: use Dim::is_invar --- src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index 05683cdf55..fc27523362 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -61,10 +61,10 @@ struct Current1vbSplit void calc_j2_split_dim_x(fields_t curr_cache, real_t qni_wni, const real_t* xm, const real_t* xp) { - const int dim = 0; + constexpr int dim = 0; int im = fint(xm[dim]); int ip = fint(xp[dim]); - if (!Dim::InvarX::value && ip != im) { + if (!Dim::is_invar(dim) && ip != im) { real_t x1[3]; calc_j2_split_along_dim(dim, im, x1, xm, xp); calc_j2_one_cell(curr_cache, qni_wni, xm, x1); @@ -77,10 +77,10 @@ struct Current1vbSplit void calc_j2_split_dim_y(fields_t curr_cache, real_t qni_wni, const real_t* xm, const real_t* xp) { - const int dim = 1; + constexpr int dim = 1; int im = fint(xm[dim]); int ip = fint(xp[dim]); - if (!Dim::InvarY::value && ip != im) { + if (!Dim::is_invar(dim) && ip != im) { real_t x1[3]; calc_j2_split_along_dim(dim, im, x1, xm, xp); calc_j2_split_dim_x(curr_cache, qni_wni, xm, x1); @@ -93,10 +93,10 @@ struct Current1vbSplit void calc_j2_split_dim_z(fields_t curr_cache, real_t qni_wni, const real_t* xm, const real_t* xp) { - const int dim = 2; + constexpr int dim = 2; int im = fint(xm[dim]); int ip = fint(xp[dim]); - if (!Dim::InvarZ::value && ip != im) { + if (!Dim::is_invar(dim) && ip != im) { real_t x1[3]; calc_j2_split_along_dim(dim, im, x1, xm, xp); calc_j2_split_dim_y(curr_cache, qni_wni, xm, x1); From 6f72fe92c785f3d43e1647f712e0b49097fb7ba1 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:06:53 -0500 Subject: [PATCH 03/16] inc_curr_1vb_split: todo if constexpr --- src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index fc27523362..9b6b065ac2 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -58,6 +58,8 @@ struct Current1vbSplit } } + // TODO c++17 combine these calc_j2_split_dim_* with if constexpr + void calc_j2_split_dim_x(fields_t curr_cache, real_t qni_wni, const real_t* xm, const real_t* xp) { From 2391f5c8023e3716675bc74b90e032effa25ab6a Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:08:10 -0500 Subject: [PATCH 04/16] inc_curr_1vb_split: pass ip --- src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index 9b6b065ac2..455da9ddb5 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -35,11 +35,10 @@ struct Current1vbSplit deposition_(curr_cache, i, qni_wni, dx, xa, Dim{}); } - static void calc_j2_split_along_dim(int dim, int im, real_t x1[3], + static void calc_j2_split_along_dim(int dim, int im, int ip, real_t x1[3], const real_t xm[3], const real_t xp[3]) { - real_t bnd = 0.f; // quell warning - int ip = fint(xp[dim]); + real_t bnd = 0.f; // quell warning if (ip == im + 1) { // crossed boundary to right bnd = ip; } else if (ip == im - 1) { // crosses boundary to left @@ -68,7 +67,7 @@ struct Current1vbSplit int ip = fint(xp[dim]); if (!Dim::is_invar(dim) && ip != im) { real_t x1[3]; - calc_j2_split_along_dim(dim, im, x1, xm, xp); + calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_one_cell(curr_cache, qni_wni, xm, x1); calc_j2_one_cell(curr_cache, qni_wni, x1, xp); } else { @@ -84,7 +83,7 @@ struct Current1vbSplit int ip = fint(xp[dim]); if (!Dim::is_invar(dim) && ip != im) { real_t x1[3]; - calc_j2_split_along_dim(dim, im, x1, xm, xp); + calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_split_dim_x(curr_cache, qni_wni, xm, x1); calc_j2_split_dim_x(curr_cache, qni_wni, x1, xp); } else { @@ -100,7 +99,7 @@ struct Current1vbSplit int ip = fint(xp[dim]); if (!Dim::is_invar(dim) && ip != im) { real_t x1[3]; - calc_j2_split_along_dim(dim, im, x1, xm, xp); + calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_split_dim_y(curr_cache, qni_wni, xm, x1); calc_j2_split_dim_y(curr_cache, qni_wni, x1, xp); } else { From 2cbb3aa4b99317192d4ae79c0eb0b743cbad62d2 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:10:03 -0500 Subject: [PATCH 05/16] inc_curr_1vb_split: swap if/else for readability --- .../psc_push_particles/inc_curr_1vb_split.cxx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index 455da9ddb5..51b4ef5301 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -65,13 +65,13 @@ struct Current1vbSplit constexpr int dim = 0; int im = fint(xm[dim]); int ip = fint(xp[dim]); - if (!Dim::is_invar(dim) && ip != im) { + if (Dim::is_invar(dim) || ip == im) { + calc_j2_one_cell(curr_cache, qni_wni, xm, xp); + } else { real_t x1[3]; calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_one_cell(curr_cache, qni_wni, xm, x1); calc_j2_one_cell(curr_cache, qni_wni, x1, xp); - } else { - calc_j2_one_cell(curr_cache, qni_wni, xm, xp); } } @@ -81,13 +81,13 @@ struct Current1vbSplit constexpr int dim = 1; int im = fint(xm[dim]); int ip = fint(xp[dim]); - if (!Dim::is_invar(dim) && ip != im) { + if (Dim::is_invar(dim) || ip == im) { + calc_j2_split_dim_x(curr_cache, qni_wni, xm, xp); + } else { real_t x1[3]; calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_split_dim_x(curr_cache, qni_wni, xm, x1); calc_j2_split_dim_x(curr_cache, qni_wni, x1, xp); - } else { - calc_j2_split_dim_x(curr_cache, qni_wni, xm, xp); } } @@ -97,13 +97,13 @@ struct Current1vbSplit constexpr int dim = 2; int im = fint(xm[dim]); int ip = fint(xp[dim]); - if (!Dim::is_invar(dim) && ip != im) { + if (Dim::is_invar(dim) || ip == im) { + calc_j2_split_dim_y(curr_cache, qni_wni, xm, xp); + } else { real_t x1[3]; calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_split_dim_y(curr_cache, qni_wni, xm, x1); calc_j2_split_dim_y(curr_cache, qni_wni, x1, xp); - } else { - calc_j2_split_dim_y(curr_cache, qni_wni, xm, xp); } } From c2d1c8eddb820e20675083b3b2e8bc3baf3b7607 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:13:57 -0500 Subject: [PATCH 06/16] inc_curr_1vb_split: bnd = std::max --- src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index 51b4ef5301..5941ac84ef 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -38,14 +38,9 @@ struct Current1vbSplit static void calc_j2_split_along_dim(int dim, int im, int ip, real_t x1[3], const real_t xm[3], const real_t xp[3]) { - real_t bnd = 0.f; // quell warning - if (ip == im + 1) { // crossed boundary to right - bnd = ip; - } else if (ip == im - 1) { // crosses boundary to left - bnd = im; - } else { - assert(0); - } + // boundary is at the lower edge of the cell with the higher index + real_t bnd = std::max(im, ip); + real_t frac = (bnd - xm[dim]) / (xp[dim] - xm[dim]); // FIXME, set d == dim value to exact boundary? for (int d = 0; d < 3; d++) { From 9a5ac4288f8798332f2740a79b9d487bc5dcc339 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:26:01 -0500 Subject: [PATCH 07/16] inc_curr_1vb_split: rm fixed fixme --- src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index 5941ac84ef..08c060448b 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -42,7 +42,6 @@ struct Current1vbSplit real_t bnd = std::max(im, ip); real_t frac = (bnd - xm[dim]) / (xp[dim] - xm[dim]); - // FIXME, set d == dim value to exact boundary? for (int d = 0; d < 3; d++) { if (d == dim) { x1[d] = bnd; From 71836c7d7ece1dfc8d1521a924e09bc9a7b15485 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:35:57 -0500 Subject: [PATCH 08/16] inc_curr_1vb_split: use Real3, Int3 --- .../psc_push_particles/inc_curr_1vb_split.cxx | 37 +++++++++---------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index 08c060448b..0c09d29f94 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -19,14 +19,13 @@ struct Current1vbSplit deposition_(real_t(grid.norm.fnqs / grid.dt) * Real3(grid.domain.dx)) {} - void calc_j2_one_cell(fields_t curr_cache, real_t qni_wni, const real_t xm[3], - const real_t xp[3]) + void calc_j2_one_cell(fields_t curr_cache, real_t qni_wni, const Real3& xm, + const Real3& xp) { - real_t dx[3] = {xp[0] - xm[0], xp[1] - xm[1], xp[2] - xm[2]}; - real_t xa[3] = {.5f * (xm[0] + xp[0]), .5f * (xm[1] + xp[1]), - .5f * (xm[2] + xp[2])}; + Real3 dx = xp - xm; + Real3 xa = real_t(0.5) * (xp + xm); - int i[3]; + Int3 i; for (int d = 0; d < 3; d++) { i[d] = fint(xa[d]); xa[d] -= i[d]; @@ -35,8 +34,8 @@ struct Current1vbSplit deposition_(curr_cache, i, qni_wni, dx, xa, Dim{}); } - static void calc_j2_split_along_dim(int dim, int im, int ip, real_t x1[3], - const real_t xm[3], const real_t xp[3]) + static void calc_j2_split_along_dim(int dim, int im, int ip, Real3& x1, + const Real3& xm, const Real3& xp) { // boundary is at the lower edge of the cell with the higher index real_t bnd = std::max(im, ip); @@ -53,8 +52,8 @@ struct Current1vbSplit // TODO c++17 combine these calc_j2_split_dim_* with if constexpr - void calc_j2_split_dim_x(fields_t curr_cache, real_t qni_wni, - const real_t* xm, const real_t* xp) + void calc_j2_split_dim_x(fields_t curr_cache, real_t qni_wni, const Real3& xm, + const Real3& xp) { constexpr int dim = 0; int im = fint(xm[dim]); @@ -62,15 +61,15 @@ struct Current1vbSplit if (Dim::is_invar(dim) || ip == im) { calc_j2_one_cell(curr_cache, qni_wni, xm, xp); } else { - real_t x1[3]; + Real3 x1; calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_one_cell(curr_cache, qni_wni, xm, x1); calc_j2_one_cell(curr_cache, qni_wni, x1, xp); } } - void calc_j2_split_dim_y(fields_t curr_cache, real_t qni_wni, - const real_t* xm, const real_t* xp) + void calc_j2_split_dim_y(fields_t curr_cache, real_t qni_wni, const Real3& xm, + const Real3& xp) { constexpr int dim = 1; int im = fint(xm[dim]); @@ -78,15 +77,15 @@ struct Current1vbSplit if (Dim::is_invar(dim) || ip == im) { calc_j2_split_dim_x(curr_cache, qni_wni, xm, xp); } else { - real_t x1[3]; + Real3 x1; calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_split_dim_x(curr_cache, qni_wni, xm, x1); calc_j2_split_dim_x(curr_cache, qni_wni, x1, xp); } } - void calc_j2_split_dim_z(fields_t curr_cache, real_t qni_wni, - const real_t* xm, const real_t* xp) + void calc_j2_split_dim_z(fields_t curr_cache, real_t qni_wni, const Real3& xm, + const Real3& xp) { constexpr int dim = 2; int im = fint(xm[dim]); @@ -94,7 +93,7 @@ struct Current1vbSplit if (Dim::is_invar(dim) || ip == im) { calc_j2_split_dim_y(curr_cache, qni_wni, xm, xp); } else { - real_t x1[3]; + Real3 x1; calc_j2_split_along_dim(dim, im, ip, x1, xm, xp); calc_j2_split_dim_y(curr_cache, qni_wni, xm, x1); calc_j2_split_dim_y(curr_cache, qni_wni, x1, xp); @@ -104,8 +103,8 @@ struct Current1vbSplit // ---------------------------------------------------------------------- // calc_j - void calc_j(fields_t curr_cache, real_t* xm, real_t* xp, int* lf, int* lg, - real_t qni_wni, real_t* vxi) + void calc_j(fields_t curr_cache, Real3& xm, Real3& xp, const Int3& lf, + const Int3& lg, real_t qni_wni, Real3& vxi) { // this way, we guarantee that the average position in invar directions // will remain in the 0th cell From 505f98db70affbceb65b487dd5c45cbc76f0a2ed Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 13:36:52 -0500 Subject: [PATCH 09/16] current_deposition: use Vec3 --- src/include/psc/current_deposition.hxx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/include/psc/current_deposition.hxx b/src/include/psc/current_deposition.hxx index f46d5a161c..7781e953eb 100644 --- a/src/include/psc/current_deposition.hxx +++ b/src/include/psc/current_deposition.hxx @@ -14,11 +14,11 @@ public: GT_INLINE CurrentDeposition1vb(real3_t fnqs) : fnqs_{fnqs} {} - GT_INLINE void operator()(Curr& curr, const int i[3], real_t qni_wni, - const real_t dx[3], const real_t xa[3], + GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni, + const real3_t& dx, const real3_t& xa, dim_xyz tag_dim) const { - real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2]; + real_t h = (1.f / real_t(12.f)) * dx.prod(); real3_t fnq = qni_wni * fnqs_; curr.add(0, i[0], i[1], i[2], @@ -49,11 +49,11 @@ public: fnq[2] * (dx[2] * (xa[0]) * (xa[1]) + h)); }; - GT_INLINE void operator()(Curr& curr, const int i[3], real_t qni_wni, - const real_t dx[3], const real_t xa[3], + GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni, + const real3_t& dx, const real3_t& xa, dim_xz tag_dim) const { - real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2]; + real_t h = (1.f / real_t(12.f)) * dx.prod(); real3_t fnq = qni_wni * fnqs_; curr.add(0, i[0], i[1], i[2], fnq[0] * (dx[0] * (1.f - xa[2]))); @@ -72,11 +72,11 @@ public: curr.add(2, i[0] + 1, i[1], i[2], fnq[2] * (dx[2] * (xa[0]))); } - GT_INLINE void operator()(Curr& curr, const int i[3], real_t qni_wni, - const real_t dx[3], const real_t xa[3], + GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni, + const real3_t& dx, const real3_t& xa, dim_yz tag_dim) const { - real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2]; + real_t h = (1.f / real_t(12.f)) * dx.prod(); real3_t fnq = qni_wni * fnqs_; curr.add(0, i[0], i[1], i[2], From 30ad921a48baf6828646d033de807b5beb685cf7 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 14:18:27 -0500 Subject: [PATCH 10/16] *: use Vec3 everywhere --- src/include/psc/current_deposition.hxx | 18 ++-- src/include/psc/deposit.hxx | 66 ++++++-------- src/include/psc/moment.hxx | 14 +-- src/libpsc/psc_push_fields/marder_impl.hxx | 8 +- src/libpsc/tests/test_current_deposition.cxx | 94 ++++++++++---------- src/libpsc/tests/test_deposit.cxx | 28 +++--- 6 files changed, 109 insertions(+), 119 deletions(-) diff --git a/src/include/psc/current_deposition.hxx b/src/include/psc/current_deposition.hxx index 7781e953eb..2ec3e391aa 100644 --- a/src/include/psc/current_deposition.hxx +++ b/src/include/psc/current_deposition.hxx @@ -10,16 +10,16 @@ class CurrentDeposition1vb public: using Curr = C; using real_t = typename Curr::real_t; - using real3_t = Vec3; + using Real3 = Vec3; - GT_INLINE CurrentDeposition1vb(real3_t fnqs) : fnqs_{fnqs} {} + GT_INLINE CurrentDeposition1vb(Real3 fnqs) : fnqs_{fnqs} {} GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni, - const real3_t& dx, const real3_t& xa, + const Real3& dx, const Real3& xa, dim_xyz tag_dim) const { real_t h = (1.f / real_t(12.f)) * dx.prod(); - real3_t fnq = qni_wni * fnqs_; + Real3 fnq = qni_wni * fnqs_; curr.add(0, i[0], i[1], i[2], fnq[0] * (dx[0] * (1.f - xa[1]) * (1.f - xa[2]) + h)); @@ -50,11 +50,11 @@ public: }; GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni, - const real3_t& dx, const real3_t& xa, + const Real3& dx, const Real3& xa, dim_xz tag_dim) const { real_t h = (1.f / real_t(12.f)) * dx.prod(); - real3_t fnq = qni_wni * fnqs_; + Real3 fnq = qni_wni * fnqs_; curr.add(0, i[0], i[1], i[2], fnq[0] * (dx[0] * (1.f - xa[2]))); curr.add(0, i[0], i[1], i[2] + 1, fnq[0] * (dx[0] * (xa[2]))); @@ -73,12 +73,12 @@ public: } GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni, - const real3_t& dx, const real3_t& xa, + const Real3& dx, const Real3& xa, dim_yz tag_dim) const { real_t h = (1.f / real_t(12.f)) * dx.prod(); - real3_t fnq = qni_wni * fnqs_; + Real3 fnq = qni_wni * fnqs_; curr.add(0, i[0], i[1], i[2], fnq[0] * (dx[0] * (1.f - xa[1]) * (1.f - xa[2]) + h)); curr.add(0, i[0], i[1] + 1, i[2], @@ -96,7 +96,7 @@ public: } private: - real3_t fnqs_; + Real3 fnqs_; }; }; // namespace psc diff --git a/src/include/psc/deposit.hxx b/src/include/psc/deposit.hxx index af0a89f3f8..f00d3a9235 100644 --- a/src/include/psc/deposit.hxx +++ b/src/include/psc/deposit.hxx @@ -26,13 +26,13 @@ class Deposit1st { public: using real_t = R; - using real3_t = gt::sarray; + using Real3 = Vec3; using dim_t = D; private: template - void operator()(F& flds, const gt::sarray& l, const real3_t& h, - real_t val, dim_yz tag) + void operator()(F& flds, const Int3& l, const Real3& h, real_t val, + dim_yz tag) { flds(0, l[1] + 0, l[2] + 0) += val * (1.f - h[1]) * (1.f - h[2]); flds(0, l[1] + 1, l[2] + 0) += val * h[1] * (1.f - h[2]); @@ -41,8 +41,8 @@ private: } template - void operator()(F& flds, const gt::sarray& l, const real3_t& h, - real_t val, dim_xyz tag) + void operator()(F& flds, const Int3& l, const Real3& h, real_t val, + dim_xyz tag) { // clang-format off flds(l[0] + 0, l[1] + 0, l[2] + 0) += val * (1.f - h[0]) * (1.f - h[1]) * (1.f - h[2]); @@ -58,8 +58,7 @@ private: public: template - void operator()(F& flds, const gt::sarray& l, const real3_t& h, - real_t val) + void operator()(F& flds, const Int3& l, const Real3& h, real_t val) { (*this)(flds, l, h, val, dim_t{}); } @@ -76,13 +75,13 @@ class Deposit2nd { public: using real_t = R; - using real3_t = gt::sarray; + using Real3 = Vec3; using dim_t = D; private: template - void operator()(F& flds, const gt::sarray& l, const real3_t& h, - real_t val, dim_yz tag) + void operator()(F& flds, const Int3& l, const Real3& h, real_t val, + dim_yz tag) { real_t gmy = .5f * (.5f + h[1]) * (.5f + h[1]); real_t gmz = .5f * (.5f + h[2]) * (.5f + h[2]); @@ -105,8 +104,8 @@ private: } template - void operator()(F& flds, const gt::sarray& l, const real3_t& h, - real_t val, dim_xyz tag) + void operator()(F& flds, const Int3& l, const Real3& h, real_t val, + dim_xyz tag) { real_t gmx = .5f * (.5f + h[0]) * (.5f + h[0]); real_t gmy = .5f * (.5f + h[1]) * (.5f + h[1]); @@ -151,8 +150,7 @@ private: public: template - void operator()(F& flds, const gt::sarray& l, const real3_t& h, - real_t val) + void operator()(F& flds, const Int3& l, const Real3& h, real_t val) { (*this)(flds, l, h, val, dim_t{}); } @@ -178,12 +176,10 @@ public: static std::string suffix() { return "_1st_nc"; } template - void operator()(F& flds, const gt::sarray& ib, - const gt::sarray& x, T val) - + void operator()(F& flds, const Int3& ib, const Vec3& x, T val) { - gt::sarray l; - gt::sarray h; + Int3 l; + Vec3 h; for (int d = 0; d < 3; d++) { l[d] = fint(x[d]); h[d] = x[d] - l[d]; @@ -201,12 +197,10 @@ public: static std::string suffix() { return "_1st_cc"; } template - void operator()(F& flds, const gt::sarray& ib, - const gt::sarray& x, T val) - + void operator()(F& flds, const Int3& ib, const Vec3& x, T val) { - gt::sarray l; - gt::sarray h; + Int3 l; + Vec3 h; for (int d = 0; d < 3; d++) { l[d] = fint(x[d] - .5f); h[d] = x[d] - .5f - l[d]; @@ -228,12 +222,10 @@ public: static std::string suffix() { return "_2nd_nc"; } template - void operator()(F& flds, const gt::sarray& ib, - const gt::sarray& x, T val) - + void operator()(F& flds, const Int3& ib, const Vec3& x, T val) { - gt::sarray l; - gt::sarray h; + Int3 l; + Vec3 h; for (int d = 0; d < 3; d++) { l[d] = nint(x[d]); h[d] = l[d] - x[d]; // negated! @@ -245,14 +237,14 @@ public: }; template -void nc(F& flds, const gt::sarray& ib, const gt::sarray& x, T val) +void nc(F& flds, const Int3& ib, const Vec3& x, T val) { Deposit1stNc deposit; deposit(flds, ib, x, val); } template -void cc(F& flds, const gt::sarray& ib, const gt::sarray& x, T val) +void cc(F& flds, const Int3& ib, const Vec3& x, T val) { Deposit1stCc deposit; deposit(flds, ib, x, val); @@ -276,27 +268,25 @@ class Deposit public: using real_t = R; using dim_t = D; - using real3_t = gt::sarray; + using Real3 = Vec3; using DepositNorm = DEPOSIT_NORM; static const centering::Centering CENTERING = DepositNorm::CENTERING; static std::string suffix() { return DepositNorm::suffix(); } - Deposit(const real3_t& dx, real_t fnqs) : dxi_{real_t(1.) / dx}, fnqs_{fnqs} - {} + Deposit(const Real3& dx, real_t fnqs) : dxi_{real_t(1.) / dx}, fnqs_{fnqs} {} template - void operator()(const F& flds, const gt::shape_type<3>& ib, const real3_t& xi, - real_t val) + void operator()(const F& flds, const Int3& ib, const Real3& xi, real_t val) { - real3_t x = xi * dxi_; + Real3 x = xi * dxi_; real_t value = fnqs_ * val; DepositNorm deposit; deposit(flds, ib, x, value); } - real3_t dxi_; + Real3 dxi_; real_t fnqs_; }; diff --git a/src/include/psc/moment.hxx b/src/include/psc/moment.hxx index 39672f7ed5..08a108714d 100644 --- a/src/include/psc/moment.hxx +++ b/src/include/psc/moment.hxx @@ -18,11 +18,11 @@ struct DepositParticlesContext using storage_type = S; using dim_t = D; using real_t = typename storage_type::value_type; - using real3_t = gt::sarray; + using Real3 = Vec3; using Deposit = DepositCode; DepositParticlesContext(storage_type& mflds_gt, const Int3& ib, - const real3_t& dx, real_t fnqs) + const Real3& dx, real_t fnqs) : mflds_gt{mflds_gt}, ib{ib}, deposit{dx, fnqs} {} @@ -35,7 +35,7 @@ struct DepositParticlesContext Int3 ib; Deposit deposit; int p; - real3_t x; + Real3 x; }; template ; + using Real3 = Vec3; using DepositCtx = DepositParticlesContext; template - void operator()(const real3_t& dx, real_t fnqs, storage_type& mflds_gt, + void operator()(const Real3& dx, real_t fnqs, storage_type& mflds_gt, const Int3& ib, const Mparticles& mprts, F&& func) { auto accessor = mprts.accessor(); @@ -70,10 +70,10 @@ template