Skip to content

Partially cleanup inc_curr_1vb_split #375

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

Open
JamesMcClung wants to merge 16 commits into
base: main
Choose a base branch
from
Open

Partially cleanup inc_curr_1vb_split #375

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
152209a
inc_curr_1vb_split: name ip
JamesMcClung Dec 31, 2025
b467fd0
inc_curr_1vb_split: use Dim::is_invar
JamesMcClung Dec 31, 2025
2d83ae5
inc_curr_1vb_split: todo if constexpr
JamesMcClung Dec 31, 2025
d05f91f
inc_curr_1vb_split: pass ip
JamesMcClung Dec 31, 2025
ee4e4c8
inc_curr_1vb_split: swap if/else for readability
JamesMcClung Dec 31, 2025
2cc2858
inc_curr_1vb_split: bnd = std::max
JamesMcClung Dec 31, 2025
e472797
inc_curr_1vb_split: rm fixed fixme
JamesMcClung Dec 31, 2025
24ac72a
inc_curr_1vb_split: use Real3, Int3
JamesMcClung Dec 31, 2025
6edcc3b
current_deposition: use Vec3
JamesMcClung Dec 31, 2025
fe30582
*: use Vec3 everywhere
JamesMcClung Dec 31, 2025
9d6caac
gtensor: rm operators
JamesMcClung Dec 31, 2025
3184c30
inc_curr_1vb_split: return x1
JamesMcClung Dec 31, 2025
e415c6a
inc_curr_1vb_split: use inv
JamesMcClung Dec 31, 2025
cb0506f
inc_curr_1vb_split: use fint
JamesMcClung Dec 31, 2025
51e673e
current_deposition: format to be obviously correct
JamesMcClung Dec 31, 2025
613303b
*: citation comments
JamesMcClung Dec 31, 2025
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
130 changes: 58 additions & 72 deletions src/include/psc/current_deposition.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,93 +10,79 @@ class CurrentDeposition1vb
public:
using Curr = C;
using real_t = typename Curr::real_t;
using real3_t = Vec3<real_t>;
using Real3 = Vec3<real_t>;

GT_INLINE CurrentDeposition1vb(real3_t fnqs) : fnqs_{fnqs} {}
GT_INLINE CurrentDeposition1vb(Real3 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& dx, const Real3& xa,
dim_xyz tag_dim) const
{
real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2];
real3_t 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],
fnq[0] * (dx[0] * (xa[1]) * (1.f - xa[2]) - h));
curr.add(0, i[0], i[1], i[2] + 1,
fnq[0] * (dx[0] * (1.f - xa[1]) * (xa[2]) - h));
curr.add(0, i[0], i[1] + 1, i[2] + 1,
fnq[0] * (dx[0] * (xa[1]) * (xa[2]) + h));

curr.add(1, i[0], i[1], i[2],
fnq[1] * (dx[1] * (1.f - xa[0]) * (1.f - xa[2]) + h));
curr.add(1, i[0] + 1, i[1], i[2],
fnq[1] * (dx[1] * (xa[0]) * (1.f - xa[2]) - h));
curr.add(1, i[0], i[1], i[2] + 1,
fnq[1] * (dx[1] * (1.f - xa[0]) * (xa[2]) - h));
curr.add(1, i[0] + 1, i[1], i[2] + 1,
fnq[1] * (dx[1] * (xa[0]) * (xa[2]) + h));

curr.add(2, i[0], i[1], i[2],
fnq[2] * (dx[2] * (1.f - xa[0]) * (1.f - xa[1]) + h));
curr.add(2, i[0] + 1, i[1], i[2],
fnq[2] * (dx[2] * (xa[0]) * (1.f - xa[1]) - h));
curr.add(2, i[0], i[1] + 1, i[2],
fnq[2] * (dx[2] * (1.f - xa[0]) * (xa[1]) - h));
curr.add(2, i[0] + 1, i[1] + 1, i[2],
fnq[2] * (dx[2] * (xa[0]) * (xa[1]) + h));
real_t h = (1.f / real_t(12.f)) * dx.prod();
Real3 fnq = qni_wni * fnqs_;

// clang-format off
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], fnq[0] * (dx[0] * ( xa[1]) * (1.f - xa[2]) - h));
curr.add(0, i[0], i[1], i[2] + 1, fnq[0] * (dx[0] * (1.f - xa[1]) * ( xa[2]) - h));
curr.add(0, i[0], i[1] + 1, i[2] + 1, fnq[0] * (dx[0] * ( xa[1]) * ( xa[2]) + h));

