Interpolate Current

src/FEM/GridPathSegmenter.h

@@ -0,0 +1,29 @@
+#ifndef IPPL_GRIDPATH_SEGMENTER_H
+#define IPPL_GRIDPATH_SEGMENTER_H
+#include <array>
+
+namespace ippl {
+
+template<unsigned Dim, typename T>
+struct Segment {
+  ippl::Vector<T,Dim> p0, p1;
+};
+
+struct DefaultCellCrossingRule {};
+
+template<unsigned Dim, typename T, typename Rule = DefaultCellCrossingRule>
+struct GridPathSegmenter {
+
+  static KOKKOS_INLINE_FUNCTION
+  std::array<Segment<Dim,T>, Dim+1>
+  split(const ippl::Vector<T, Dim>& A,
+        const ippl::Vector<T, Dim>& B,
+        const ippl::Vector<T, Dim>& origin,
+        const ippl::Vector<T, Dim>& h);
+};
+
+} // namespace ippl 
+
+#include "GridPathSegmenter.hpp"
+
+#endif

src/FEM/GridPathSegmenter.hpp

@@ -0,0 +1,123 @@
+#include <array>
+#include <cmath>
+#include <type_traits>
+
+namespace ippl {
+
+template<int I, int End, class F>
+KOKKOS_INLINE_FUNCTION
+constexpr void static_for(F&& f) {
+  if constexpr (I < End) {
+    f(std::integral_constant<int, I>{});
+    static_for<I+1, End>(std::forward<F>(f));
+  }
+}
+
+template<typename T>
+KOKKOS_INLINE_FUNCTION
+void sort2(T& a, T& b) { if (a > b) { T t=a; a=b; b=t; } }
+
+template<typename T>
+KOKKOS_INLINE_FUNCTION
+void sort3(T& a, T& b, T& c) { sort2(a,b); sort2(b,c); sort2(a,b); }
+
+template<unsigned Dim, typename T>
+KOKKOS_INLINE_FUNCTION
+ippl::Vector<T,Dim> lerp_point(const ippl::Vector<T,Dim>& A,
+                             const ippl::Vector<T,Dim>& B, T t) {
+  ippl::Vector<T,Dim> out{}; 
+  for (unsigned a=0; a<Dim; ++a) out[a] = A[a] + (B[a]-A[a]) * t;
+  return out;
+}
+
+// ---------- per-axis cut times (for DefaultCellCrossingRule) ----------
+template<unsigned Dim, typename T>
+struct CutTimes { std::array<T,Dim> t; };
+
+template<unsigned Dim, typename T>
+KOKKOS_INLINE_FUNCTION
+CutTimes<Dim,T> compute_axis_cuts_default(
+  const ippl::Vector<T,Dim>& A,
+  const ippl::Vector<T,Dim>& B,
+  const ippl::Vector<T,Dim>& origin,
+  const ippl::Vector<T,Dim>& h)
+{
+  CutTimes<Dim,T> cuts;
+  for (unsigned a=0; a<Dim; ++a) cuts.t[a] = (T)2; // sentinel (>1)
+
+  const T eps1 = (T)1e-12, one = (T)1;
+
+  auto axis_cut = [&](auto Ax) {
+    constexpr int a = Ax;
+    const T d   = B[a] - A[a];
+    const T eps = eps1 * h[a];
+    if (std::fabs(d) <= eps) return; // no motion → no crossing
+
+    const T nA  = (A[a] - origin[a]) / h[a]; // in cell units
+    // nearest plane index in direction toward B; bias off-plane a hair
+    const T k   = (d > 0) ? std::ceil(nA - eps1) : std::floor(nA + eps1);
+    const T pa  = origin[a] + k * h[a];
+    const T t   = (pa - A[a]) / d;
+
+    if (t > eps1 && t < one - eps1) cuts.t[a] = t;
+  };
+
+  static_for<0,Dim>(axis_cut);
+  return cuts;
+}
+
+// ---------- endpoints builder: [0, t_sorted..., 1] (keeps duplicates) ----------
+template<unsigned Dim, typename T>
+KOKKOS_INLINE_FUNCTION
+void make_endpoints_fixed(const CutTimes<Dim,T>& cuts,
+                          std::array<T,Dim+2>& Tcuts /*out*/)
+{
+  T t0 = cuts.t[0];
+  T t1 = (Dim>=2) ? cuts.t[1] : (T)2;
+  T t2 = (Dim==3) ? cuts.t[2] : (T)2;
+
+  if constexpr (Dim==2) {
+    sort2(t0, t1);
+  } else { // Dim==3
+    sort3(t0, t1, t2);
+  }
+
+  const T one = (T)1;
+  auto clamp_or_one = [&](T v)->T { return (v >= (T)1.5) ? one : (v < one ? v : one); };
+
+  Tcuts[0] = (T)0;
+  Tcuts[1] = clamp_or_one(t0);
+  if constexpr (Dim>=2) Tcuts[2] = clamp_or_one(t1);
+  if constexpr (Dim==3) Tcuts[3] = clamp_or_one(t2);
+  Tcuts[Dim+1] = one;
+}
+
+// ---------------------------------------------------------------------------------
+// DefaultCellCrossingRule specialization 
+// ---------------------------------------------------------------------------------
+
+template<unsigned Dim, typename T, typename Rule>
+KOKKOS_INLINE_FUNCTION
+std::array<Segment<Dim,T>, Dim+1>
+GridPathSegmenter<Dim,T,Rule>::split(
+    const ippl::Vector<T,Dim>& A,
+    const ippl::Vector<T,Dim>& B,
+    const ippl::Vector<T,Dim>& origin,
+    const ippl::Vector<T,Dim>& h)
+{
+  const auto cuts = compute_axis_cuts_default<Dim,T>(A,B,origin,h);
+
+  std::array<T,Dim+2> Tcuts{};
+  make_endpoints_fixed<Dim,T>(cuts, Tcuts);
+
+  std::array<Segment<Dim,T>, Dim+1> segs{};
+  for (unsigned i = 0; i < Dim + 1; ++i) {
+    const T ta = Tcuts[i];
+    const T tb = Tcuts[i+1];
+    segs[i].p0 = lerp_point<Dim,T>(A,B,ta);
+    segs[i].p1 = lerp_point<Dim,T>(A,B,tb);
+  }
+  return segs;
+}
+
+} // namespace ippl 

