Mirror pyccel modules as CUDA namespaces

CUDA_STRATEGY.md

@@ -104,6 +104,7 @@ Conventions:
 - The folder name, the pyccel function name, the CUDA `extern "C" __global__` function name and the name of the declared `Kernel` are all the same.
 - The pyccel file is `<name>_kernels.py` (so `struphy compile` picks it up); the CUDA file is `<name>_cuda.cu`.
 - Shared CUDA code (device helper functions, argument structs) lives in `.cuh` headers next to the pyccel module it mirrors, e.g. `bsplines/bsplines_kernels.cuh` for `bsplines/bsplines_kernels.py`. A `.cu` file includes them relative to the source root, e.g. `#include "struphy/kernel_arguments/pusher_args.cuh"`.
+- The device functions of a header live in the namespace `struphy_cuda::<module>`, where `<module>` is the name of the pyccel module they port: `struphy_cuda::bsplines_kernels` (`bsplines/bsplines_kernels.cuh`), `struphy_cuda::evaluation_kernels` (`geometry/evaluation_kernels.cuh`), `struphy_cuda::cuboid_kernels` (`geometry/domains/cuboid/cuboid_cuda.cuh`, which ports `cuboid_kernels.py`), and so on. CUDA-only helpers and constants go into the namespace of the module whose calculation they extract (`bsplines_kernels::MAX_SPLINE_DEGREE`, `evaluation_kernels_3d::SplineScratch`). Calls into another module are qualified with the module name, as the pyccel code calls `linalg_kernels.matrix_inv`: `linalg_kernels::matrix_inv(...)`; calls within the same module are not. A `.cu` file writes `using namespace struphy_cuda;` after its includes and qualifies every helper call (`evaluation_kernels_3d::get_spans(...)`). Only `struphy_cuda::pi` (`geometry/domains/constants_cuda.cuh`) and the generated argument structs (`MarkerArgs`, `DomainArgs`, `DerhamArgs`, `LocalProjectorsArgs`, which are global because cunumpy matches the kernel signature against their names) are outside a module namespace.
 - Every function in a `.cu` or `.cuh` file has a documentation comment (`/** ... */`, the C++ equivalent of a docstring) immediately above its definition. Identify the pyccel counterpart and explain the operation, parameters, outputs or return value, array shapes/layouts, and any CUDA-specific limits or storage requirements. Document CUDA-only helpers by identifying the pyccel calculation they extract.
 - CUDA parameter, local variable and scratch-field names exactly match the corresponding pyccel names, including spelling and capitalization (e.g. `eta1`, `args_markers`, `args_domain`, `args_derham`, `df_out`, `dfm`, `dfinv`, `v_logical`, `tn`, `pn`, `pd`, `il`, `det_a`, `span1`, `bn1`, `bd1`). Do not abbreviate them to `x`, `m`, `d`, `a`, or `out` when the pyccel counterpart uses a different name. Additional CUDA-only variables or arguments, such as pointer lengths, thread indices or per-thread scratch storage, must have descriptive names and documented purposes.
 - Kernels read and write markers through the `Array2D<double>` view: `args.markers(ip, j)`. The view carries the existing device pointer, shape and element strides; no array data is copied. The column count is `args.markers.shape[1]`.

src/struphy/bsplines/bsplines_kernels.cuh

@@ -1,5 +1,5 @@
 #pragma once
-namespace struphy_cuda {
+namespace struphy_cuda::bsplines_kernels {
 // Fixed capacity of the per-thread spline scratch arrays.
 constexpr int MAX_SPLINE_DEGREE = 8;
 
@@ -168,4 +168,4 @@ __device__ inline void b_der_splines_slim(const double* tn, int pn, double eta,
     for (int j = 0; j <= pn; ++j) bn[j] = values[pn][j];
     der[pn] = saved;
 }
-}
+}  // namespace struphy_cuda::bsplines_kernels