curr.add(1, i[0], i[1], i[2], fnq[1] * (dx[1] * (1.f - xa[2]) * (1.f - xa[0]) + h));
curr.add(1, i[0], i[1], i[2] + 1, fnq[1] * (dx[1] * ( xa[2]) * (1.f - xa[0]) - h));
curr.add(1, i[0] + 1, i[1], i[2], fnq[1] * (dx[1] * (1.f - xa[2]) * ( xa[0]) - h));
curr.add(1, i[0] + 1, i[1], i[2] + 1, fnq[1] * (dx[1] * ( xa[2]) * ( xa[0]) + h));

curr.add(2, i[0], i[1], i[2], fnq[2] * (dx[2] * (1.f - xa[0]) * (1.f - xa[1]) + h));
curr.add(2, i[0] + 1, i[1], i[2], fnq[2] * (dx[2] * ( xa[0]) * (1.f - xa[1]) - h));
curr.add(2, i[0], i[1] + 1, i[2], fnq[2] * (dx[2] * (1.f - xa[0]) * ( xa[1]) - h));
curr.add(2, i[0] + 1, i[1] + 1, i[2], fnq[2] * (dx[2] * ( xa[0]) * ( xa[1]) + h));
// clang-format on
};

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& dx, const Real3& xa,
dim_xz tag_dim) const
{
real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2];
real3_t 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])));

curr.add(1, i[0], i[1], i[2],
fnq[1] * (dx[1] * (1.f - xa[0]) * (1.f - xa[2]) + h));
curr.add(1, i[0] + 1, i[1], i[2],
fnq[1] * (dx[1] * (xa[0]) * (1.f - xa[2]) - h));
curr.add(1, i[0], i[1], i[2] + 1,
fnq[1] * (dx[1] * (1.f - xa[0]) * (xa[2]) - h));
curr.add(1, i[0] + 1, i[1], i[2] + 1,
fnq[1] * (dx[1] * (xa[0]) * (xa[2]) + h));

curr.add(2, i[0], i[1], i[2], fnq[2] * (dx[2] * (1.f - xa[0])));
curr.add(2, i[0] + 1, i[1], i[2], fnq[2] * (dx[2] * (xa[0])));
real_t h = (1.f / real_t(12.f)) * dx.prod();
Real3 fnq = qni_wni * fnqs_;

// clang-format off
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]) ));

curr.add(1, i[0], i[1], i[2], fnq[1] * (dx[1] * (1.f - xa[2]) * (1.f - xa[0]) + h));
curr.add(1, i[0], i[1], i[2] + 1, fnq[1] * (dx[1] * ( xa[2]) * (1.f - xa[0]) - h));
curr.add(1, i[0] + 1, i[1], i[2], fnq[1] * (dx[1] * (1.f - xa[2]) * ( xa[0]) - h));
curr.add(1, i[0] + 1, i[1], i[2] + 1, fnq[1] * (dx[1] * ( xa[2]) * ( xa[0]) + h));

curr.add(2, i[0], i[1], i[2], fnq[2] * (dx[2] * (1.f - xa[0]) ));
curr.add(2, i[0] + 1, i[1], i[2], fnq[2] * (dx[2] * ( xa[0]) ));
// clang-format on
}

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& dx, const Real3& xa,
dim_yz tag_dim) const
{
real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2];

real3_t 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],
fnq[0] * (dx[0] * (xa[1]) * (1.f - xa[2]) - h));
curr.add(0, i[0], i[1], i[2] + 1,
fnq[0] * (dx[0] * (1.f - xa[1]) * (xa[2]) - h));
curr.add(0, i[0], i[1] + 1, i[2] + 1,
fnq[0] * (dx[0] * (xa[1]) * (xa[2]) + h));

curr.add(1, i[0], i[1], i[2], fnq[1] * (dx[1] * (1.f - xa[2])));
curr.add(1, i[0], i[1], i[2] + 1, fnq[1] * (dx[1] * (xa[2])));

curr.add(2, i[0], i[1], i[2], fnq[2] * (dx[2] * (1.f - xa[1])));
curr.add(2, i[0], i[1] + 1, i[2], fnq[2] * (dx[2] * (xa[1])));
real_t h = (1.f / real_t(12.f)) * dx.prod();
Real3 fnq = qni_wni * fnqs_;