src/FEM/ProjectCurrent.hpp

@@ -0,0 +1,42 @@
+#ifndef IPPL_PROJECT_CURRENT_H
+#define IPPL_PROJECT_CURRENT_H
+
+namespace ippl {
+
+
+template <typename Mesh,
+          typename ChargeAttrib,
+          typename PosAttrib,
+          typename FEMVector,
+          typename NedelecSpace,
+          typename policy_type = Kokkos::RangePolicy<>>
+inline void assemble_current_whitney1(const Mesh& mesh,
+                                      const ChargeAttrib& q_attrib,
+                                      const PosAttrib& X0,
+                                      const PosAttrib& X1,
+                                      FEMVector& fem_vector,
+                                      const NedelecSpace& space,
+                                      policy_type iteration_policy)
+{
+  using T = Mesh::value_type;
+
+  const auto origin = mesh.getOrigin();
+  const auto h   = mesh.getMeshSpacing();
+  constexpr unsigned Dim = Mesh::Dimension;
+
+  Kokkos::parallel_for("assemble_current_whitney1_make_segments", iteration_policy,
+    KOKKOS_LAMBDA(const std::size_t p) {
+
+    auto segs = ippl::GridPathSegmenter<Dim, T, ippl::DefaultCellCrossingRule>
+                    ::split(X0(p), X1(p), origin, h);
+
+
+
+
+  });
+
+
+}
+
+}
+#endif

src/IpplCore.h

@@ -56,5 +56,7 @@
 
 // FEM Operators
 #include "FEM/FEMInterpolate.hpp"
+#include "FEM/GridPathSegmenter.h"
+#include "FEM/ProjectCurrent.hpp"
 
 #endif

unit_tests/FEM/AssembleCurrentRHS.cpp