src/struphy/bsplines/evaluation_kernels_2d.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "cunumpy/array_view.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::evaluation_kernels_2d {
 /**
  * Sum the non-zero contributions of a 2d spline, as in evaluation_kernels_2d.evaluation_kernel_2d.
  *
@@ -26,4 +26,4 @@ __device__ inline double evaluation_kernel_2d(int p1, int p2, const double* basi
     }
     return spline_value;
 }
-}
+}  // namespace struphy_cuda::evaluation_kernels_2d

src/struphy/bsplines/evaluation_kernels_3d.cuh

@@ -2,12 +2,15 @@
 #include "struphy/kernel_arguments/pusher_args.cuh"
 #include "struphy/kernel_arguments/spline_args.cuh"
 #include "struphy/bsplines/bsplines_kernels.cuh"
-namespace struphy_cuda {
+#include "cunumpy/array_view.cuh"
+namespace struphy_cuda::evaluation_kernels_3d {
 /** Per-thread outputs of get_spans; spline names match DerhamArguments. */
 struct SplineScratch {
     int span1, span2, span3;
-    double bn1[MAX_SPLINE_DEGREE + 1], bn2[MAX_SPLINE_DEGREE + 1], bn3[MAX_SPLINE_DEGREE + 1];
-    double bd1[MAX_SPLINE_DEGREE], bd2[MAX_SPLINE_DEGREE], bd3[MAX_SPLINE_DEGREE];
+    double bn1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], bn2[bsplines_kernels::MAX_SPLINE_DEGREE + 1],
+        bn3[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
+    double bd1[bsplines_kernels::MAX_SPLINE_DEGREE], bd2[bsplines_kernels::MAX_SPLINE_DEGREE],
+        bd3[bsplines_kernels::MAX_SPLINE_DEGREE];
 };
 
 /**
@@ -25,20 +28,17 @@ struct SplineScratch {
  */
 __device__ inline void get_spans(double eta1, double eta2, double eta3,
                                 const DerhamArgs& args_derham, SplineScratch& scratch) {
-    int span1 = find_span(args_derham.tn1, args_derham.nt1, args_derham.pn[0], eta1);
-    int span2 = find_span(args_derham.tn2, args_derham.nt2, args_derham.pn[1], eta2);
-    int span3 = find_span(args_derham.tn3, args_derham.nt3, args_derham.pn[2], eta3);
+    int span1 = bsplines_kernels::find_span(args_derham.tn1, args_derham.nt1, args_derham.pn[0], eta1);
+    int span2 = bsplines_kernels::find_span(args_derham.tn2, args_derham.nt2, args_derham.pn[1], eta2);
+    int span3 = bsplines_kernels::find_span(args_derham.tn3, args_derham.nt3, args_derham.pn[2], eta3);
     scratch.span1 = span1;
     scratch.span2 = span2;
     scratch.span3 = span3;
-    b_d_splines_slim(args_derham.tn1, args_derham.pn[0], eta1, span1, scratch.bn1, scratch.bd1);
-    b_d_splines_slim(args_derham.tn2, args_derham.pn[1], eta2, span2, scratch.bn2, scratch.bd2);
-    b_d_splines_slim(args_derham.tn3, args_derham.pn[2], eta3, span3, scratch.bn3, scratch.bd3);
-}
+    bsplines_kernels::b_d_splines_slim(args_derham.tn1, args_derham.pn[0], eta1, span1, scratch.bn1, scratch.bd1);
+    bsplines_kernels::b_d_splines_slim(args_derham.tn2, args_derham.pn[1], eta2, span2, scratch.bn2, scratch.bd2);
+    bsplines_kernels::b_d_splines_slim(args_derham.tn3, args_derham.pn[2], eta3, span3, scratch.bn3, scratch.bd3);
 }
 