// clang-format off
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], fnq[0] * (dx[0] * ( xa[1]) * (1.f - xa[2]) - h));
curr.add(0, i[0], i[1], i[2] + 1, fnq[0] * (dx[0] * (1.f - xa[1]) * ( xa[2]) - h));
curr.add(0, i[0], i[1] + 1, i[2] + 1, fnq[0] * (dx[0] * ( xa[1]) * ( xa[2]) + h));

curr.add(1, i[0], i[1], i[2], fnq[1] * (dx[1] * (1.f - xa[2]) ));
curr.add(1, i[0], i[1], i[2] + 1, fnq[1] * (dx[1] * ( xa[2]) ));

curr.add(2, i[0], i[1], i[2], fnq[2] * (dx[2] * (1.f - xa[1]) ));
curr.add(2, i[0], i[1] + 1, i[2], fnq[2] * (dx[2] * ( xa[1]) ));
// clang-format on
}

private:
real3_t fnqs_;
Real3 fnqs_;
};

}; // namespace psc
66 changes: 28 additions & 38 deletions src/include/psc/deposit.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,13 +26,13 @@ class Deposit1st
{
public:
using real_t = R;
using real3_t = gt::sarray<R, 3>;
using Real3 = Vec3<real_t>;
using dim_t = D;

private:
template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& 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]);
Expand All @@ -41,8 +41,8 @@ private:
}

template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& 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]);
Expand All @@ -58,8 +58,7 @@ private:

public:
template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& 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{});
}
Expand All @@ -76,13 +75,13 @@ class Deposit2nd
{
public:
using real_t = R;
using real3_t = gt::sarray<R, 3>;
using Real3 = Vec3<R>;
using dim_t = D;

private:
template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& 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]);
Expand All @@ -105,8 +104,8 @@ private:
}

template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& 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]);
Expand Down Expand Up @@ -151,8 +150,7 @@ private:

public:
template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& 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{});
}
Expand All @@ -178,12 +176,10 @@ public:
static std::string suffix() { return "_1st_nc"; }

template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& ib,
const gt::sarray<T, 3>& x, T val)

void operator()(F& flds, const Int3& ib, const Vec3<T>& x, T val)
{
gt::sarray<int, 3> l;
gt::sarray<T, 3> h;
Int3 l;
Vec3<T> h;
for (int d = 0; d < 3; d++) {
l[d] = fint(x[d]);
h[d] = x[d] - l[d];
Expand All @@ -201,12 +197,10 @@ public:
static std::string suffix() { return "_1st_cc"; }

template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& ib,
const gt::sarray<T, 3>& x, T val)

void operator()(F& flds, const Int3& ib, const Vec3<T>& x, T val)
{
gt::sarray<int, 3> l;
gt::sarray<T, 3> h;
Int3 l;
Vec3<T> h;
for (int d = 0; d < 3; d++) {
l[d] = fint(x[d] - .5f);
h[d] = x[d] - .5f - l[d];
Expand All @@ -228,12 +222,10 @@ public:
static std::string suffix() { return "_2nd_nc"; }

template <typename F>
void operator()(F& flds, const gt::sarray<int, 3>& ib,
const gt::sarray<T, 3>& x, T val)

void operator()(F& flds, const Int3& ib, const Vec3<T>& x, T val)
{
gt::sarray<int, 3> l;
gt::sarray<T, 3> h;
Int3 l;
Vec3<T> h;
for (int d = 0; d < 3; d++) {
l[d] = nint(x[d]);
h[d] = l[d] - x[d]; // negated!
Expand All @@ -245,14 +237,14 @@ public:
};

template <typename D, typename F, typename T>
void nc(F& flds, const gt::sarray<int, 3>& ib, const gt::sarray<T, 3>& x, T val)
void nc(F& flds, const Int3& ib, const Vec3<T>& x, T val)
{
Deposit1stNc<T, D> deposit;
deposit(flds, ib, x, val);
}

template <typename D, typename F, typename T>
void cc(F& flds, const gt::sarray<int, 3>& ib, const gt::sarray<T, 3>& x, T val)
void cc(F& flds, const Int3& ib, const Vec3<T>& x, T val)
{
Deposit1stCc<T, D> deposit;
deposit(flds, ib, x, val);
Expand All @@ -276,27 +268,25 @@ class Deposit
public:
using real_t = R;
using dim_t = D;
using real3_t = gt::sarray<real_t, 3>;
using Real3 = Vec3<real_t>;
using DepositNorm = DEPOSIT_NORM<real_t, dim_t>;
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 <typename F>
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_;
};

Expand Down
Loading