@@ -0,0 +1,171 @@
+#include "Ippl.h"
+#include "TestUtils.h"
+#include "gtest/gtest.h"
+
+template <class PLayout>
+struct Bunch : public ippl::ParticleBase<PLayout> {
+    Bunch(PLayout& playout)
+        : ippl::ParticleBase<PLayout>(playout) {
+        this->addAttribute(Q);
+    }
+
+    ~Bunch() {}
+
+    typedef ippl::ParticleAttrib<double> charge_container_type;
+    charge_container_type Q;
+};
+
+
+
+template <unsigned Dim, typename T> struct ElemSelector;
+
+template <typename T>
+struct ElemSelector<1, T> {
+  using Elem = ippl::EdgeElement<T>;
+  using Quad = ippl::MidpointQuadrature<T, 1, Elem>;
+
+  static Elem make_elem() { return Elem{}; }
+  static Quad make_quad(const Elem& e) { return Quad(e); }
+};
+
+template <typename T>
+struct ElemSelector<2, T> {
+  using Elem = ippl::QuadrilateralElement<T>;
+  using Quad = ippl::MidpointQuadrature<T, 1, Elem>;
+
+  static Elem make_elem() { return Elem{}; }
+  static Quad make_quad(const Elem& e) { return Quad(e); }
+};
+
+template <typename T>
+struct ElemSelector<3, T> {
+  using Elem = ippl::HexahedralElement<T>;
+  using Quad = ippl::MidpointQuadrature<T, 1, Elem>;
+
+  static Elem make_elem() { return Elem{}; }
+  static Quad make_quad(const Elem& e) { return Quad(e); }
+};
+
+
+template <typename>
+class AssembleCurrentTest;
+
+template <typename T, unsigned Dim>
+class AssembleCurrentTest<Parameters<T, Rank<Dim>>> : public ::testing::Test {
+public:
+  using value_type = T;
+  static constexpr unsigned dim = Dim;
+
+  using Mesh_t      = ippl::UniformCartesian<T, dim>;
+  using Centering_t = typename Mesh_t::DefaultCentering;
+
+  using ElemSel     = ElemSelector<dim, T>;
+  using Elem        = typename ElemSel::Elem;
+  using Quad        = typename ElemSel::Quad;
+  // Somehow the field type is completely irrelevant for the fem space??
+  using FieldType   = ippl::Field<T, Dim, Mesh_t, typename Mesh_t::DefaultCentering>;
+  using Layout      = ippl::FieldLayout<Dim>;
+
+  using NedelecType = ippl::NedelecSpace<T, Dim, 1, Elem, Quad, FieldType>;
+
+  using playout_t   = ippl::ParticleSpatialLayout<T, dim>;
+  using bunch_t     = Bunch<playout_t>;
+
+  static ippl::NDIndex<dim> make_owned_nd(int nx) {
+    ippl::Index I0(nx);
+    if constexpr (dim == 1)
+      return ippl::NDIndex<1>(I0);
+    else if constexpr (dim == 2)
+      return ippl::NDIndex<2>(I0, I0);
+    else
+      return ippl::NDIndex<3>(I0, I0, I0);
+  }
+
+  static Layout make_layout(const ippl::NDIndex<dim>& owned) {
+    std::array<bool, dim> par{}; par.fill(true);
+    return ippl::FieldLayout<dim>(MPI_COMM_WORLD, owned, par);
+  }
+
+  static Mesh_t make_mesh(const ippl::NDIndex<dim>& owned,
+                          const ippl::Vector<T, dim>& h,
+                          const ippl::Vector<T, dim>& origin) {
+    return Mesh_t(owned, h, origin);
+  }
+
+  static NedelecType make_space(Mesh_t& mesh, Layout l) {
+    Elem e = ElemSel::make_elem();
+    Quad q = ElemSel::make_quad(e);
+    return NedelecType(mesh, e, q, l);
+  }
+};
+
+using Precisions = TestParams::Precisions;
+using Ranks      = TestParams::Ranks<2, 3>;
+using Tests      = TestForTypes<CreateCombinations<Precisions, Ranks>::type>::type;
+
+TYPED_TEST_SUITE(AssembleCurrentTest, Tests);
+
+// ------------------ Actual minimal test ------------------
+
+TYPED_TEST(AssembleCurrentTest, SingleParticle_Smoke) {
+  using T   = typename TestFixture::value_type;
+  constexpr unsigned Dim = TestFixture::dim;
+
+  using bunch_t    = typename TestFixture::bunch_t;
+  using playout_t  = typename TestFixture::playout_t;
+
+  int nx = 4;
+
+  ippl::Vector<T, Dim> origin(0.0);
+  ippl::Vector<T, Dim> h(1.0);
+
+  auto owned  = TestFixture::make_owned_nd(nx);
+  auto layout = TestFixture::make_layout(owned);
+  auto mesh   = TestFixture::make_mesh(owned, h, origin);
+  auto space  = TestFixture::make_space(mesh, layout);
+
+  playout_t playout(layout, mesh);
+  bunch_t   bunch(playout);
+
+  // --- create 1 particle ---
+  bunch.create(1);
+  {
+    auto R_host = bunch.R.getHostMirror();
+    auto Q_host = bunch.Q.getHostMirror();
+    for (unsigned d=0; d<Dim; ++d)
+      R_host(0)[d] = T(0.5) * h[d];  // inside domain
+    Q_host(0) = 1.0;
+    Kokkos::deep_copy(bunch.R.getView(), R_host);
+    Kokkos::deep_copy(bunch.Q.getView(), Q_host);
+    bunch.update();
+  }
+
+  // --- mock "next" positions (slightly moved) ---
+  bunch_t bunch_next(playout);
+  bunch_next.create(1);
+  {
+    auto Rn_host = bunch_next.R.getHostMirror();
+    for (unsigned d=0; d<Dim; ++d)
+      Rn_host(0)[d] = T(0.6) * h[d];
+    Kokkos::deep_copy(bunch_next.R.getView(), Rn_host);
+    bunch_next.update();
+  }
+
+  auto policy = Kokkos::RangePolicy<>(0, bunch.getLocalNum());
+  auto fem_vector = space.createFEMVector();
+
+  // Call function (no assertions yet)
+  ippl::assemble_current_whitney1(mesh, bunch.Q, bunch.R, bunch_next.R, fem_vector, space, policy);
+
+
+  SUCCEED() << "assemble_current_whitney1 ran without error for 1 particle in "
+            << Dim << "D.";
+}
+
+int main(int argc, char* argv[]) {
+  ippl::initialize(argc, argv);
+  ::testing::InitGoogleTest(&argc, argv);
+  int result = RUN_ALL_TESTS();
+  ippl::finalize();
+  return result;
+}

unit_tests/FEM/CMakeLists.txt

@@ -1,6 +1,8 @@
 file(RELATIVE_PATH _relPath "${PROJECT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}")
 message(STATUS "Adding unit tests found in ${_relPath}")
 
+add_ippl_test(AssembleCurrentRHS)
+add_ippl_test(GenerateSegments)
 add_ippl_test(EdgeElement)
 add_ippl_test(QuadrilateralElement)
 add_ippl_test(HexahedralElement)