From 3b033a37666c20ecba761db85bc31d1c14711513 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 21 Oct 2025 12:17:04 -0400 Subject: [PATCH 01/15] rng: rename params --- src/include/rng.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/include/rng.hxx b/src/include/rng.hxx index 9b8c7236e..56155baee 100644 --- a/src/include/rng.hxx +++ b/src/include/rng.hxx @@ -81,11 +81,11 @@ struct Normal Normal() : Normal(0, 1) {} Real get() { return dist(*gen); } - // FIXME remove me, or make standalone func - Real get(Real mean, Real stdev) + + Real get(Real mean_override, Real stdev_override) { // should be possible to pass params to existing dist - return std::normal_distribution(mean, stdev)(*gen); + return std::normal_distribution(mean_override, stdev_override)(*gen); } private: From d203893ea1e71d99c07d03f73da9a869c6407c3f Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 21 Oct 2025 12:17:27 -0400 Subject: [PATCH 02/15] rng: +Uniform::get overload --- src/include/rng.hxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/include/rng.hxx b/src/include/rng.hxx index 56155baee..2842c5710 100644 --- a/src/include/rng.hxx +++ b/src/include/rng.hxx @@ -59,6 +59,12 @@ struct Uniform Real get() { return dist(*gen); } + Real get(Real min_override, Real max_override) + { + return std::uniform_real_distribution(min_override, + max_override)(*gen); + } + private: detail::GeneratorPtr gen; std::uniform_real_distribution dist; From abcc7eef3cf95e20ff140aaadffc7154dabc2c96 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 30 Oct 2025 15:03:57 -0400 Subject: [PATCH 03/15] Vec3: +util methods --- src/kg/include/kg/Vec3.h | 59 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index aff5ba926..4c842ef83 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -9,6 +9,7 @@ #include "cuda_compat.h" #include #include +#include "psc_bits.h" namespace kg { @@ -102,6 +103,64 @@ struct Vec : gt::sarray return *this; } + KG_INLINE T sum() const + { + T sum{0}; + for (size_t i = 0; i < N; i++) { + sum += (*this)[i]; + } + return sum; + } + + KG_INLINE T prod() const + { + T prod{1}; + for (size_t i = 0; i < N; i++) { + prod *= (*this)[i]; + } + return prod; + } + + KG_INLINE T dot(const Vec& w) const + { + T sum{0}; + for (size_t i = 0; i < N; i++) { + sum += (*this)[i] * w[i]; + } + return sum; + } + + KG_INLINE Vec cross(const Vec& w) const + { + return {(*this)[1] * w[2] - (*this)[2] * w[1], + (*this)[2] * w[0] - (*this)[0] * w[2], + (*this)[0] * w[1] - (*this)[1] * w[0]}; + } + + KG_INLINE T mag2() const { return this->dot(this); } + + KG_INLINE T mag() const { return sqrt(this->mag2()); } + + KG_INLINE T max() const + { + auto max = (*this)[0]; + for (size_t i = 1; i < N; i++) { + max = std::max(max, (*this)[i]); + } + return max; + } + + KG_INLINE Vec inv() const { return T(1) / *this; } + + KG_INLINE Vec fint() const + { + Vec res; + for (size_t i = 0; i < N; i++) { + res[i] = ::fint((*this)[i]); + } + return res; + } + // conversion to pointer KG_INLINE operator const T*() const { return this->data(); } From 0f09ddfbd4bf031d06c51430a6e05a79c724c55c Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 13:59:05 -0500 Subject: [PATCH 04/15] Vec3: +missing div operators --- src/kg/include/kg/Vec3.h | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index 4c842ef83..b8cf3362a 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -103,6 +103,14 @@ struct Vec : gt::sarray return *this; } + KG_INLINE Vec& operator/=(T s) + { + for (size_t i = 0; i < N; i++) { + (*this)[i] /= s; + } + return *this; + } + KG_INLINE T sum() const { T sum{0}; @@ -212,6 +220,24 @@ KG_INLINE Vec operator*(T s, const Vec& v) return res; } +template +KG_INLINE Vec operator/(T s, const Vec& v) +{ + Vec res; + for (size_t i = 0; i < N; i++) { + res[i] = s / v[i]; + } + return res; +} + +template +KG_INLINE Vec operator/(const Vec& v, T s) +{ + Vec res = v; + res /= s; + return res; +} + template KG_INLINE Vec operator/(const Vec& v, const Vec& w) { From 67a71ded7a5f0a363a83a6571ae3437549542604 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 10:59:58 -0500 Subject: [PATCH 05/15] dim: constexpr is_invar --- src/include/dim.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/dim.hxx b/src/include/dim.hxx index 089097a27..8a076d0f8 100644 --- a/src/include/dim.hxx +++ b/src/include/dim.hxx @@ -22,7 +22,7 @@ struct Invar return {!InvarX::value, !InvarY::value, !InvarZ::value}; } - static bool is_invar(int dim) + static constexpr bool is_invar(int dim) { switch (dim) { case 0: return InvarX::value; From 4ce3384e6dbd4c7a1efbc4ffac5108b99a77d9de Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 11:03:42 -0500 Subject: [PATCH 06/15] dim: todo consteval --- src/include/dim.hxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/include/dim.hxx b/src/include/dim.hxx index 8a076d0f8..4e4453728 100644 --- a/src/include/dim.hxx +++ b/src/include/dim.hxx @@ -22,6 +22,7 @@ struct Invar return {!InvarX::value, !InvarY::value, !InvarZ::value}; } + // TODO c++20 make consteval? static constexpr bool is_invar(int dim) { switch (dim) { From 8743326b5a36db66717fbd2b0e1225fb2336c5cf Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 31 Dec 2025 15:41:38 -0500 Subject: [PATCH 07/15] push_particles_1vb: use inv() --- src/libpsc/psc_push_particles/push_particles_1vb.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index 912ad48ad..32d32e337 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -24,7 +24,7 @@ struct PushParticlesVb static void push_mprts(Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); - Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx); + Real3 dxi = grid.domain.dx.inv(); real_t dq_kind[MAX_NR_KINDS]; auto& kinds = grid.kinds; assert(kinds.size() <= MAX_NR_KINDS); From 2234ac0c90a8c3ad55c6dd320b9e1707a700fd10 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 6 Jan 2026 11:51:11 -0500 Subject: [PATCH 08/15] push_particles_1vb: -find_idx_pos_1st_rel --- src/libpsc/psc_push_particles/push_particles_1vb.hxx | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index 32d32e337..a09841a43 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -60,10 +60,8 @@ struct PushParticlesVb auto v = advance.calc_v(prt.u()); advance.push_x(prt.x(), v); - Int3 final_index; - Real3 final_pos_normalized; - find_idx_pos_1st_rel(prt.x(), dxi, final_index, final_pos_normalized, - real_t(0.)); + Real3 final_pos_normalized = prt.x() * dxi; + Int3 final_index = final_pos_normalized.fint(); // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt int lg[3]; From 599b35ab46a58d7fbd47d65000291609c80ae54a Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 8 Jan 2026 12:02:30 -0500 Subject: [PATCH 09/15] push_particles_1vb: manually convert to Real3 --- src/libpsc/psc_push_particles/push_particles_1vb.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index a09841a43..18c8ff180 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -24,7 +24,7 @@ struct PushParticlesVb static void push_mprts(Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); - Real3 dxi = grid.domain.dx.inv(); + Real3 dxi = Real3(grid.domain.dx).inv(); real_t dq_kind[MAX_NR_KINDS]; auto& kinds = grid.kinds; assert(kinds.size() <= MAX_NR_KINDS); From b3b6ba9419f5a548474adb6b0ae372013848dcca Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 14 Jan 2026 15:06:06 -0500 Subject: [PATCH 10/15] Vec3: fix mag2 --- src/kg/include/kg/Vec3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index b8cf3362a..19bfc073b 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -145,7 +145,7 @@ struct Vec : gt::sarray (*this)[0] * w[1] - (*this)[1] * w[0]}; } - KG_INLINE T mag2() const { return this->dot(this); } + KG_INLINE T mag2() const { return this->dot(*this); } KG_INLINE T mag() const { return sqrt(this->mag2()); } From 50b37d94e7fed4a6e8075f0d21a3289682fe538a Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 14 Jan 2026 15:16:25 -0500 Subject: [PATCH 11/15] Vec3: +min --- src/kg/include/kg/Vec3.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index 19bfc073b..94401ff6f 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -151,11 +151,12 @@ struct Vec : gt::sarray KG_INLINE T max() const { - auto max = (*this)[0]; - for (size_t i = 1; i < N; i++) { - max = std::max(max, (*this)[i]); - } - return max; + return *std::max_element(this->begin(), this->end()); + } + + KG_INLINE T min() const + { + return *std::min_element(this->begin(), this->end()); } KG_INLINE Vec inv() const { return T(1) / *this; } From de25be04e721f08a1696f51a5f42cabd63736757 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 14 Jan 2026 15:33:08 -0500 Subject: [PATCH 12/15] Vec3: +min, max with arg --- src/kg/include/kg/Vec3.h | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index 94401ff6f..2c4beee8b 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -159,6 +159,24 @@ struct Vec : gt::sarray return *std::min_element(this->begin(), this->end()); } + static KG_INLINE Vec max(const Vec& v, const Vec& w) + { + Vec res; + for (int i = 0; i < N; i++) { + res[i] = std::max(v[i], w[i]); + } + return res; + } + + static KG_INLINE Vec min(const Vec& v, const Vec& w) + { + Vec res; + for (int i = 0; i < N; i++) { + res[i] = std::min(v[i], w[i]); + } + return res; + } + KG_INLINE Vec inv() const { return T(1) / *this; } KG_INLINE Vec fint() const From 7e525d07f36c6dc81c57253ef6d124c86b86a207 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 14 Jan 2026 15:38:49 -0500 Subject: [PATCH 13/15] Vec3: +scalar +, - ops --- src/kg/include/kg/Vec3.h | 41 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index 2c4beee8b..907755ef5 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -71,6 +71,14 @@ struct Vec : gt::sarray return *this; } + KG_INLINE Vec& operator+=(T s) + { + for (size_t i = 0; i < N; i++) { + (*this)[i] += s; + } + return *this; + } + KG_INLINE Vec& operator-=(const Vec& w) { for (size_t i = 0; i < N; i++) { @@ -79,6 +87,14 @@ struct Vec : gt::sarray return *this; } + KG_INLINE Vec& operator-=(T s) + { + for (size_t i = 0; i < N; i++) { + (*this)[i] -= s; + } + return *this; + } + KG_INLINE Vec& operator*=(const Vec& w) { for (size_t i = 0; i < N; i++) { @@ -215,6 +231,14 @@ KG_INLINE Vec operator+(const Vec& v, const Vec& w) return res; } +template +KG_INLINE Vec operator+(T s, const Vec& v) +{ + Vec res = v; + res += s; + return res; +} + template KG_INLINE Vec operator-(const Vec& v, const Vec& w) { @@ -223,6 +247,23 @@ KG_INLINE Vec operator-(const Vec& v, const Vec& w) return res; } +template +KG_INLINE Vec operator-(T s, const Vec& v) +{ + Vec res; + for (size_t i = 0; i < N; i++) { + res[i] = s - v[i]; + } + return res; +} +template +KG_INLINE Vec operator-(const Vec& v, T s) +{ + Vec res = v; + res -= s; + return res; +} + template KG_INLINE Vec operator*(const Vec& v, const Vec& w) { From 47ecfb2408f1388632688001fbbb444bd95401db Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 14 Jan 2026 16:40:40 -0500 Subject: [PATCH 14/15] +VecRange --- src/kg/include/kg/VecRange.hxx | 80 ++++++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 src/kg/include/kg/VecRange.hxx diff --git a/src/kg/include/kg/VecRange.hxx b/src/kg/include/kg/VecRange.hxx new file mode 100644 index 000000000..e8590ea55 --- /dev/null +++ b/src/kg/include/kg/VecRange.hxx @@ -0,0 +1,80 @@ +#pragma once + +#include "Vec3.h" + +/** + * @brief Facilitates iteration over `N`-dimensional spaces in row-major order. + * + * For example, + * ```c++ + * for (Int3 i3 : VecRange(mins, maxs)) { + * // ... + * } + * ``` + * is equivalent to + * ```c++ + * for (int ix = mins[0]; ix < maxs[0]; ix++) { + * for (int iy = mins[1]; iy < maxs[1]; iy++) { + * for (int iz = mins[2]; iz < maxs[2]; iz++) { + * Int3 i3{ix, iy, iz}; + * // ... + * } + * } + * } + * ``` + * @tparam T scalar type (e.g. `int`) + * @tparam N number of dimensions + */ +template +class VecRange +{ +public: + using Vec = kg::Vec; + + class Iterator + { + public: + Iterator(const Vec& starts, const Vec& stops, bool end = false) + : starts_(starts), stops_(stops), current_(starts), end_(end) + {} + + Iterator& operator++() + { + for (int d = N - 1; d >= 0; d--) { + ++current_[d]; + + if (current_[d] < stops_[d]) { + break; + } else { + current_[d] = starts_[d]; + + if (d == 0) { + end_ = true; + } + } + } + return *this; + } + + Vec operator*() { return current_; } + + bool operator!=(const Iterator& other) const { return end_ != other.end_; } + + private: + const Vec& starts_; + const Vec& stops_; + Vec current_; + bool end_; + }; + + VecRange(const Vec& starts, const Vec& stops) : starts_(starts), stops_(stops) + {} + + Iterator begin() const { return Iterator(starts_, stops_, false); } + + Iterator end() const { return Iterator(starts_, stops_, true); } + +private: + const Vec starts_; + const Vec stops_; +}; \ No newline at end of file From 69c0f0d3f20f7778b2bef072ae3e1373a17504b2 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 15 Jan 2026 09:34:47 -0500 Subject: [PATCH 15/15] Vec3: add comments --- src/kg/include/kg/Vec3.h | 57 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index 907755ef5..b77395e3b 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -127,6 +127,10 @@ struct Vec : gt::sarray return *this; } + /** + * @brief Calculates the sum of this vec's elements. + * @return the sum + */ KG_INLINE T sum() const { T sum{0}; @@ -136,6 +140,10 @@ struct Vec : gt::sarray return sum; } + /** + * @brief Calculates the product of this vec's elements. + * @return the product + */ KG_INLINE T prod() const { T prod{1}; @@ -145,6 +153,11 @@ struct Vec : gt::sarray return prod; } + /** + * @brief Calculates the dot product of this and another vec. + * @param w the other vec + * @return the dot product + */ KG_INLINE T dot(const Vec& w) const { T sum{0}; @@ -154,6 +167,12 @@ struct Vec : gt::sarray return sum; } + /** + * @brief Calculates the cross product of this and another vec. This method + * assumes `N` is 3. + * @param w the right-hand-side vec + * @return the cross product + */ KG_INLINE Vec cross(const Vec& w) const { return {(*this)[1] * w[2] - (*this)[2] * w[1], @@ -161,20 +180,43 @@ struct Vec : gt::sarray (*this)[0] * w[1] - (*this)[1] * w[0]}; } + /** + * @brief Calculates the square of the magnitude of this vec. + * @return the squared magnitude + */ KG_INLINE T mag2() const { return this->dot(*this); } + /** + * @brief Calculates the magnitude of this vec by taking the square root of + * the squared magnitude. + * @return the magnitude + */ KG_INLINE T mag() const { return sqrt(this->mag2()); } + /** + * @brief Determines the maximal value of any of this vec's elements. + * @return the maximal value + */ KG_INLINE T max() const { return *std::max_element(this->begin(), this->end()); } + /** + * @brief Determines the minimal value of any of this vec's elements. + * @return the minimal value + */ KG_INLINE T min() const { return *std::min_element(this->begin(), this->end()); } + /** + * @brief Calculates the elementwise maximum between two vecs. + * @param v the first vec + * @param w the second vec + * @return a vec containing the maximal values + */ static KG_INLINE Vec max(const Vec& v, const Vec& w) { Vec res; @@ -184,6 +226,12 @@ struct Vec : gt::sarray return res; } + /** + * @brief Calculates the elementwise minimum between two vecs. + * @param v the first vec + * @param w the second vec + * @return a vec containing the minimal values + */ static KG_INLINE Vec min(const Vec& v, const Vec& w) { Vec res; @@ -193,8 +241,17 @@ struct Vec : gt::sarray return res; } + /** + * @brief Calculates the (multiplicative) inverse of each of this vec's + * elements. + * @return a vec of the inverses + */ KG_INLINE Vec inv() const { return T(1) / *this; } + /** + * @brief Calculates the floor of each of this vec's elements. + * @return an integer vec of the floors + */ KG_INLINE Vec fint() const { Vec res;