-#include "cunumpy/array_view.cuh"
-namespace struphy_cuda {
 /**
  * Sum the non-zero contributions of a distributed spline, as in evaluation_kernels_3d.eval_spline_mpi_kernel.
  *
@@ -139,7 +139,6 @@ __device__ inline void eval_2form_spline_mpi(int span1, int span2, int span3, co
     out[2] = eval_spline_mpi_kernel(args_derham.pn[0] - 1, args_derham.pn[1] - 1, args_derham.pn[2], scratch.bd1,
                                     scratch.bd2, scratch.bn3, span1, span2, span3, form_coeffs_3, args_derham.starts);
 }
-}
 
 // Point-wise evaluation shared by the eval_spline_mpi_* kernels (bsplines/kernels/<name>/<name>_cuda.cu).
 /**
@@ -158,31 +157,33 @@ __device__ inline void eval_2form_spline_mpi(int span1, int span2, int span3, co
  */
 __device__ inline double eval_spline_mpi(double eta1, double eta2, double eta3, Array3D<double> _data,
                                          const SplineArgs& args_spline) {
-    struphy_cuda::SplineScratch scratch;
+    SplineScratch scratch;
 
     // get spline values at eta
-    scratch.span1 = struphy_cuda::find_span(args_spline.tn1.data, args_spline.tn1.shape[0], args_spline.pn[0], eta1);
-    scratch.span2 = struphy_cuda::find_span(args_spline.tn2.data, args_spline.tn2.shape[0], args_spline.pn[1], eta2);
-    scratch.span3 = struphy_cuda::find_span(args_spline.tn3.data, args_spline.tn3.shape[0], args_spline.pn[2], eta3);
-    struphy_cuda::b_d_splines_slim(args_spline.tn1.data, args_spline.pn[0], eta1, scratch.span1, scratch.bn1,
-                                   scratch.bd1);
-    struphy_cuda::b_d_splines_slim(args_spline.tn2.data, args_spline.pn[1], eta2, scratch.span2, scratch.bn2,
-                                   scratch.bd2);
-    struphy_cuda::b_d_splines_slim(args_spline.tn3.data, args_spline.pn[2], eta3, scratch.span3, scratch.bn3,
-                                   scratch.bd3);
+    scratch.span1 =
+        bsplines_kernels::find_span(args_spline.tn1.data, args_spline.tn1.shape[0], args_spline.pn[0], eta1);
+    scratch.span2 =
+        bsplines_kernels::find_span(args_spline.tn2.data, args_spline.tn2.shape[0], args_spline.pn[1], eta2);
+    scratch.span3 =
+        bsplines_kernels::find_span(args_spline.tn3.data, args_spline.tn3.shape[0], args_spline.pn[2], eta3);
+    bsplines_kernels::b_d_splines_slim(args_spline.tn1.data, args_spline.pn[0], eta1, scratch.span1, scratch.bn1,
+                                       scratch.bd1);
+    bsplines_kernels::b_d_splines_slim(args_spline.tn2.data, args_spline.pn[1], eta2, scratch.span2, scratch.bn2,
+                                       scratch.bd2);
+    bsplines_kernels::b_d_splines_slim(args_spline.tn3.data, args_spline.pn[2], eta3, scratch.span3, scratch.bn3,
+                                       scratch.bd3);
 
     const double* b1 = args_spline.kind[0] == 0 ? scratch.bn1 : scratch.bd1;
     const double* b2 = args_spline.kind[1] == 0 ? scratch.bn2 : scratch.bd2;
     const double* b3 = args_spline.kind[2] == 0 ? scratch.bn3 : scratch.bd3;
 
-    double value = struphy_cuda::eval_spline_mpi_kernel(
+    double value = eval_spline_mpi_kernel(
         args_spline.pn[0] - args_spline.kind[0], args_spline.pn[1] - args_spline.kind[1],
         args_spline.pn[2] - args_spline.kind[2], b1, b2, b3, scratch.span1, scratch.span2, scratch.span3, _data,
         args_spline.starts);
     return value;
 }
 
-namespace struphy_cuda {
 /**
  * Sum the non-zero contributions of a spline with global coefficients, as in evaluation_kernels_3d.evaluation_kernel_3d.
  *
@@ -214,4 +215,4 @@ __device__ inline double evaluation_kernel_3d(int p1, int p2, int p3, const doub
     }
     return spline_value;
 }
-}
+}  // namespace struphy_cuda::evaluation_kernels_3d

src/struphy/bsplines/kernels/eval_spline_mpi_markers/eval_spline_mpi_markers_cuda.cu

@@ -1,6 +1,8 @@
 #include "struphy/bsplines/evaluation_kernels_3d.cuh"
 #include "struphy/kernel_arguments/spline_args.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in eval_spline_mpi_markers_kernels.py. Array views carry
 // dimensions and strides, so no CUDA-only lengths are needed.
 
@@ -26,5 +28,6 @@ extern "C" __global__ void eval_spline_mpi_markers(Array2D<double> markers, Arra
     // point not in process domain
     if (markers(ip, 0) == -1.) return;
 
-    values(ip) = eval_spline_mpi(markers(ip, 0), markers(ip, 1), markers(ip, 2), _data, args_spline);
+    values(ip) = evaluation_kernels_3d::eval_spline_mpi(markers(ip, 0), markers(ip, 1), markers(ip, 2), _data,
+                                                        args_spline);
 }

src/struphy/bsplines/kernels/eval_spline_mpi_matrix/eval_spline_mpi_matrix_cuda.cu

@@ -1,6 +1,8 @@
 #include "struphy/bsplines/evaluation_kernels_3d.cuh"
 #include "struphy/kernel_arguments/spline_args.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in eval_spline_mpi_matrix_kernels.py. Array views carry
 // dimensions and strides, so no CUDA-only lengths are needed.
 
@@ -30,5 +32,6 @@ extern "C" __global__ void eval_spline_mpi_matrix(Array3D<double> eta1, Array3D<
     // point not in process domain
     if (eta1(i, j, k) == -1. || eta2(i, j, k) == -1. || eta3(i, j, k) == -1.) return;
 
-    values(i, j, k) = eval_spline_mpi(eta1(i, j, k), eta2(i, j, k), eta3(i, j, k), _data, args_spline);
+    values(i, j, k) = evaluation_kernels_3d::eval_spline_mpi(eta1(i, j, k), eta2(i, j, k), eta3(i, j, k), _data,
+                                                             args_spline);
 }

src/struphy/bsplines/kernels/eval_spline_mpi_sparse_meshgrid/eval_spline_mpi_sparse_meshgrid_cuda.cu

@@ -1,6 +1,8 @@
 #include "struphy/bsplines/evaluation_kernels_3d.cuh"
 #include "struphy/kernel_arguments/spline_args.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in eval_spline_mpi_sparse_meshgrid_kernels.py. Array views
 // carry dimensions and strides, so no CUDA-only lengths are needed.
 
@@ -30,5 +32,6 @@ extern "C" __global__ void eval_spline_mpi_sparse_meshgrid(Array3D<double> eta1,
     // point not in process domain
     if (eta1(i, 0, 0) == -1. || eta2(0, j, 0) == -1. || eta3(0, 0, k) == -1.) return;
 
-    values(i, j, k) = eval_spline_mpi(eta1(i, 0, 0), eta2(0, j, 0), eta3(0, 0, k), _data, args_spline);
+    values(i, j, k) = evaluation_kernels_3d::eval_spline_mpi(eta1(i, 0, 0), eta2(0, j, 0), eta3(0, 0, k), _data,
+                                                             args_spline);
 }

src/struphy/bsplines/kernels/eval_spline_mpi_tensor_product_fixed/eval_spline_mpi_tensor_product_fixed_cuda.cu

@@ -1,5 +1,7 @@
 #include "struphy/bsplines/evaluation_kernels_3d.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in eval_spline_mpi_tensor_product_fixed_kernels.py. Array
 // views carry dimensions and strides, so no CUDA-only lengths are needed.
 
@@ -42,12 +44,13 @@ extern "C" __global__ void eval_spline_mpi_tensor_product_fixed(Array1D<long lon
     long long j = (ijk / nk) % nj;
     long long i = ijk / (nj * nk);
 
-    double b1[struphy_cuda::MAX_SPLINE_DEGREE + 1], b2[struphy_cuda::MAX_SPLINE_DEGREE + 1],
-        b3[struphy_cuda::MAX_SPLINE_DEGREE + 1];
+    double b1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], b2[bsplines_kernels::MAX_SPLINE_DEGREE + 1],
+        b3[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
     for (int il = 0; il <= pn[0] - kind[0]; ++il) b1[il] = b1s(i, il);
     for (int il = 0; il <= pn[1] - kind[1]; ++il) b2[il] = b2s(j, il);
     for (int il = 0; il <= pn[2] - kind[2]; ++il) b3[il] = b3s(k, il);
 
-    values(i, j, k) = struphy_cuda::eval_spline_mpi_kernel(pn[0] - kind[0], pn[1] - kind[1], pn[2] - kind[2], b1, b2,
-                                                           b3, span1s(i), span2s(j), span3s(k), _data, starts);
+    values(i, j, k) =
+        evaluation_kernels_3d::eval_spline_mpi_kernel(pn[0] - kind[0], pn[1] - kind[1], pn[2] - kind[2], b1, b2, b3,
+                                                      span1s(i), span2s(j), span3s(k), _data, starts);
 }

src/struphy/geometry/domains/colella/colella_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::colella_kernels {
 /**
  * Evaluate the Colella mapping, as in colella_kernels.colella.
  *
@@ -46,4 +46,4 @@ __device__ inline void colella_df(double eta1, double eta2, double lx, double ly
     df_out[7] = 0.0;
     df_out[8] = lz;
 }
-}
+}  // namespace struphy_cuda::colella_kernels

src/struphy/geometry/domains/cuboid/cuboid_cuda.cuh

@@ -1,5 +1,5 @@
 #pragma once
-namespace struphy_cuda {
+namespace struphy_cuda::cuboid_kernels {
 /**
  * Evaluate the cuboid mapping, as in cuboid_kernels.cuboid.
  *
@@ -45,4 +45,4 @@ __device__ inline void cuboid_df(double l1, double r1, double l2, double r2, dou
     df_out[7] = 0.;
     df_out[8] = r3 - l3;
 }
-}
+}  // namespace struphy_cuda::cuboid_kernels

src/struphy/geometry/domains/hollow_cylinder/hollow_cylinder_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::hollow_cylinder_kernels {
 /**
  * Evaluate the hollow cylinder mapping, as in hollow_cylinder_kernels.hollow_cyl.
  *
@@ -48,4 +48,4 @@ __device__ inline void hollow_cyl_df(double eta1, double eta2, double a1, double
     df_out[7] = 0.0;
     df_out[8] = lz;
 }
-}
+}  // namespace struphy_cuda::hollow_cylinder_kernels

src/struphy/geometry/domains/hollow_torus/hollow_torus_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::hollow_torus_kernels {
 /**
  * Evaluate the hollow torus mapping, as in hollow_torus_kernels.hollow_torus.
  *
@@ -98,4 +98,4 @@ __device__ inline void hollow_torus_df(double eta1, double eta2, double eta3, do
         df_out[8] = 0.0;
     }
 }
-}
+}  // namespace struphy_cuda::hollow_torus_kernels

src/struphy/geometry/domains/orthogonal/orthogonal_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::orthogonal_kernels {
 /**
  * Evaluate the orthogonal mapping, as in orthogonal_kernels.orthogonal.
  *
@@ -45,4 +45,4 @@ __device__ inline void orthogonal_df(double eta1, double eta2, double lx, double
     df_out[7] = 0.0;
     df_out[8] = lz;
 }
-}
+}  // namespace struphy_cuda::orthogonal_kernels

src/struphy/geometry/domains/powered_elliptic_cylinder/powered_elliptic_cylinder_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::powered_elliptic_cylinder_kernels {
 /**
  * Evaluate the powered elliptic cylinder mapping, as in powered_elliptic_cylinder_kernels.powered_ellipse.
  *
@@ -47,4 +47,4 @@ __device__ inline void powered_ellipse_df(double eta1, double eta2, double eta3,
     df_out[7] = 0.0;
     df_out[8] = lz;
 }
-}
+}  // namespace struphy_cuda::powered_elliptic_cylinder_kernels

src/struphy/geometry/domains/shafranov_dshaped_cylinder/shafranov_dshaped_cylinder_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::shafranov_dshaped_cylinder_kernels {
 /**
  * Evaluate the D-shaped Shafranov mapping, as in shafranov_dshaped_cylinder_kernels.shafranov_dshaped.
  *
@@ -59,4 +59,4 @@ __device__ inline void shafranov_dshaped_df(double eta1, double eta2, double eta
     df_out[7] = 0.0;
     df_out[8] = lz;
 }
-}  // namespace struphy_cuda
+}  // namespace struphy_cuda::shafranov_dshaped_cylinder_kernels

src/struphy/geometry/domains/shafranov_shift_cylinder/shafranov_shift_cylinder_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::shafranov_shift_cylinder_kernels {
 /**
  * Evaluate the Shafranov shift mapping, as in shafranov_shift_cylinder_kernels.shafranov_shift.
  *
@@ -46,4 +46,4 @@ __device__ inline void shafranov_shift_df(double eta1, double eta2, double eta3,
     df_out[7] = 0.0;
     df_out[8] = lz;
 }
-}
+}  // namespace struphy_cuda::shafranov_shift_cylinder_kernels

src/struphy/geometry/domains/shafranov_sqrt_cylinder/shafranov_sqrt_cylinder_cuda.cuh

@@ -1,6 +1,6 @@
 #pragma once
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::shafranov_sqrt_cylinder_kernels {
 /**
  * Evaluate the Shafranov sqrt mapping, as in shafranov_sqrt_cylinder_kernels.shafranov_sqrt.
  *
@@ -46,4 +46,4 @@ __device__ inline void shafranov_sqrt_df(double eta1, double eta2, double eta3,
     df_out[7] = 0.0;
     df_out[8] = lz;
 }
-}
+}  // namespace struphy_cuda::shafranov_sqrt_cylinder_kernels