Interpolate Current #437

JoMee · 2025-10-28T14:34:02Z

This PR adds the segment splitting as well as the current state of the interpolation scheme.

…efine that function

mohsensadr · 2025-10-29T14:51:10Z

src/FEM/GridPathSegmenter.h

+
+  static KOKKOS_INLINE_FUNCTION
+  std::array<Segment<Dim,T>, Dim+1>
+  split(const ippl::Vector<T, Dim>& A,


Isn't here the templates of GridPathSegmenter missing?

Shouldn't it be

GridPathSegmenter<Dim,T,Rule>::split(

instead?

In the implementation .hpp file, you have

template<unsigned Dim, typename T, typename Rule> KOKKOS_INLINE_FUNCTION std::array<Segment<Dim,T>, Dim+1> GridPathSegmenter<Dim,T,Rule>::split(

mohsensadr · 2025-10-29T14:51:17Z

src/FEM/GridPathSegmenter.hpp

+
+template<typename T>
+KOKKOS_INLINE_FUNCTION
+void sort2(T& a, T& b) { if (a > b) { T t=a; a=b; b=t; } }


would it make sense to allow compile-time evaluation with constexpr? same for sort3, and lerp_point?

mohsensadr · 2025-10-29T14:51:22Z

src/FEM/GridPathSegmenter.hpp

+void make_endpoints_fixed(const CutTimes<Dim,T>& cuts,
+                          std::array<T,Dim+2>& Tcuts /*out*/)
+{
+  T t0 = cuts.t[0];


what if Dim==1? If that's never the case, then consider having a check. Because with Dim=1, I think the else in line 82 would be out of bound.

mohsensadr · 2025-10-29T14:52:33Z

src/FEM/ProjectCurrent.hpp

+                                      const NedelecSpace& space,
+                                      policy_type iteration_policy)
+{
+  using T = Mesh::value_type;


Shouldn't it be?

using T = typename Mesh::value_type;

mohsensadr · 2025-10-29T14:59:26Z

src/FEM/ProjectCurrent.hpp

+  const auto h   = mesh.getMeshSpacing();
+  constexpr unsigned Dim = Mesh::Dimension;
+
+  Kokkos::parallel_for("assemble_current_whitney1_make_segments", iteration_policy,


Are we sure that h, origin, etc used in the Kokkos loop are not references to host-only memory?

Also, there are too many empty lines here somehow.

s-mayani · 2025-10-30T10:00:20Z

src/FEM/GridPathSegmenter.h

+
+template<unsigned Dim, typename T>
+struct Segment {
+  ippl::Vector<T,Dim> p0, p1;


Since you're in namespace ippl, you don't need to specify ippl:: each time for the Vector<T,Dim> type. Comment applies to all instances of this in the file.

s-mayani · 2025-10-30T10:07:05Z

src/FEM/GridPathSegmenter.hpp

+    f(std::integral_constant<int, I>{});
+    static_for<I+1, End>(std::forward<F>(f));


Will these std:: functions work on GPU? Or does it not matter since it's marked constexpr?

I also don't think that standard C++'s STL, <cmath>, <algorithm>, etc run on the GPUs, unless proven otherwise :) One has to be very careful, and double check CUDA's manual for each of them.

I see your point on constexpr, that it pushes the compiler to compute the expressions on the host at compile time, and therefore does not compute it on GPU, so it may not cause any problem as it is not deployed on the GPU. But is that for sure? I see comments saying that expressions marked with constexpr maybe are computed at compile time depending on compiler being able to find all the known values at the compiled time or not. Otherwise, it is a runtime calculation regardless of having constexpr.

If this code is supposed to run on the host only, then why do we have KOKKOS_INLINE_FUNCTION?

s-mayani · 2025-10-30T10:17:14Z

src/FEM/GridPathSegmenter.hpp

+
+// ---------- per-axis cut times (for DefaultCellCrossingRule) ----------
+template<unsigned Dim, typename T>
+struct CutTimes { std::array<T,Dim> t; };


Is this array only going to be accessed on host? I don't think std::array is accessible within Kokkos kernels on device. One would need to use: https://kokkos.org/kokkos-core-wiki/API/core/stl-compat/Array.html, or ippl Vectors.

s-mayani · 2025-10-30T10:17:42Z

src/FEM/GridPathSegmenter.hpp

+  const ippl::Vector<T,Dim>& A,
+  const ippl::Vector<T,Dim>& B,
+  const ippl::Vector<T,Dim>& origin,
+  const ippl::Vector<T,Dim>& h)


again, since we are in the ippl namespace, no need for ippl::

s-mayani · 2025-10-30T10:28:53Z

src/FEM/GridPathSegmenter.hpp

+
+  auto axis_cut = [&](auto Ax) {
+    constexpr int a = Ax;
+    const T d   = B[a] - A[a];


Maybe rename d to dist to avoid confusion with d = Dim (which is the case in many places in ippl)?

s-mayani · 2025-10-30T10:29:57Z

src/FEM/GridPathSegmenter.hpp

+
+    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);


std::ceil and std::floor should be replaced by Kokkos::ceil and Kokkos::floor.

s-mayani · 2025-10-30T10:32:46Z

src/FEM/GridPathSegmenter.hpp

+    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;


add a comment explaining why we need this here

s-mayani · 2025-10-30T10:35:36Z

src/FEM/GridPathSegmenter.hpp

+  const ippl::Vector<T,Dim>& h)
+{
+  CutTimes<Dim,T> cuts;
+  for (unsigned a=0; a<Dim; ++a) cuts.t[a] = (T)2; // sentinel (>1)


maybe make this comment more explicit (e.g. "initialize to 2 since value >1 means no crossing")

s-mayani · 2025-10-30T10:43:50Z

src/FEM/GridPathSegmenter.hpp

+template<unsigned Dim, typename T>
+KOKKOS_INLINE_FUNCTION
+void make_endpoints_fixed(const CutTimes<Dim,T>& cuts,
+                          std::array<T,Dim+2>& Tcuts /*out*/)


same comment as before for std::array

s-mayani · 2025-10-30T12:11:55Z

src/FEM/GridPathSegmenter.hpp

+  }
+
+  const T one = (T)1;
+  auto clamp_or_one = [&](T v)->T { return (v >= (T)1.5) ? one : (v < one ? v : one); };


Isn't this effectively the same as just doing (v >= 1) ? one: v?

s-mayani

Have reviewed the grid path segmenter, and the preliminary structure of the current projecter looks good.

JoMee added 4 commits October 17, 2025 09:16

Implement function to split up particle through multiple cells

8ff4da7

Add test skeleton for the current interpolation function as well as d…

a87c646

…efine that function

Change type of value to interpolate to ippl::Vector

7771196

wrap up

bdcc36d

aaadelmann requested review from aaadelmann, mohsensadr, s-mayani and srikrrish October 29, 2025 00:55

aaadelmann added the gitlab-mirror label Oct 29, 2025

mohsensadr requested changes Oct 29, 2025

View reviewed changes

s-mayani reviewed Oct 30, 2025

View reviewed changes

		f(std::integral_constant<int, I>{});
		static_for<I+1, End>(std::forward<F>(f));

Interpolate Current #437

Are you sure you want to change the base?

Interpolate Current #437

Uh oh!

Conversation

JoMee commented Oct 28, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mohsensadr Oct 30, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

s-mayani Oct 30, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

s-mayani left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

mohsensadr Oct 30, 2025 •

edited

Loading

s-mayani Oct 30, 2025 •

edited

Loading