diff --git a/CUDA_STRATEGY.md b/CUDA_STRATEGY.md
index a3b5c4671..a7724a594 100644
--- a/CUDA_STRATEGY.md
+++ b/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]`.
diff --git a/src/struphy/bsplines/bsplines_kernels.cuh b/src/struphy/bsplines/bsplines_kernels.cuh
index 7527b0731..46f199cc3 100644
--- a/src/struphy/bsplines/bsplines_kernels.cuh
+++ b/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
diff --git a/src/struphy/bsplines/evaluation_kernels_2d.cuh b/src/struphy/bsplines/evaluation_kernels_2d.cuh
index f53a384bf..bb4e5a0f2 100644
--- a/src/struphy/bsplines/evaluation_kernels_2d.cuh
+++ b/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
diff --git a/src/struphy/bsplines/evaluation_kernels_3d.cuh b/src/struphy/bsplines/evaluation_kernels_3d.cuh
index ede3e5f40..02c704a00 100644
--- a/src/struphy/bsplines/evaluation_kernels_3d.cuh
+++ b/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
diff --git a/src/struphy/bsplines/kernels/eval_spline_mpi_markers/eval_spline_mpi_markers_cuda.cu b/src/struphy/bsplines/kernels/eval_spline_mpi_markers/eval_spline_mpi_markers_cuda.cu
index c39ccfdda..08314d224 100644
--- a/src/struphy/bsplines/kernels/eval_spline_mpi_markers/eval_spline_mpi_markers_cuda.cu
+++ b/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);
 }
diff --git a/src/struphy/bsplines/kernels/eval_spline_mpi_matrix/eval_spline_mpi_matrix_cuda.cu b/src/struphy/bsplines/kernels/eval_spline_mpi_matrix/eval_spline_mpi_matrix_cuda.cu
index 61cceedc0..553e8ecea 100644
--- a/src/struphy/bsplines/kernels/eval_spline_mpi_matrix/eval_spline_mpi_matrix_cuda.cu
+++ b/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);
 }
diff --git a/src/struphy/bsplines/kernels/eval_spline_mpi_sparse_meshgrid/eval_spline_mpi_sparse_meshgrid_cuda.cu b/src/struphy/bsplines/kernels/eval_spline_mpi_sparse_meshgrid/eval_spline_mpi_sparse_meshgrid_cuda.cu
index 1fde8a9a3..322da4fa8 100644
--- a/src/struphy/bsplines/kernels/eval_spline_mpi_sparse_meshgrid/eval_spline_mpi_sparse_meshgrid_cuda.cu
+++ b/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);
 }
diff --git a/src/struphy/bsplines/kernels/eval_spline_mpi_tensor_product_fixed/eval_spline_mpi_tensor_product_fixed_cuda.cu b/src/struphy/bsplines/kernels/eval_spline_mpi_tensor_product_fixed/eval_spline_mpi_tensor_product_fixed_cuda.cu
index fa7716720..36736ca29 100644
--- a/src/struphy/bsplines/kernels/eval_spline_mpi_tensor_product_fixed/eval_spline_mpi_tensor_product_fixed_cuda.cu
+++ b/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);
 }
diff --git a/src/struphy/geometry/domains/colella/colella_cuda.cuh b/src/struphy/geometry/domains/colella/colella_cuda.cuh
index 74211e514..82c3c2136 100644
--- a/src/struphy/geometry/domains/colella/colella_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/cuboid/cuboid_cuda.cuh b/src/struphy/geometry/domains/cuboid/cuboid_cuda.cuh
index 352d10ee9..e830524dd 100644
--- a/src/struphy/geometry/domains/cuboid/cuboid_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/hollow_cylinder/hollow_cylinder_cuda.cuh b/src/struphy/geometry/domains/hollow_cylinder/hollow_cylinder_cuda.cuh
index bf1c7ce5a..9b8724b61 100644
--- a/src/struphy/geometry/domains/hollow_cylinder/hollow_cylinder_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/hollow_torus/hollow_torus_cuda.cuh b/src/struphy/geometry/domains/hollow_torus/hollow_torus_cuda.cuh
index a3db4b672..53896cdfc 100644
--- a/src/struphy/geometry/domains/hollow_torus/hollow_torus_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/orthogonal/orthogonal_cuda.cuh b/src/struphy/geometry/domains/orthogonal/orthogonal_cuda.cuh
index 769284fdb..03e41b623 100644
--- a/src/struphy/geometry/domains/orthogonal/orthogonal_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/powered_elliptic_cylinder/powered_elliptic_cylinder_cuda.cuh b/src/struphy/geometry/domains/powered_elliptic_cylinder/powered_elliptic_cylinder_cuda.cuh
index 4fe528025..c7812a7e5 100644
--- a/src/struphy/geometry/domains/powered_elliptic_cylinder/powered_elliptic_cylinder_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/shafranov_dshaped_cylinder/shafranov_dshaped_cylinder_cuda.cuh b/src/struphy/geometry/domains/shafranov_dshaped_cylinder/shafranov_dshaped_cylinder_cuda.cuh
index 5d2b993f3..b11ec974c 100644
--- a/src/struphy/geometry/domains/shafranov_dshaped_cylinder/shafranov_dshaped_cylinder_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/shafranov_shift_cylinder/shafranov_shift_cylinder_cuda.cuh b/src/struphy/geometry/domains/shafranov_shift_cylinder/shafranov_shift_cylinder_cuda.cuh
index a2aa07b7b..a108d6ee7 100644
--- a/src/struphy/geometry/domains/shafranov_shift_cylinder/shafranov_shift_cylinder_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/domains/shafranov_sqrt_cylinder/shafranov_sqrt_cylinder_cuda.cuh b/src/struphy/geometry/domains/shafranov_sqrt_cylinder/shafranov_sqrt_cylinder_cuda.cuh
index beebcfef8..45aa3cb0d 100644
--- a/src/struphy/geometry/domains/shafranov_sqrt_cylinder/shafranov_sqrt_cylinder_cuda.cuh
+++ b/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
diff --git a/src/struphy/geometry/evaluation_kernels.cuh b/src/struphy/geometry/evaluation_kernels.cuh
index d27ed0ac6..eafebf6b9 100644
--- a/src/struphy/geometry/evaluation_kernels.cuh
+++ b/src/struphy/geometry/evaluation_kernels.cuh
@@ -16,7 +16,7 @@
 #include "struphy/geometry/domains/shafranov_shift_cylinder/shafranov_shift_cylinder_cuda.cuh"
 #include "struphy/geometry/domains/shafranov_sqrt_cylinder/shafranov_sqrt_cylinder_cuda.cuh"
 #include "struphy/geometry/domains/shafranov_dshaped_cylinder/shafranov_dshaped_cylinder_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::evaluation_kernels {
 /**
  * Point-wise evaluation of (x, y, z) = F(eta1, eta2, eta3), as in geometry.evaluation_kernels.f.
  *
@@ -32,43 +32,53 @@ namespace struphy_cuda {
 __device__ inline void f(double eta1, double eta2, double eta3, const DomainArgs& args, double* f_out) {
     switch (args.kind_map) {
         case 0:
-            spline_3d(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args.ind3, args, f_out);
+            spline_mappings_kernels::spline_3d(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args.ind3, args,
+                                               f_out);
             return;
         case 1:
-            spline_2d_straight(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args, args.params[0], f_out);
+            spline_mappings_kernels::spline_2d_straight(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args,
+                                                        args.params[0], f_out);
             return;
         case 2:
-            spline_2d_torus(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args, args.params[0], f_out);
+            spline_mappings_kernels::spline_2d_torus(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args,
+                                                     args.params[0], f_out);
             return;
         case 10:
-            cuboid(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], args.params[4],
-                   args.params[5], f_out);
+            cuboid_kernels::cuboid(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
+                                   args.params[4], args.params[5], f_out);
             return;
         case 11:
-            orthogonal(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], f_out);
+            orthogonal_kernels::orthogonal(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2],
+                                           args.params[3], f_out);
             return;
         case 12:
-            colella(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], f_out);
+            colella_kernels::colella(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
+                                     f_out);
             return;
         case 20:
-            hollow_cyl(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], f_out);
+            hollow_cylinder_kernels::hollow_cyl(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2],
+                                                args.params[3], f_out);
             return;
         case 21:
-            powered_ellipse(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], f_out);
+            powered_elliptic_cylinder_kernels::powered_ellipse(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                               args.params[2], args.params[3], f_out);
             return;
         case 22:
-            hollow_torus(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
-                         args.params[4], args.params[5], f_out);
+            hollow_torus_kernels::hollow_torus(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2],
+                                               args.params[3], args.params[4], args.params[5], f_out);
             return;
         case 30:
-            shafranov_shift(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], f_out);
+            shafranov_shift_cylinder_kernels::shafranov_shift(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                              args.params[2], args.params[3], f_out);
             return;
         case 31:
-            shafranov_sqrt(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], f_out);
+            shafranov_sqrt_cylinder_kernels::shafranov_sqrt(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                            args.params[2], args.params[3], f_out);
             return;
         case 32:
-            shafranov_dshaped(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
-                              args.params[4], args.params[5], args.params[6], f_out);
+            shafranov_dshaped_cylinder_kernels::shafranov_dshaped(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                                  args.params[2], args.params[3], args.params[4],
+                                                                  args.params[5], args.params[6], f_out);
             return;
         default:
             asm("trap;");  // unknown kind_map
@@ -90,45 +100,53 @@ __device__ inline void f(double eta1, double eta2, double eta3, const DomainArgs
 __device__ inline void df(double eta1, double eta2, double eta3, const DomainArgs& args, double* df_out) {
     switch (args.kind_map) {
         case 0:
-            spline_3d_df(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args.ind3, args, df_out);
+            spline_mappings_kernels::spline_3d_df(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args.ind3, args,
+                                                  df_out);
             return;
         case 1:
-            spline_2d_straight_df(eta1, eta2, args.degree, args.ind1, args.ind2, args, args.params[0], df_out);
+            spline_mappings_kernels::spline_2d_straight_df(eta1, eta2, args.degree, args.ind1, args.ind2, args,
+                                                           args.params[0], df_out);
             return;
         case 2:
-            spline_2d_torus_df(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args, args.params[0], df_out);
+            spline_mappings_kernels::spline_2d_torus_df(eta1, eta2, eta3, args.degree, args.ind1, args.ind2, args,
+                                                        args.params[0], df_out);
             return;
         case 10:
-            cuboid_df(args.params[0], args.params[1], args.params[2], args.params[3], args.params[4], args.params[5],
-                      df_out);
+            cuboid_kernels::cuboid_df(args.params[0], args.params[1], args.params[2], args.params[3], args.params[4],
+                                      args.params[5], df_out);
             return;
         case 11:
-            orthogonal_df(eta1, eta2, args.params[0], args.params[1], args.params[2], args.params[3], df_out);
+            orthogonal_kernels::orthogonal_df(eta1, eta2, args.params[0], args.params[1], args.params[2],
+                                              args.params[3], df_out);
             return;
         case 12:
-            colella_df(eta1, eta2, args.params[0], args.params[1], args.params[2], args.params[3], df_out);
+            colella_kernels::colella_df(eta1, eta2, args.params[0], args.params[1], args.params[2], args.params[3],
+                                        df_out);
             return;
         case 20:
-            hollow_cyl_df(eta1, eta2, args.params[0], args.params[1], args.params[2], args.params[3], df_out);
+            hollow_cylinder_kernels::hollow_cyl_df(eta1, eta2, args.params[0], args.params[1], args.params[2],
+                                                   args.params[3], df_out);
             return;
         case 21:
-            powered_ellipse_df(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
-                               df_out);
+            powered_elliptic_cylinder_kernels::powered_ellipse_df(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                                  args.params[2], args.params[3], df_out);
             return;
         case 22:
-            hollow_torus_df(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
-                            args.params[4], args.params[5], df_out);
+            hollow_torus_kernels::hollow_torus_df(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2],
+                                                  args.params[3], args.params[4], args.params[5], df_out);
             return;
         case 30:
-            shafranov_shift_df(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
-                               df_out);
+            shafranov_shift_cylinder_kernels::shafranov_shift_df(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                                 args.params[2], args.params[3], df_out);
             return;
         case 31:
-            shafranov_sqrt_df(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3], df_out);
+            shafranov_sqrt_cylinder_kernels::shafranov_sqrt_df(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                               args.params[2], args.params[3], df_out);
             return;
         case 32:
-            shafranov_dshaped_df(eta1, eta2, eta3, args.params[0], args.params[1], args.params[2], args.params[3],
-                                 args.params[4], args.params[5], args.params[6], df_out);
+            shafranov_dshaped_cylinder_kernels::shafranov_dshaped_df(eta1, eta2, eta3, args.params[0], args.params[1],
+                                                                     args.params[2], args.params[3], args.params[4],
+                                                                     args.params[5], args.params[6], df_out);
             return;
         default:
             asm("trap;");  // unknown kind_map
@@ -147,7 +165,7 @@ __device__ inline void df(double eta1, double eta2, double eta3, const DomainArg
  */
 __device__ inline double det_df(double eta1, double eta2, double eta3, const DomainArgs& args, double* tmp1) {
     df(eta1, eta2, eta3, args, tmp1);
-    double detdf = det(tmp1);
+    double detdf = linalg_kernels::det(tmp1);
     return detdf;
 }
 
@@ -165,7 +183,7 @@ __device__ inline double det_df(double eta1, double eta2, double eta3, const Dom
 __device__ inline void df_inv(double eta1, double eta2, double eta3, const DomainArgs& args, double* tmp1,
                               bool avoid_round_off, double* dfinv_out) {
     df(eta1, eta2, eta3, args, tmp1);
-    matrix_inv(tmp1, dfinv_out);
+    linalg_kernels::matrix_inv(tmp1, dfinv_out);
 
     // set known (analytical) zero components manually to zero to avoid round-off error remainders!
     if (avoid_round_off) {
@@ -253,8 +271,8 @@ __device__ inline void df_inv(double eta1, double eta2, double eta3, const Domai
 __device__ inline void g(double eta1, double eta2, double eta3, const DomainArgs& args, double* tmp1, double* tmp2,
                          bool avoid_round_off, double* g_out) {
     df(eta1, eta2, eta3, args, tmp1);
-    transpose(tmp1, tmp2);
-    matrix_matrix(tmp2, tmp1, g_out);
+    linalg_kernels::transpose(tmp1, tmp2);
+    linalg_kernels::matrix_matrix(tmp2, tmp1, g_out);
 
     // set known (analytical) zero components manually to zero to avoid round-off error remainders!
     if (avoid_round_off) {
@@ -363,7 +381,7 @@ __device__ inline void g(double eta1, double eta2, double eta3, const DomainArgs
 __device__ inline void g_inv(double eta1, double eta2, double eta3, const DomainArgs& args, double* tmp1, double* tmp2,
                              double* tmp3, bool avoid_round_off, double* ginv_out) {
     g(eta1, eta2, eta3, args, tmp1, tmp2, avoid_round_off, tmp3);
-    matrix_inv(tmp3, ginv_out);
+    linalg_kernels::matrix_inv(tmp3, ginv_out);
 
     // set known (analytical) zero components manually to zero to avoid round-off error remainders!
     if (avoid_round_off) {
@@ -509,4 +527,4 @@ __device__ inline void select_metric_coeff(double eta1, double eta2, double eta3
         g_inv(eta1, eta2, eta3, args, tmp1, tmp2, tmp3, avoid_round_off, out);
     }
 }
-}  // namespace struphy_cuda
+}  // namespace struphy_cuda::evaluation_kernels
diff --git a/src/struphy/geometry/kernels/kernel_evaluate/kernel_evaluate_cuda.cu b/src/struphy/geometry/kernels/kernel_evaluate/kernel_evaluate_cuda.cu
index 6aff717a4..e6b9515d7 100644
--- a/src/struphy/geometry/kernels/kernel_evaluate/kernel_evaluate_cuda.cu
+++ b/src/struphy/geometry/kernels/kernel_evaluate/kernel_evaluate_cuda.cu
@@ -1,5 +1,7 @@
 #include "struphy/geometry/evaluation_kernels.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in kernel_evaluate_kernels.py. mat_f has five
 // dimensions, more than cunumpy's array views (Array1D-Array4D), so it is passed as a C-contiguous pointer.
 
@@ -47,7 +49,7 @@ extern "C" __global__ void kernel_evaluate(Array3D<double> eta1, Array3D<double>
 
     for (int k = 0; k < 9; ++k) out[k] = mat_f_point[k];
 
-    struphy_cuda::select_metric_coeff(e1, e2, e3, kind_coeff, args, tmp0, tmp1, tmp2, tmp3, avoid_round_off, out);
+    evaluation_kernels::select_metric_coeff(e1, e2, e3, kind_coeff, args, tmp0, tmp1, tmp2, tmp3, avoid_round_off, out);
 
     for (int k = 0; k < 9; ++k) mat_f_point[k] = out[k];
 }
diff --git a/src/struphy/geometry/kernels/kernel_evaluate_pic/kernel_evaluate_pic_cuda.cu b/src/struphy/geometry/kernels/kernel_evaluate_pic/kernel_evaluate_pic_cuda.cu
index 523ad1f6c..9644a2fca 100644
--- a/src/struphy/geometry/kernels/kernel_evaluate_pic/kernel_evaluate_pic_cuda.cu
+++ b/src/struphy/geometry/kernels/kernel_evaluate_pic/kernel_evaluate_pic_cuda.cu
@@ -1,5 +1,7 @@
 #include "struphy/geometry/evaluation_kernels.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in kernel_evaluate_pic_kernels.py.
 
 /**
@@ -70,8 +72,8 @@ extern "C" __global__ void kernel_evaluate_pic(Array2D<double> markers, int kind
         for (int k = 0; k < 3; ++k)
             for (int l = 0; l < 3; ++l) out[3 * k + l] = mat_f(counter, k, l);
 
-        struphy_cuda::select_metric_coeff(e1, e2, e3, kind_coeff, args, tmp0, tmp1, tmp2, tmp3, avoid_round_off,
-                                          out);
+        evaluation_kernels::select_metric_coeff(e1, e2, e3, kind_coeff, args, tmp0, tmp1, tmp2, tmp3, avoid_round_off,
+                                                out);
 
         for (int k = 0; k < 3; ++k)
             for (int l = 0; l < 3; ++l) mat_f(counter, k, l) = out[3 * k + l];
diff --git a/src/struphy/geometry/kernels/kernel_pullpush/kernel_pullpush_cuda.cu b/src/struphy/geometry/kernels/kernel_pullpush/kernel_pullpush_cuda.cu
index 952f979da..607eb312b 100644
--- a/src/struphy/geometry/kernels/kernel_pullpush/kernel_pullpush_cuda.cu
+++ b/src/struphy/geometry/kernels/kernel_pullpush/kernel_pullpush_cuda.cu
@@ -1,5 +1,7 @@
 #include "struphy/geometry/transform_kernels.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in kernel_pullpush_kernels.py.
 
 /**
@@ -50,11 +52,11 @@ extern "C" __global__ void kernel_pullpush(Array4D<double> a, Array3D<double> et
     for (long long k = 0; k < n_tmp2; ++k) tmp2[k] = out(i1, i2, i3, k);
 
     if (kind_transform == 0)
-        struphy_cuda::pull(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
+        transform_kernels::pull(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
     else if (kind_transform == 1)
-        struphy_cuda::push(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
+        transform_kernels::push(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
     else
-        struphy_cuda::tran(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
+        transform_kernels::tran(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
 
     for (long long k = 0; k < n_tmp2; ++k) out(i1, i2, i3, k) = tmp2[k];
 }
diff --git a/src/struphy/geometry/kernels/kernel_pullpush_pic/kernel_pullpush_pic_cuda.cu b/src/struphy/geometry/kernels/kernel_pullpush_pic/kernel_pullpush_pic_cuda.cu
index b58239686..466ce645f 100644
--- a/src/struphy/geometry/kernels/kernel_pullpush_pic/kernel_pullpush_pic_cuda.cu
+++ b/src/struphy/geometry/kernels/kernel_pullpush_pic/kernel_pullpush_pic_cuda.cu
@@ -1,5 +1,7 @@
 #include "struphy/geometry/transform_kernels.cuh"
 
+using namespace struphy_cuda;
+
 // Same arguments, in the same order, as the pyccel kernel in kernel_pullpush_pic_kernels.py.
 
 /**
@@ -83,11 +85,11 @@ extern "C" __global__ void kernel_pullpush_pic(Array2D<double> a, Array2D<double
         for (long long k = 0; k < n_tmp2; ++k) tmp2[k] = out(counter_o, k);
 
         if (kind_transform == 0)
-            struphy_cuda::pull(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
+            transform_kernels::pull(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
         else if (kind_transform == 1)
-            struphy_cuda::push(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
+            transform_kernels::push(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
         else
-            struphy_cuda::tran(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
+            transform_kernels::tran(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2);
 
         for (long long k = 0; k < n_tmp2; ++k) out(counter_o, k) = tmp2[k];
     }
diff --git a/src/struphy/geometry/spline_mappings_kernels.cuh b/src/struphy/geometry/spline_mappings_kernels.cuh
index d37581dad..8290a517d 100644
--- a/src/struphy/geometry/spline_mappings_kernels.cuh
+++ b/src/struphy/geometry/spline_mappings_kernels.cuh
@@ -5,7 +5,7 @@
 #include "struphy/bsplines/evaluation_kernels_2d.cuh"
 #include "struphy/bsplines/evaluation_kernels_3d.cuh"
 #include "struphy/geometry/domains/constants_cuda.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::spline_mappings_kernels {
 /**
  * Row i of a 2d array view, as the Pyccel slice ind[i, :] (no copy).
  *
@@ -53,25 +53,29 @@ __device__ inline void spline_3d(double eta1, double eta2, double eta3, const lo
                                  const Array2D<long long>& ind1, const Array2D<long long>& ind2,
                                  const Array2D<long long>& ind3, const DomainArgs& args, double* f_out) {
     // mapping spans
-    int span1 = find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
-    int span2 = find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
-    int span3 = find_span(args.t3.data, args.t3.shape[0], degree[2], eta3);
+    int span1 = bsplines_kernels::find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
+    int span2 = bsplines_kernels::find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
+    int span3 = bsplines_kernels::find_span(args.t3.data, args.t3.shape[0], degree[2], eta3);
 
     // degree + 1 non-zero mapping splines
-    double b1[MAX_SPLINE_DEGREE + 1], b2[MAX_SPLINE_DEGREE + 1], b3[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];
 
-    b_splines_slim(args.t1.data, degree[0], eta1, span1, b1);
-    b_splines_slim(args.t2.data, degree[1], eta2, span2, b2);
-    b_splines_slim(args.t3.data, degree[2], eta3, span3, b3);
+    bsplines_kernels::b_splines_slim(args.t1.data, degree[0], eta1, span1, b1);
+    bsplines_kernels::b_splines_slim(args.t2.data, degree[1], eta2, span2, b2);
+    bsplines_kernels::b_splines_slim(args.t3.data, degree[2], eta3, span3, b3);
 
     // Evaluate spline mapping
     Array1D<long long> tmp1 = spline_index_row(ind1, span1 - degree[0]);
     Array1D<long long> tmp2 = spline_index_row(ind2, span2 - degree[1]);
     Array1D<long long> tmp3 = spline_index_row(ind3, span3 - degree[2]);
 
-    f_out[0] = evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, b3, tmp1, tmp2, tmp3, args.cx);
-    f_out[1] = evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, b3, tmp1, tmp2, tmp3, args.cy);
-    f_out[2] = evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, b3, tmp1, tmp2, tmp3, args.cz);
+    f_out[0] = evaluation_kernels_3d::evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, b3, tmp1, tmp2,
+                                                           tmp3, args.cx);
+    f_out[1] = evaluation_kernels_3d::evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, b3, tmp1, tmp2,
+                                                           tmp3, args.cy);
+    f_out[2] = evaluation_kernels_3d::evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, b3, tmp1, tmp2,
+                                                           tmp3, args.cz);
 }
 
 /**
@@ -93,17 +97,19 @@ __device__ inline void spline_3d_df(double eta1, double eta2, double eta3, const
                                     const Array2D<long long>& ind1, const Array2D<long long>& ind2,
                                     const Array2D<long long>& ind3, const DomainArgs& args, double* df_out) {
     // mapping spans
-    int span1 = find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
-    int span2 = find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
-    int span3 = find_span(args.t3.data, args.t3.shape[0], degree[2], eta3);
+    int span1 = bsplines_kernels::find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
+    int span2 = bsplines_kernels::find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
+    int span3 = bsplines_kernels::find_span(args.t3.data, args.t3.shape[0], degree[2], eta3);
 
     // non-zero splines of mapping, and derivatives
-    double b1[MAX_SPLINE_DEGREE + 1], b2[MAX_SPLINE_DEGREE + 1], b3[MAX_SPLINE_DEGREE + 1];
-    double der1[MAX_SPLINE_DEGREE + 1], der2[MAX_SPLINE_DEGREE + 1], der3[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];
+    double der1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], der2[bsplines_kernels::MAX_SPLINE_DEGREE + 1],
+        der3[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
 
-    b_der_splines_slim(args.t1.data, degree[0], eta1, span1, b1, der1);
-    b_der_splines_slim(args.t2.data, degree[1], eta2, span2, b2, der2);
-    b_der_splines_slim(args.t3.data, degree[2], eta3, span3, b3, der3);
+    bsplines_kernels::b_der_splines_slim(args.t1.data, degree[0], eta1, span1, b1, der1);
+    bsplines_kernels::b_der_splines_slim(args.t2.data, degree[1], eta2, span2, b2, der2);
+    bsplines_kernels::b_der_splines_slim(args.t3.data, degree[2], eta3, span3, b3, der3);
 
     // Evaluation of Jacobian
     Array1D<long long> tmp1 = spline_index_row(ind1, span1 - degree[0]);
@@ -112,9 +118,12 @@ __device__ inline void spline_3d_df(double eta1, double eta2, double eta3, const
 
     const Array3D<double>* c[3] = {&args.cx, &args.cy, &args.cz};
     for (int n = 0; n < 3; ++n) {
-        df_out[3 * n + 0] = evaluation_kernel_3d(degree[0], degree[1], degree[2], der1, b2, b3, tmp1, tmp2, tmp3, *c[n]);
-        df_out[3 * n + 1] = evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, der2, b3, tmp1, tmp2, tmp3, *c[n]);
-        df_out[3 * n + 2] = evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, der3, tmp1, tmp2, tmp3, *c[n]);
+        df_out[3 * n + 0] = evaluation_kernels_3d::evaluation_kernel_3d(degree[0], degree[1], degree[2], der1, b2, b3,
+                                                                        tmp1, tmp2, tmp3, *c[n]);
+        df_out[3 * n + 1] = evaluation_kernels_3d::evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, der2, b3,
+                                                                        tmp1, tmp2, tmp3, *c[n]);
+        df_out[3 * n + 2] = evaluation_kernels_3d::evaluation_kernel_3d(degree[0], degree[1], degree[2], b1, b2, der3,
+                                                                        tmp1, tmp2, tmp3, *c[n]);
     }
 }
 
@@ -144,21 +153,21 @@ __device__ inline void spline_2d_straight(double eta1, double eta2, double eta3,
     Array2D<double> cy = first_plane(args.cy);
 
     // mapping spans
-    int span1 = find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
-    int span2 = find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
+    int span1 = bsplines_kernels::find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
+    int span2 = bsplines_kernels::find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
 
     // degree + 1 non-zero mapping splines
-    double b1[MAX_SPLINE_DEGREE + 1], b2[MAX_SPLINE_DEGREE + 1];
+    double b1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], b2[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
 
-    b_splines_slim(args.t1.data, degree[0], eta1, span1, b1);
-    b_splines_slim(args.t2.data, degree[1], eta2, span2, b2);
+    bsplines_kernels::b_splines_slim(args.t1.data, degree[0], eta1, span1, b1);
+    bsplines_kernels::b_splines_slim(args.t2.data, degree[1], eta2, span2, b2);
 
     // Evaluate mapping
     Array1D<long long> tmp1 = spline_index_row(ind1, span1 - degree[0]);
     Array1D<long long> tmp2 = spline_index_row(ind2, span2 - degree[1]);
 
-    f_out[0] = evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx);
-    f_out[1] = evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cy);
+    f_out[0] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx);
+    f_out[1] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cy);
     f_out[2] = lz * eta3;
 
     if (eta1 == 0. && cx(0, 0) == cx(0, 1)) f_out[0] = cx(0, 0);
@@ -187,25 +196,25 @@ __device__ inline void spline_2d_straight_df(double eta1, double eta2, const lon
     Array2D<double> cy = first_plane(args.cy);
 
     // mapping spans
-    int span1 = find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
-    int span2 = find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
+    int span1 = bsplines_kernels::find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
+    int span2 = bsplines_kernels::find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
 
     // non-zero splines of mapping, and derivatives
-    double b1[MAX_SPLINE_DEGREE + 1], b2[MAX_SPLINE_DEGREE + 1];
-    double der1[MAX_SPLINE_DEGREE + 1], der2[MAX_SPLINE_DEGREE + 1];
+    double b1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], b2[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
+    double der1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], der2[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
 
-    b_der_splines_slim(args.t1.data, degree[0], eta1, span1, b1, der1);
-    b_der_splines_slim(args.t2.data, degree[1], eta2, span2, b2, der2);
+    bsplines_kernels::b_der_splines_slim(args.t1.data, degree[0], eta1, span1, b1, der1);
+    bsplines_kernels::b_der_splines_slim(args.t2.data, degree[1], eta2, span2, b2, der2);
 
     // Evaluation of Jacobian
     Array1D<long long> tmp1 = spline_index_row(ind1, span1 - degree[0]);
     Array1D<long long> tmp2 = spline_index_row(ind2, span2 - degree[1]);
 
-    df_out[0] = evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cx);
-    df_out[1] = evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cx);
+    df_out[0] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cx);
+    df_out[1] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cx);
     df_out[2] = 0.;
-    df_out[3] = evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cy);
-    df_out[4] = evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cy);
+    df_out[3] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cy);
+    df_out[4] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cy);
     df_out[5] = 0.;
     df_out[6] = 0.;
     df_out[7] = 0.;
@@ -239,23 +248,24 @@ __device__ inline void spline_2d_torus(double eta1, double eta2, double eta3, co
     Array2D<double> cy = first_plane(args.cy);
 
     // mapping spans
-    int span1 = find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
-    int span2 = find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
+    int span1 = bsplines_kernels::find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
+    int span2 = bsplines_kernels::find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
 
     // degree + 1 non-zero mapping splines
-    double b1[MAX_SPLINE_DEGREE + 1], b2[MAX_SPLINE_DEGREE + 1];
+    double b1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], b2[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
 
-    b_splines_slim(args.t1.data, degree[0], eta1, span1, b1);
-    b_splines_slim(args.t2.data, degree[1], eta2, span2, b2);
+    bsplines_kernels::b_splines_slim(args.t1.data, degree[0], eta1, span1, b1);
+    bsplines_kernels::b_splines_slim(args.t2.data, degree[1], eta2, span2, b2);
 
     // Evaluate mapping
     Array1D<long long> tmp1 = spline_index_row(ind1, span1 - degree[0]);
     Array1D<long long> tmp2 = spline_index_row(ind2, span2 - degree[1]);
 
-    f_out[0] = evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) * cos(2 * pi * eta3 / tor_period);
-    f_out[1] = evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) * (-1) *
+    f_out[0] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) *
+               cos(2 * pi * eta3 / tor_period);
+    f_out[1] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) * (-1) *
                sin(2 * pi * eta3 / tor_period);
-    f_out[2] = evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cy);
+    f_out[2] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cy);
 
     if (eta1 == 0. && cx(0, 0) == cx(0, 1)) {
         f_out[0] = cx(0, 0) * cos(2 * pi * eta3 / tor_period);
@@ -285,31 +295,33 @@ __device__ inline void spline_2d_torus_df(double eta1, double eta2, double eta3,
     Array2D<double> cy = first_plane(args.cy);
 
     // mapping spans
-    int span1 = find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
-    int span2 = find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
+    int span1 = bsplines_kernels::find_span(args.t1.data, args.t1.shape[0], degree[0], eta1);
+    int span2 = bsplines_kernels::find_span(args.t2.data, args.t2.shape[0], degree[1], eta2);
 
     // non-zero splines of mapping, and derivatives
-    double b1[MAX_SPLINE_DEGREE + 1], b2[MAX_SPLINE_DEGREE + 1];
-    double der1[MAX_SPLINE_DEGREE + 1], der2[MAX_SPLINE_DEGREE + 1];
+    double b1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], b2[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
+    double der1[bsplines_kernels::MAX_SPLINE_DEGREE + 1], der2[bsplines_kernels::MAX_SPLINE_DEGREE + 1];
 
-    b_der_splines_slim(args.t1.data, degree[0], eta1, span1, b1, der1);
-    b_der_splines_slim(args.t2.data, degree[1], eta2, span2, b2, der2);
+    bsplines_kernels::b_der_splines_slim(args.t1.data, degree[0], eta1, span1, b1, der1);
+    bsplines_kernels::b_der_splines_slim(args.t2.data, degree[1], eta2, span2, b2, der2);
 
     Array1D<long long> tmp1 = spline_index_row(ind1, span1 - degree[0]);
     Array1D<long long> tmp2 = spline_index_row(ind2, span2 - degree[1]);
 
-    df_out[0] = evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cx) * cos(2 * pi * eta3 / tor_period);
-    df_out[1] = evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cx) * cos(2 * pi * eta3 / tor_period);
-    df_out[2] = evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) * sin(2 * pi * eta3 / tor_period) *
-                (-2 * pi / tor_period);
-    df_out[3] = evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cx) * (-1) *
+    df_out[0] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cx) *
+                cos(2 * pi * eta3 / tor_period);
+    df_out[1] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cx) *
+                cos(2 * pi * eta3 / tor_period);
+    df_out[2] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) *
+                sin(2 * pi * eta3 / tor_period) * (-2 * pi / tor_period);
+    df_out[3] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cx) * (-1) *
                 sin(2 * pi * eta3 / tor_period);
-    df_out[4] = evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cx) * (-1) *
+    df_out[4] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cx) * (-1) *
                 sin(2 * pi * eta3 / tor_period);
-    df_out[5] = evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) * (-1) *
+    df_out[5] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, b2, tmp1, tmp2, cx) * (-1) *
                 cos(2 * pi * eta3 / tor_period) * 2 * pi / tor_period;
-    df_out[6] = evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cy);
-    df_out[7] = evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cy);
+    df_out[6] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], der1, b2, tmp1, tmp2, cy);
+    df_out[7] = evaluation_kernels_2d::evaluation_kernel_2d(degree[0], degree[1], b1, der2, tmp1, tmp2, cy);
     df_out[8] = 0.;
 
     if (eta1 == 0. && cx(0, 0) == cx(0, 1)) {
@@ -319,4 +331,4 @@ __device__ inline void spline_2d_torus_df(double eta1, double eta2, double eta3,
 
     if (eta1 == 0. && cy(0, 0) == cy(0, 1)) df_out[7] = 0.;
 }
-}
+}  // namespace struphy_cuda::spline_mappings_kernels
diff --git a/src/struphy/geometry/tests/spline_mapping_cases.py b/src/struphy/geometry/tests/spline_mapping_cases.py
index a85819142..4f5dddde9 100644
--- a/src/struphy/geometry/tests/spline_mapping_cases.py
+++ b/src/struphy/geometry/tests/spline_mapping_cases.py
@@ -16,21 +16,21 @@
     double f_out[3], df_out[9];
     switch (args_domain.kind_map) {
         case 0:
-            struphy_cuda::spline_3d(eta1, eta2, eta3, args_domain.degree, args_domain.ind1, args_domain.ind2,
+            struphy_cuda::spline_mappings_kernels::spline_3d(eta1, eta2, eta3, args_domain.degree, args_domain.ind1, args_domain.ind2,
                                     args_domain.ind3, args_domain, f_out);
-            struphy_cuda::spline_3d_df(eta1, eta2, eta3, args_domain.degree, args_domain.ind1, args_domain.ind2,
+            struphy_cuda::spline_mappings_kernels::spline_3d_df(eta1, eta2, eta3, args_domain.degree, args_domain.ind1, args_domain.ind2,
                                        args_domain.ind3, args_domain, df_out);
             break;
         case 1:
-            struphy_cuda::spline_2d_straight(eta1, eta2, eta3, args_domain.degree, args_domain.ind1,
+            struphy_cuda::spline_mappings_kernels::spline_2d_straight(eta1, eta2, eta3, args_domain.degree, args_domain.ind1,
                                              args_domain.ind2, args_domain, args_domain.params[0], f_out);
-            struphy_cuda::spline_2d_straight_df(eta1, eta2, args_domain.degree, args_domain.ind1, args_domain.ind2,
+            struphy_cuda::spline_mappings_kernels::spline_2d_straight_df(eta1, eta2, args_domain.degree, args_domain.ind1, args_domain.ind2,
                                                 args_domain, args_domain.params[0], df_out);
             break;
         case 2:
-            struphy_cuda::spline_2d_torus(eta1, eta2, eta3, args_domain.degree, args_domain.ind1, args_domain.ind2,
+            struphy_cuda::spline_mappings_kernels::spline_2d_torus(eta1, eta2, eta3, args_domain.degree, args_domain.ind1, args_domain.ind2,
                                           args_domain, args_domain.params[0], f_out);
-            struphy_cuda::spline_2d_torus_df(eta1, eta2, eta3, args_domain.degree, args_domain.ind1,
+            struphy_cuda::spline_mappings_kernels::spline_2d_torus_df(eta1, eta2, eta3, args_domain.degree, args_domain.ind1,
                                              args_domain.ind2, args_domain, args_domain.params[0], df_out);
             break;
         default:
@@ -148,10 +148,10 @@ def host_domain(name):
 extern "C" __global__ void evaluate_der_splines(const double* t, int nt, int p, const double* x, double* out, int n) {
     int i = blockDim.x * blockIdx.x + threadIdx.x;
     if (i >= n) return;
-    int span = struphy_cuda::find_span(t, nt, p, x[i]);
+    int span = struphy_cuda::bsplines_kernels::find_span(t, nt, p, x[i]);
     double* row = out + i * (3 * p + 3);
-    struphy_cuda::b_splines_slim(t, p, x[i], span, row);
-    struphy_cuda::b_der_splines_slim(t, p, x[i], span, row + p + 1, row + 2 * p + 2);
+    struphy_cuda::bsplines_kernels::b_splines_slim(t, p, x[i], span, row);
+    struphy_cuda::bsplines_kernels::b_der_splines_slim(t, p, x[i], span, row + p + 1, row + 2 * p + 2);
 }
 """
 
diff --git a/src/struphy/geometry/transform_kernels.cuh b/src/struphy/geometry/transform_kernels.cuh
index 204535bbd..31115e5e1 100644
--- a/src/struphy/geometry/transform_kernels.cuh
+++ b/src/struphy/geometry/transform_kernels.cuh
@@ -5,7 +5,7 @@
 // Vector-valued kinds (kind_fun >= 10) read and write three entries of a and out; scalar kinds only a[0], out[0].
 #include "struphy/geometry/evaluation_kernels.cuh"
 #include "struphy/linear_algebra/linalg_kernels.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::transform_kernels {
 /**
  * Pull-back of a Cartesian scalar/vector field to a differential p-form, as in transform_kernels.pull.
  *
@@ -24,8 +24,8 @@ __device__ inline void pull(const double* a, double eta1, double eta2, double et
 
     // evaluate Jacobian matrix and its determinant
     if (kind_fun > 0) {
-        df(eta1, eta2, eta3, args_domain, dfmat1);
-        detdf = det(dfmat1);
+        evaluation_kernels::df(eta1, eta2, eta3, args_domain, dfmat1);
+        detdf = linalg_kernels::det(dfmat1);
     }
 
     // 0-form
@@ -38,21 +38,21 @@ __device__ inline void pull(const double* a, double eta1, double eta2, double et
     }
     // 1-form
     else if (kind_fun == 10) {
-        transpose(dfmat1, dfmat2);
-        matrix_vector(dfmat2, a, out);
+        linalg_kernels::transpose(dfmat1, dfmat2);
+        linalg_kernels::matrix_vector(dfmat2, a, out);
     }
     // 2-form
     else if (kind_fun == 11) {
-        matrix_inv_with_det(dfmat1, 1.0, dfmat2);
-        matrix_vector(dfmat2, a, out);
+        linalg_kernels::matrix_inv_with_det(dfmat1, 1.0, dfmat2);
+        linalg_kernels::matrix_vector(dfmat2, a, out);
 
         if (detdf < 0.0)
             for (int i = 0; i < 3; ++i) out[i] = -out[i];
     }
     // vector
     else if (kind_fun == 12) {
-        matrix_inv(dfmat1, dfmat2);
-        matrix_vector(dfmat2, a, out);
+        linalg_kernels::matrix_inv(dfmat1, dfmat2);
+        linalg_kernels::matrix_vector(dfmat2, a, out);
     }
 }
 
@@ -74,8 +74,8 @@ __device__ inline void push(const double* a, double eta1, double eta2, double et
 
     // evaluate Jacobian matrix and its determinant
     if (kind_fun > 0) {
-        df(eta1, eta2, eta3, args_domain, dfmat1);
-        detdf = det(dfmat1);
+        evaluation_kernels::df(eta1, eta2, eta3, args_domain, dfmat1);
+        detdf = linalg_kernels::det(dfmat1);
     }
 
     // 0-form
@@ -88,18 +88,18 @@ __device__ inline void push(const double* a, double eta1, double eta2, double et
     }
     // 1-form
     else if (kind_fun == 10) {
-        matrix_inv_with_det(dfmat1, detdf, dfmat2);
-        transpose(dfmat2, dfmat3);
-        matrix_vector(dfmat3, a, out);
+        linalg_kernels::matrix_inv_with_det(dfmat1, detdf, dfmat2);
+        linalg_kernels::transpose(dfmat2, dfmat3);
+        linalg_kernels::matrix_vector(dfmat3, a, out);
     }
     // 2-form
     else if (kind_fun == 11) {
-        matrix_vector(dfmat1, a, out);
+        linalg_kernels::matrix_vector(dfmat1, a, out);
         for (int i = 0; i < 3; ++i) out[i] = out[i] / fabs(detdf);
     }
     // vector
     else if (kind_fun == 12) {
-        matrix_vector(dfmat1, a, out);
+        linalg_kernels::matrix_vector(dfmat1, a, out);
     }
 }
 
@@ -135,8 +135,8 @@ __device__ inline void tran(const double* a, double eta1, double eta2, double et
     double vec1[3], vec2[3];
 
     // evaluate Jacobian matrix and its determinant
-    df(eta1, eta2, eta3, args_domain, dfmat1);
-    double detdf = det(dfmat1);
+    evaluation_kernels::df(eta1, eta2, eta3, args_domain, dfmat1);
+    double detdf = linalg_kernels::det(dfmat1);
 
     // 0-form to 3-form
     if (kind_fun == 0) {
@@ -148,17 +148,17 @@ __device__ inline void tran(const double* a, double eta1, double eta2, double et
     }
     // 1-form to 2-form (a^2 = G^(-1) * a^1 * |det(DF)|)
     else if (kind_fun == 10) {
-        matrix_inv_with_det(dfmat1, detdf, dfmat2);
-        transpose(dfmat2, dfmat3);
-        matrix_vector(dfmat3, a, vec1);
-        matrix_vector(dfmat2, vec1, out);
+        linalg_kernels::matrix_inv_with_det(dfmat1, detdf, dfmat2);
+        linalg_kernels::transpose(dfmat2, dfmat3);
+        linalg_kernels::matrix_vector(dfmat3, a, vec1);
+        linalg_kernels::matrix_vector(dfmat2, vec1, out);
         for (int i = 0; i < 3; ++i) out[i] = out[i] * fabs(detdf);
     }
     // 2-form to 1-form (a^1 = G * a^2 / |det(DF)|)
     else if (kind_fun == 11) {
-        transpose(dfmat1, dfmat2);
-        matrix_vector(dfmat1, a, vec1);
-        matrix_vector(dfmat2, vec1, out);
+        linalg_kernels::transpose(dfmat1, dfmat2);
+        linalg_kernels::matrix_vector(dfmat1, a, vec1);
+        linalg_kernels::matrix_vector(dfmat2, vec1, out);
         for (int i = 0; i < 3; ++i) out[i] = out[i] / fabs(detdf);
     }
     // norm vector to vector
@@ -168,9 +168,9 @@ __device__ inline void tran(const double* a, double eta1, double eta2, double et
     // norm vector to 1-form (a^1 = G * a)
     else if (kind_fun == 13) {
         for (int i = 0; i < 3; ++i) vec1[i] = a[i] / column_norm(dfmat1, i);
-        transpose(dfmat1, dfmat2);
-        matrix_vector(dfmat1, vec1, vec2);
-        matrix_vector(dfmat2, vec2, out);
+        linalg_kernels::transpose(dfmat1, dfmat2);
+        linalg_kernels::matrix_vector(dfmat1, vec1, vec2);
+        linalg_kernels::matrix_vector(dfmat2, vec2, out);
     }
     // norm vector to 2-form (a^2 = |det(DF)| * a)
     else if (kind_fun == 14) {
@@ -179,9 +179,9 @@ __device__ inline void tran(const double* a, double eta1, double eta2, double et
     }
     // vector to 1-form (a^1 = G * a)
     else if (kind_fun == 15) {
-        transpose(dfmat1, dfmat2);
-        matrix_vector(dfmat1, a, vec1);
-        matrix_vector(dfmat2, vec1, out);
+        linalg_kernels::transpose(dfmat1, dfmat2);
+        linalg_kernels::matrix_vector(dfmat1, a, vec1);
+        linalg_kernels::matrix_vector(dfmat2, vec1, out);
     }
     // vector to 2-form (a^2 = |det(DF)| * a)
     else if (kind_fun == 16) {
@@ -189,10 +189,10 @@ __device__ inline void tran(const double* a, double eta1, double eta2, double et
     }
     // 1-form to vector (a = G^(-1) * a^1)
     else if (kind_fun == 17) {
-        matrix_inv_with_det(dfmat1, detdf, dfmat2);
-        transpose(dfmat2, dfmat3);
-        matrix_vector(dfmat3, a, vec1);
-        matrix_vector(dfmat2, vec1, out);
+        linalg_kernels::matrix_inv_with_det(dfmat1, detdf, dfmat2);
+        linalg_kernels::transpose(dfmat2, dfmat3);
+        linalg_kernels::matrix_vector(dfmat3, a, vec1);
+        linalg_kernels::matrix_vector(dfmat2, vec1, out);
     }
     // 2-form to vector (a = a^2 / |det(DF)|)
     else if (kind_fun == 18) {
@@ -200,10 +200,10 @@ __device__ inline void tran(const double* a, double eta1, double eta2, double et
     }
     // 1-form to normalized vector
     else if (kind_fun == 19) {
-        matrix_inv_with_det(dfmat1, detdf, dfmat2);
-        transpose(dfmat2, dfmat3);
-        matrix_vector(dfmat3, a, vec1);
-        matrix_vector(dfmat2, vec1, out);
+        linalg_kernels::matrix_inv_with_det(dfmat1, detdf, dfmat2);
+        linalg_kernels::transpose(dfmat2, dfmat3);
+        linalg_kernels::matrix_vector(dfmat3, a, vec1);
+        linalg_kernels::matrix_vector(dfmat2, vec1, out);
         for (int i = 0; i < 3; ++i) out[i] = out[i] * column_norm(dfmat1, i);
     }
     // 2-form to normalized vector
@@ -216,4 +216,4 @@ __device__ inline void tran(const double* a, double eta1, double eta2, double et
         for (int i = 0; i < 3; ++i) out[i] = a[i] * column_norm(dfmat1, i);
     }
 }
-}
+}  // namespace struphy_cuda::transform_kernels
diff --git a/src/struphy/linear_algebra/linalg_kernels.cuh b/src/struphy/linear_algebra/linalg_kernels.cuh
index 4636f2ea3..af429c5cd 100644
--- a/src/struphy/linear_algebra/linalg_kernels.cuh
+++ b/src/struphy/linear_algebra/linalg_kernels.cuh
@@ -1,5 +1,5 @@
 #pragma once
-namespace struphy_cuda {
+namespace struphy_cuda::linalg_kernels {
 /**
  * Multiply a 3x3 matrix by a vector, as in linalg_kernels.matrix_vector.
  *
@@ -111,4 +111,4 @@ __device__ inline void outer(const double* a, const double* b, double* c) {
     for (int i = 0; i < 3; ++i)
         for (int j = 0; j < 3; ++j) c[3 * i + j] = a[i] * b[j];
 }
-}
+}  // namespace struphy_cuda::linalg_kernels
diff --git a/src/struphy/pic/accumulation/filler_kernels.cuh b/src/struphy/pic/accumulation/filler_kernels.cuh
index 8a33d7773..49355028d 100644
--- a/src/struphy/pic/accumulation/filler_kernels.cuh
+++ b/src/struphy/pic/accumulation/filler_kernels.cuh
@@ -3,7 +3,7 @@
 #pragma once
 #include "cunumpy/array_view.cuh"
 #include "cunumpy/atomic.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::filler_kernels {
 /**
  * Add basis functions times filling to a stencil vector, as in filler_kernels.fill_vec.
  *
@@ -161,4 +161,4 @@ __device__ inline void fill_mat_vec(int pi1, int pi2, int pi3, int pj1, int pj2,
         }
     }
 }
-}
+}  // namespace struphy_cuda::filler_kernels
diff --git a/src/struphy/pic/accumulation/kernels/charge_density_0form/charge_density_0form_cuda.cu b/src/struphy/pic/accumulation/kernels/charge_density_0form/charge_density_0form_cuda.cu
index 014093a1b..b03038867 100644
--- a/src/struphy/pic/accumulation/kernels/charge_density_0form/charge_density_0form_cuda.cu
+++ b/src/struphy/pic/accumulation/kernels/charge_density_0form/charge_density_0form_cuda.cu
@@ -1,5 +1,7 @@
 #include "struphy/pic/accumulation/particle_to_mat_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Accumulate the marker weights into a 0-form vector, as in charge_density_0form_kernels.charge_density_0form.
  *
@@ -32,5 +34,5 @@ extern "C" __global__ void charge_density_0form(MarkerArgs args_markers, DerhamA
     // filling is just the weights
     double filling = args_markers.markers(ip, weight_idx);
 
-    struphy_cuda::vec_fill_b_v0(args_derham, eta1, eta2, eta3, vec, filling);
+    particle_to_mat_kernels::vec_fill_b_v0(args_derham, eta1, eta2, eta3, vec, filling);
 }
diff --git a/src/struphy/pic/accumulation/kernels/linear_vlasov_ampere/linear_vlasov_ampere_cuda.cu b/src/struphy/pic/accumulation/kernels/linear_vlasov_ampere/linear_vlasov_ampere_cuda.cu
index 387013af2..e4ce83e46 100644
--- a/src/struphy/pic/accumulation/kernels/linear_vlasov_ampere/linear_vlasov_ampere_cuda.cu
+++ b/src/struphy/pic/accumulation/kernels/linear_vlasov_ampere/linear_vlasov_ampere_cuda.cu
@@ -2,6 +2,8 @@
 #include "struphy/linear_algebra/linalg_kernels.cuh"
 #include "struphy/pic/accumulation/particle_to_mat_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Accumulate the matrix and vector of the linear Vlasov-Ampere system (delta-f) into V1, as in
  * linear_vlasov_ampere_kernels.linear_vlasov_ampere.
@@ -67,16 +69,16 @@ extern "C" __global__ void linear_vlasov_ampere(MarkerArgs args_markers, DerhamA
     v[2] = args_markers.markers(ip, 5);
 
     // evaluate Jacobian, result in dfm
-    struphy_cuda::df(eta1, eta2, eta3, args_domain, dfm);
+    evaluation_kernels::df(eta1, eta2, eta3, args_domain, dfm);
 
     // invert Jacobian matrix
-    struphy_cuda::matrix_inv(dfm, df_inv);
+    linalg_kernels::matrix_inv(dfm, df_inv);
 
     // compute DF^{-1} v
-    struphy_cuda::matrix_vector(df_inv, v, df_inv_v);
+    linalg_kernels::matrix_vector(df_inv, v, df_inv_v);
 
     // filling_m = alpha^2 * kappa^2 * f0 / (N * s_0 * v_th^2) * (DF^{-1} v_p)_mu * (DF^{-1} v_p)_nu
-    struphy_cuda::outer(df_inv_v, df_inv_v, filling_m);
+    linalg_kernels::outer(df_inv_v, df_inv_v, filling_m);
     double factor = f0_values[ip] / (args_markers.markers(ip, 7) * n_markers_tot);
     for (int k = 0; k < 9; ++k) filling_m[k] *= factor;
 
@@ -84,7 +86,8 @@ extern "C" __global__ void linear_vlasov_ampere(MarkerArgs args_markers, DerhamA
     for (int k = 0; k < 3; ++k) filling_v[k] = args_markers.markers(ip, 6) * df_inv_v[k];
 
     // call the appropriate matvec filler
-    struphy_cuda::m_v_fill_b_v1_symm(args_derham, eta1, eta2, eta3, mat11, mat12, mat13, mat22, mat23, mat33,
-                                     filling_m[0], filling_m[1], filling_m[2], filling_m[4], filling_m[5],
-                                     filling_m[8], vec1, vec2, vec3, filling_v[0], filling_v[1], filling_v[2]);
+    particle_to_mat_kernels::m_v_fill_b_v1_symm(args_derham, eta1, eta2, eta3, mat11, mat12, mat13, mat22, mat23,
+                                                mat33, filling_m[0], filling_m[1], filling_m[2], filling_m[4],
+                                                filling_m[5], filling_m[8], vec1, vec2, vec3, filling_v[0],
+                                                filling_v[1], filling_v[2]);
 }
diff --git a/src/struphy/pic/accumulation/kernels/vlasov_maxwell/vlasov_maxwell_cuda.cu b/src/struphy/pic/accumulation/kernels/vlasov_maxwell/vlasov_maxwell_cuda.cu
index e6abf9bae..880538d1e 100644
--- a/src/struphy/pic/accumulation/kernels/vlasov_maxwell/vlasov_maxwell_cuda.cu
+++ b/src/struphy/pic/accumulation/kernels/vlasov_maxwell/vlasov_maxwell_cuda.cu
@@ -2,6 +2,8 @@
 #include "struphy/linear_algebra/linalg_kernels.cuh"
 #include "struphy/pic/accumulation/particle_to_mat_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Accumulate the matrix and vector of the full-orbit Vlasov-Maxwell (Vlasov-Ampere) coupling into V1, as in
  * vlasov_maxwell_kernels.vlasov_maxwell.
@@ -52,7 +54,7 @@ extern "C" __global__ void vlasov_maxwell(MarkerArgs args_markers, DerhamArgs ar
     double eta3 = args_markers.markers(ip, 2);
 
     // evaluate Jacobian, result in dfm
-    struphy_cuda::df(eta1, eta2, eta3, args_domain, dfm);
+    evaluation_kernels::df(eta1, eta2, eta3, args_domain, dfm);
 
     // compute shifted and stretched velocity
     v[0] = args_markers.markers(ip, 3);
@@ -60,10 +62,10 @@ extern "C" __global__ void vlasov_maxwell(MarkerArgs args_markers, DerhamArgs ar
     v[2] = args_markers.markers(ip, 5);
 
     // filling functions
-    struphy_cuda::matrix_inv(dfm, df_inv);
-    struphy_cuda::transpose(df_inv, df_inv_t);
-    struphy_cuda::matrix_matrix(df_inv, df_inv_t, g_inv);
-    struphy_cuda::matrix_vector(df_inv, v, df_inv_times_v);
+    linalg_kernels::matrix_inv(dfm, df_inv);
+    linalg_kernels::transpose(df_inv, df_inv_t);
+    linalg_kernels::matrix_matrix(df_inv, df_inv_t, g_inv);
+    linalg_kernels::matrix_vector(df_inv, v, df_inv_times_v);
 
     // filling_m = w_p * DF^{-1} * DF^{-T}
     for (int k = 0; k < 9; ++k) filling_m[k] = args_markers.markers(ip, 6) * g_inv[k];
@@ -72,7 +74,8 @@ extern "C" __global__ void vlasov_maxwell(MarkerArgs args_markers, DerhamArgs ar
     for (int k = 0; k < 3; ++k) filling_v[k] = args_markers.markers(ip, 6) * df_inv_times_v[k];
 
     // call the appropriate matvec filler
-    struphy_cuda::m_v_fill_b_v1_symm(args_derham, eta1, eta2, eta3, mat11, mat12, mat13, mat22, mat23, mat33,
-                                     filling_m[0], filling_m[1], filling_m[2], filling_m[4], filling_m[5],
-                                     filling_m[8], vec1, vec2, vec3, filling_v[0], filling_v[1], filling_v[2]);
+    particle_to_mat_kernels::m_v_fill_b_v1_symm(args_derham, eta1, eta2, eta3, mat11, mat12, mat13, mat22, mat23,
+                                                mat33, filling_m[0], filling_m[1], filling_m[2], filling_m[4],
+                                                filling_m[5], filling_m[8], vec1, vec2, vec3, filling_v[0],
+                                                filling_v[1], filling_v[2]);
 }
diff --git a/src/struphy/pic/accumulation/particle_to_mat_kernels.cuh b/src/struphy/pic/accumulation/particle_to_mat_kernels.cuh
index ade422b97..c587e40cd 100644
--- a/src/struphy/pic/accumulation/particle_to_mat_kernels.cuh
+++ b/src/struphy/pic/accumulation/particle_to_mat_kernels.cuh
@@ -2,7 +2,7 @@
 #pragma once
 #include "struphy/bsplines/evaluation_kernels_3d.cuh"
 #include "struphy/pic/accumulation/filler_kernels.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::particle_to_mat_kernels {
 /**
  * Add one marker's contribution to a 0-form vector, as in particle_to_mat_kernels.vec_fill_b_v0.
  *
@@ -21,12 +21,12 @@ __device__ inline void vec_fill_b_v0(const DerhamArgs& args_derham, double eta1,
     int pn3 = args_derham.pn[2];
 
     // CUDA-only scratch holds the spans and the spline values pyccel keeps in args_derham
-    SplineScratch scratch;
-    get_spans(eta1, eta2, eta3, args_derham, scratch);
+    evaluation_kernels_3d::SplineScratch scratch;
+    evaluation_kernels_3d::get_spans(eta1, eta2, eta3, args_derham, scratch);
     int span1 = scratch.span1, span2 = scratch.span2, span3 = scratch.span3;
 
-    fill_vec(pn1, pn2, pn3, scratch.bn1, scratch.bn2, scratch.bn3, span1, span2, span3, args_derham.starts, vec,
-             fill);
+    filler_kernels::fill_vec(pn1, pn2, pn3, scratch.bn1, scratch.bn2, scratch.bn3, span1, span2, span3,
+                             args_derham.starts, vec, fill);
 }
 
 /**
@@ -74,27 +74,33 @@ __device__ inline void m_v_fill_b_v1_symm(const DerhamArgs& args_derham, double
     int pd3 = pn3 - 1;
 
     // CUDA-only scratch holds the spans and the spline values pyccel keeps in args_derham
-    SplineScratch scratch;
-    get_spans(eta1, eta2, eta3, args_derham, scratch);
+    evaluation_kernels_3d::SplineScratch scratch;
+    evaluation_kernels_3d::get_spans(eta1, eta2, eta3, args_derham, scratch);
     int span1 = scratch.span1, span2 = scratch.span2, span3 = scratch.span3;
 
     // fill matrix entries
-    fill_mat_vec(pd1, pn2, pn3, pd1, pn2, pn3, scratch.bd1, scratch.bn2, scratch.bn3, scratch.bd1, scratch.bn2,
-                 scratch.bn3, span1, span2, span3, args_derham.starts, args_derham.pn, mat11, fill11, vec1, fill1);
+    filler_kernels::fill_mat_vec(pd1, pn2, pn3, pd1, pn2, pn3, scratch.bd1, scratch.bn2, scratch.bn3, scratch.bd1,
+                                 scratch.bn2, scratch.bn3, span1, span2, span3, args_derham.starts, args_derham.pn,
+                                 mat11, fill11, vec1, fill1);
 
-    fill_mat_vec(pn1, pd2, pn3, pn1, pd2, pn3, scratch.bn1, scratch.bd2, scratch.bn3, scratch.bn1, scratch.bd2,
-                 scratch.bn3, span1, span2, span3, args_derham.starts, args_derham.pn, mat22, fill22, vec2, fill2);
+    filler_kernels::fill_mat_vec(pn1, pd2, pn3, pn1, pd2, pn3, scratch.bn1, scratch.bd2, scratch.bn3, scratch.bn1,
+                                 scratch.bd2, scratch.bn3, span1, span2, span3, args_derham.starts, args_derham.pn,
+                                 mat22, fill22, vec2, fill2);
 
-    fill_mat_vec(pn1, pn2, pd3, pn1, pn2, pd3, scratch.bn1, scratch.bn2, scratch.bd3, scratch.bn1, scratch.bn2,
-                 scratch.bd3, span1, span2, span3, args_derham.starts, args_derham.pn, mat33, fill33, vec3, fill3);
+    filler_kernels::fill_mat_vec(pn1, pn2, pd3, pn1, pn2, pd3, scratch.bn1, scratch.bn2, scratch.bd3, scratch.bn1,
+                                 scratch.bn2, scratch.bd3, span1, span2, span3, args_derham.starts, args_derham.pn,
+                                 mat33, fill33, vec3, fill3);
 
-    fill_mat(pd1, pn2, pn3, pn1, pd2, pn3, scratch.bd1, scratch.bn2, scratch.bn3, scratch.bn1, scratch.bd2,
-             scratch.bn3, span1, span2, span3, args_derham.starts, args_derham.pn, mat12, fill12);
+    filler_kernels::fill_mat(pd1, pn2, pn3, pn1, pd2, pn3, scratch.bd1, scratch.bn2, scratch.bn3, scratch.bn1,
+                             scratch.bd2, scratch.bn3, span1, span2, span3, args_derham.starts, args_derham.pn, mat12,
+                             fill12);
 
-    fill_mat(pd1, pn2, pn3, pn1, pn2, pd3, scratch.bd1, scratch.bn2, scratch.bn3, scratch.bn1, scratch.bn2,
-             scratch.bd3, span1, span2, span3, args_derham.starts, args_derham.pn, mat13, fill13);
+    filler_kernels::fill_mat(pd1, pn2, pn3, pn1, pn2, pd3, scratch.bd1, scratch.bn2, scratch.bn3, scratch.bn1,
+                             scratch.bn2, scratch.bd3, span1, span2, span3, args_derham.starts, args_derham.pn, mat13,
+                             fill13);
 
-    fill_mat(pn1, pd2, pn3, pn1, pn2, pd3, scratch.bn1, scratch.bd2, scratch.bn3, scratch.bn1, scratch.bn2,
-             scratch.bd3, span1, span2, span3, args_derham.starts, args_derham.pn, mat23, fill23);
-}
+    filler_kernels::fill_mat(pn1, pd2, pn3, pn1, pn2, pd3, scratch.bn1, scratch.bd2, scratch.bn3, scratch.bn1,
+                             scratch.bn2, scratch.bd3, span1, span2, span3, args_derham.starts, args_derham.pn, mat23,
+                             fill23);
 }
+}  // namespace struphy_cuda::particle_to_mat_kernels
diff --git a/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_cuda.cu b/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_cuda.cu
index 7a473ed1b..7a62cb4c8 100644
--- a/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_cuda.cu
@@ -2,6 +2,8 @@
 #include "struphy/geometry/evaluation_kernels.cuh"
 #include "struphy/linear_algebra/linalg_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Evaluate B*_parallel as a 3-form, as in bstar_parallel_3form_kernels.bstar_parallel_3form.
  *
@@ -52,20 +54,21 @@ extern "C" __global__ void bstar_parallel_3form(const double* alpha, const long
     double v = alpha[3] * v_k + (1.0 - alpha[3]) * v_n;
 
     // evaluate Jacobian, result in dfm
-    struphy_cuda::df(eta[0], eta[1], eta[2], args_domain, dfm);
+    evaluation_kernels::df(eta[0], eta[1], eta[2], args_domain, dfm);
 
-    double det_df = struphy_cuda::det(dfm);
+    double det_df = linalg_kernels::det(dfm);
 
     // spline evaluation; CUDA-only scratch holds the spline values pyccel keeps in args_derham
-    struphy_cuda::SplineScratch scratch;
-    struphy_cuda::get_spans(eta[0], eta[1], eta[2], args_derham, scratch);
+    evaluation_kernels_3d::SplineScratch scratch;
+    evaluation_kernels_3d::get_spans(eta[0], eta[1], eta[2], args_derham, scratch);
     int span1 = scratch.span1, span2 = scratch.span2, span3 = scratch.span3;
 
     // compute B*_parallel
-    double B_dot_b = struphy_cuda::eval_0form_spline_mpi(span1, span2, span3, args_derham, scratch, B_dot_b_coeffs);
+    double B_dot_b =
+        evaluation_kernels_3d::eval_0form_spline_mpi(span1, span2, span3, args_derham, scratch, B_dot_b_coeffs);
 
     double b_star_parallel =
-        struphy_cuda::eval_0form_spline_mpi(span1, span2, span3, args_derham, scratch, curl_unit_b_dot_b0);
+        evaluation_kernels_3d::eval_0form_spline_mpi(span1, span2, span3, args_derham, scratch, curl_unit_b_dot_b0);
 
     b_star_parallel *= epsilon * v;
     b_star_parallel += B_dot_b;
diff --git a/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu b/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu
index 0f5a1938c..9a6c24d52 100644
--- a/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu
@@ -1,5 +1,7 @@
 #include "struphy/pic/pushing/pusher_utilities_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * One stage of an explicit Runge-Kutta position push, as in push_eta_stage_kernels.push_eta_stage.
  *
@@ -39,13 +41,13 @@ extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_
     for (int j = 0; j < 3; ++j) v[j] = args_markers.markers(ip, 3 + j);
 
     // evaluate Jacobian, result in dfm
-    struphy_cuda::df(e1, e2, e3, args_domain, dfm);
+    evaluation_kernels::df(e1, e2, e3, args_domain, dfm);
 
     // evaluate inverse Jacobian matrix
-    struphy_cuda::matrix_inv(dfm, dfinv);
+    linalg_kernels::matrix_inv(dfm, dfinv);
 
     // pull-back of velocity
-    struphy_cuda::matrix_vector(dfinv, v, k);
+    linalg_kernels::matrix_vector(dfinv, v, k);
 
     for (int j = 0; j < 3; ++j) {
         // accumulation for last stage
@@ -56,5 +58,5 @@ extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_
     }
 
     // apply kinetic boundary conditions
-    struphy_cuda::apply_kinetic_bc_marker(ip, args_markers, args_domain, false);
+    pusher_utilities_kernels::apply_kinetic_bc_marker(ip, args_markers, args_domain, false);
 }
diff --git a/src/struphy/pic/pushing/kernels/push_v_with_efield/push_v_with_efield_cuda.cu b/src/struphy/pic/pushing/kernels/push_v_with_efield/push_v_with_efield_cuda.cu
index 31a6eff1c..e5ede6bf2 100644
--- a/src/struphy/pic/pushing/kernels/push_v_with_efield/push_v_with_efield_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_v_with_efield/push_v_with_efield_cuda.cu
@@ -2,6 +2,8 @@
 #include "struphy/geometry/evaluation_kernels.cuh"
 #include "struphy/linear_algebra/linalg_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Velocity update in an electric field, as in push_v_with_efield_kernels.push_v_with_efield.
  *
@@ -36,22 +38,22 @@ extern "C" __global__ void push_v_with_efield(double dt, int stage, MarkerArgs a
     double dfm[9], dfinv[9], dfinvt[9], e_form[3], e_cart[3];
 
     // evaluate Jacobian, result in dfm
-    struphy_cuda::df(eta1, eta2, eta3, args_domain, dfm);
+    evaluation_kernels::df(eta1, eta2, eta3, args_domain, dfm);
 
     // metric coeffs
-    struphy_cuda::matrix_inv(dfm, dfinv);
-    struphy_cuda::transpose(dfinv, dfinvt);
+    linalg_kernels::matrix_inv(dfm, dfinv);
+    linalg_kernels::transpose(dfinv, dfinvt);
 
     // spline evaluation; CUDA-only scratch holds the spline values pyccel keeps in args_derham
-    struphy_cuda::SplineScratch scratch;
-    struphy_cuda::get_spans(eta1, eta2, eta3, args_derham, scratch);
+    evaluation_kernels_3d::SplineScratch scratch;
+    evaluation_kernels_3d::get_spans(eta1, eta2, eta3, args_derham, scratch);
     int span1 = scratch.span1, span2 = scratch.span2, span3 = scratch.span3;
 
     // electric field: 1-form components
-    struphy_cuda::eval_1form_spline_mpi(span1, span2, span3, args_derham, scratch, e1_1, e1_2, e1_3, e_form);
+    evaluation_kernels_3d::eval_1form_spline_mpi(span1, span2, span3, args_derham, scratch, e1_1, e1_2, e1_3, e_form);
 
     // electric field: Cartesian components
-    struphy_cuda::matrix_vector(dfinvt, e_form, e_cart);
+    linalg_kernels::matrix_vector(dfinvt, e_form, e_cart);
 
     // update velocities
     for (int j = 0; j < 3; ++j) args_markers.markers(ip, 3 + j) += dt * const_factor * e_cart[j];
diff --git a/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu b/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu
index d14f81ec6..7d2a270ec 100644
--- a/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu
@@ -3,6 +3,8 @@
 #include "struphy/geometry/evaluation_kernels.cuh"
 #include "struphy/linear_algebra/linalg_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Exact rotation of the velocity about the magnetic field, as in push_vxb_analytic_kernels.push_vxb_analytic.
  *
@@ -41,21 +43,21 @@ extern "C" __global__ void push_vxb_analytic(double dt, int stage, MarkerArgs ar
     for (int j = 0; j < 3; ++j) v[j] = args_markers.markers(ip, 3 + j);
 
     // evaluate Jacobian, result in dfm
-    struphy_cuda::df(e1, e2, e3, args_domain, dfm);
+    evaluation_kernels::df(e1, e2, e3, args_domain, dfm);
 
     // metric coeffs
-    double det_df = struphy_cuda::det(dfm);
+    double det_df = linalg_kernels::det(dfm);
 
     // spline evaluation; CUDA-only scratch holds the spline values pyccel keeps in args_derham
-    struphy_cuda::SplineScratch scratch;
-    struphy_cuda::get_spans(e1, e2, e3, args_derham, scratch);
+    evaluation_kernels_3d::SplineScratch scratch;
+    evaluation_kernels_3d::get_spans(e1, e2, e3, args_derham, scratch);
     int span1 = scratch.span1, span2 = scratch.span2, span3 = scratch.span3;
 
     // magnetic field 2-form
-    struphy_cuda::eval_2form_spline_mpi(span1, span2, span3, args_derham, scratch, b2_1, b2_2, b2_3, b_form);
+    evaluation_kernels_3d::eval_2form_spline_mpi(span1, span2, span3, args_derham, scratch, b2_1, b2_2, b2_3, b_form);
 
     // magnetic field: Cartesian components
-    struphy_cuda::matrix_vector(dfm, b_form, b_cart);
+    linalg_kernels::matrix_vector(dfm, b_form, b_cart);
     for (int j = 0; j < 3; ++j) b_cart[j] = b_cart[j] / det_df;
 
     // magnetic field: magnitude
@@ -68,14 +70,14 @@ extern "C" __global__ void push_vxb_analytic(double dt, int stage, MarkerArgs ar
     for (int j = 0; j < 3; ++j) b_norm[j] = b_cart[j] / b_abs;
 
     // parallel velocity v.b_norm
-    double vpar = struphy_cuda::scalar_dot(v, b_norm);
+    double vpar = linalg_kernels::scalar_dot(v, b_norm);
 
     // first component of perpendicular velocity
-    struphy_cuda::cross(v, b_norm, vxb_norm);
-    struphy_cuda::cross(b_norm, vxb_norm, vperp);
+    linalg_kernels::cross(v, b_norm, vxb_norm);
+    linalg_kernels::cross(b_norm, vxb_norm, vperp);
 
     // second component of perpendicular velocity
-    struphy_cuda::cross(b_norm, vperp, b_normxvperp);
+    linalg_kernels::cross(b_norm, vperp, b_normxvperp);
 
     // analytic rotation
     for (int j = 0; j < 3; ++j)
diff --git a/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu b/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu
index 8fbb4a7e2..71ef3a461 100644
--- a/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu
@@ -2,6 +2,8 @@
 #include "struphy/geometry/evaluation_kernels.cuh"
 #include "struphy/linear_algebra/linalg_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Implicit (Crank-Nicolson) rotation of the velocity, as in push_vxb_implicit_kernels.push_vxb_implicit.
  *
@@ -42,21 +44,21 @@ extern "C" __global__ void push_vxb_implicit(double dt, int stage, MarkerArgs ar
     for (int j = 0; j < 3; ++j) v[j] = args_markers.markers(ip, 3 + j);
 
     // evaluate Jacobian, result in dfm
-    struphy_cuda::df(e1, e2, e3, args_domain, dfm);
+    evaluation_kernels::df(e1, e2, e3, args_domain, dfm);
 
     // metric coeffs
-    double det_df = struphy_cuda::det(dfm);
+    double det_df = linalg_kernels::det(dfm);
 
     // spline evaluation; CUDA-only scratch holds the spline values pyccel keeps in args_derham
-    struphy_cuda::SplineScratch scratch;
-    struphy_cuda::get_spans(e1, e2, e3, args_derham, scratch);
+    evaluation_kernels_3d::SplineScratch scratch;
+    evaluation_kernels_3d::get_spans(e1, e2, e3, args_derham, scratch);
     int span1 = scratch.span1, span2 = scratch.span2, span3 = scratch.span3;
 
     // magnetic field 2-form
-    struphy_cuda::eval_2form_spline_mpi(span1, span2, span3, args_derham, scratch, b2_1, b2_2, b2_3, b_form);
+    evaluation_kernels_3d::eval_2form_spline_mpi(span1, span2, span3, args_derham, scratch, b2_1, b2_2, b2_3, b_form);
 
     // magnetic field: Cartesian components
-    struphy_cuda::matrix_vector(dfm, b_form, b_cart);
+    linalg_kernels::matrix_vector(dfm, b_form, b_cart);
     for (int j = 0; j < 3; ++j) b_cart[j] = b_cart[j] / det_df;
 
     // cross-product matrix of b_cart (row-major)
@@ -73,10 +75,10 @@ extern "C" __global__ void push_vxb_implicit(double dt, int stage, MarkerArgs ar
         lhs[j] = identity[j] - dt / 2 * b_prod[j];
     }
 
-    struphy_cuda::matrix_inv(lhs, lhs_inv);
+    linalg_kernels::matrix_inv(lhs, lhs_inv);
 
-    struphy_cuda::matrix_vector(rhs, v, vec);
-    struphy_cuda::matrix_vector(lhs_inv, vec, res);
+    linalg_kernels::matrix_vector(rhs, v, vec);
+    linalg_kernels::matrix_vector(lhs_inv, vec, res);
 
     for (int j = 0; j < 3; ++j) args_markers.markers(ip, 3 + j) = res[j];
 }
diff --git a/src/struphy/pic/pushing/kernels/push_weights_with_efield_lin_va/push_weights_with_efield_lin_va_cuda.cu b/src/struphy/pic/pushing/kernels/push_weights_with_efield_lin_va/push_weights_with_efield_lin_va_cuda.cu
index a54881333..f3dfb6e00 100644
--- a/src/struphy/pic/pushing/kernels/push_weights_with_efield_lin_va/push_weights_with_efield_lin_va_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_weights_with_efield_lin_va/push_weights_with_efield_lin_va_cuda.cu
@@ -2,6 +2,8 @@
 #include "struphy/geometry/evaluation_kernels.cuh"
 #include "struphy/linear_algebra/linalg_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Weight update of the linear Vlasov-Ampere system (delta-f), as in
  * push_weights_with_efield_lin_va_kernels.push_weights_with_efield_lin_va.
@@ -48,21 +50,21 @@ extern "C" __global__ void push_weights_with_efield_lin_va(double dt, int stage,
     for (int j = 0; j < 3; ++j) v[j] = args_markers.markers(ip, 3 + j);
 
     // spline evaluation; CUDA-only scratch holds the spline values pyccel keeps in args_derham
-    struphy_cuda::SplineScratch scratch;
-    struphy_cuda::get_spans(eta1, eta2, eta3, args_derham, scratch);
+    evaluation_kernels_3d::SplineScratch scratch;
+    evaluation_kernels_3d::get_spans(eta1, eta2, eta3, args_derham, scratch);
     int span1 = scratch.span1, span2 = scratch.span2, span3 = scratch.span3;
 
     // Compute Jacobian matrix
-    struphy_cuda::df(eta1, eta2, eta3, args_domain, dfm);
+    evaluation_kernels::df(eta1, eta2, eta3, args_domain, dfm);
 
     // invert Jacobian matrix
-    struphy_cuda::matrix_inv(dfm, df_inv);
+    linalg_kernels::matrix_inv(dfm, df_inv);
 
     // compute DF^{-1} v
-    struphy_cuda::matrix_vector(df_inv, v, df_inv_v);
+    linalg_kernels::matrix_vector(df_inv, v, df_inv_v);
 
     // E-field (1-form)
-    struphy_cuda::eval_1form_spline_mpi(span1, span2, span3, args_derham, scratch, e1_1, e1_2, e1_3, e_vec);
+    evaluation_kernels_3d::eval_1form_spline_mpi(span1, span2, span3, args_derham, scratch, e1_1, e1_2, e1_3, e_vec);
 
     // w_{n+1} = w_n + kappa * dt / (2 * N * s_0 * v_th^2) * f_0 * ( DF^{-1} v_p ) \cdot ( e_{n+1} + e_n )
     double update = (df_inv_v[0] * e_vec[0] + df_inv_v[1] * e_vec[1] + df_inv_v[2] * e_vec[2]) * f0_values[ip] *
diff --git a/src/struphy/pic/pushing/kernels/reflect/reflect_cuda.cu b/src/struphy/pic/pushing/kernels/reflect/reflect_cuda.cu
index e33c2f57e..341a3197b 100644
--- a/src/struphy/pic/pushing/kernels/reflect/reflect_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/reflect/reflect_cuda.cu
@@ -1,5 +1,7 @@
 #include "struphy/pic/pushing/pusher_utilities_kernels.cuh"
 
+using namespace struphy_cuda;
+
 /**
  * Reflect velocities with the same arguments as the pyccel reflect (reflect_kernels.py).
  *
@@ -15,5 +17,5 @@ extern "C" __global__ void reflect(Array2D<double> markers, DomainArgs args_doma
                                    Array1D<long long> outside_inds, int axis) {
     long long i = (long long)blockDim.x * blockIdx.x + threadIdx.x;
     if (i >= outside_inds.shape[0]) return;
-    struphy_cuda::reflect_velocity(outside_inds(i), markers, args_domain, axis);
+    pusher_utilities_kernels::reflect_velocity(outside_inds(i), markers, args_domain, axis);
 }
diff --git a/src/struphy/pic/pushing/pusher_utilities_kernels.cuh b/src/struphy/pic/pushing/pusher_utilities_kernels.cuh
index 48b06ef19..45f3e37e2 100644
--- a/src/struphy/pic/pushing/pusher_utilities_kernels.cuh
+++ b/src/struphy/pic/pushing/pusher_utilities_kernels.cuh
@@ -2,7 +2,7 @@
 // NVRTC provides the device math function floor without a host math.h header.
 #include "struphy/geometry/evaluation_kernels.cuh"
 #include "struphy/linear_algebra/linalg_kernels.cuh"
-namespace struphy_cuda {
+namespace struphy_cuda::pusher_utilities_kernels {
 /**
  * Reflect one marker's velocity using the Pyccel boundary helper's calculation.
  *
@@ -19,13 +19,12 @@ namespace struphy_cuda {
 __device__ inline void reflect_velocity(long long ip, Array2D<double> markers,
                                        const DomainArgs& args_domain, int axis) {
     double dfm[9], dfinv[9], v[3], v_logical[3];
-    df(markers(ip, 0), markers(ip, 1), markers(ip, 2),
-       args_domain, dfm);
-    matrix_inv(dfm, dfinv);
+    evaluation_kernels::df(markers(ip, 0), markers(ip, 1), markers(ip, 2), args_domain, dfm);
+    linalg_kernels::matrix_inv(dfm, dfinv);
     for (int j = 0; j < 3; ++j) v[j] = markers(ip, 3 + j);
-    matrix_vector(dfinv, v, v_logical);
+    linalg_kernels::matrix_vector(dfinv, v, v_logical);
     v_logical[axis] *= -1.;
-    matrix_vector(dfm, v_logical, v);
+    linalg_kernels::matrix_vector(dfm, v_logical, v);
     for (int j = 0; j < 3; ++j) markers(ip, 3 + j) = v[j];
 }
 
@@ -97,4 +96,4 @@ __device__ inline void apply_kinetic_bc_marker(int ip, const MarkerArgs& args_ma
         if (reflected[axis]) reflect_velocity(ip, args_markers.markers, args_domain, axis);
     }
 }
-}
+}  // namespace struphy_cuda::pusher_utilities_kernels
diff --git a/src/struphy/pic/tests/test_device_helpers.py b/src/struphy/pic/tests/test_device_helpers.py
index b771d8efc..d0bfaaaa7 100644
--- a/src/struphy/pic/tests/test_device_helpers.py
+++ b/src/struphy/pic/tests/test_device_helpers.py
@@ -25,11 +25,11 @@ def test_spline_helpers(degree):
 extern "C" __global__ void evaluate(const double* t, int nt, int p, const double* x, double* out, int n) {
     int i=blockDim.x*blockIdx.x+threadIdx.x;
     if(i>=n) return;
-    int span=struphy_cuda::find_span(t,nt,p,x[i]);
+    int span=struphy_cuda::bsplines_kernels::find_span(t,nt,p,x[i]);
     double* row=out+i*(3*p+3);
     row[0]=span;
-    struphy_cuda::basis_funs(t,p,x[i],span,row+1);
-    struphy_cuda::b_d_splines_slim(t,p,x[i],span,row+p+2,row+2*p+3);
+    struphy_cuda::bsplines_kernels::basis_funs(t,p,x[i],span,row+1);
+    struphy_cuda::bsplines_kernels::b_d_splines_slim(t,p,x[i],span,row+p+2,row+2*p+3);
 }
 """
     out = cp.empty((len(points), 3 * degree + 3))
@@ -59,8 +59,8 @@ def test_matrix_helpers():
 extern "C" __global__ void evaluate(const double* a, const double* v, double* inv, double* out, int n) {
     int i=blockDim.x*blockIdx.x+threadIdx.x;
     if(i>=n) return;
-    struphy_cuda::matrix_inv(a+9*i,inv+9*i);
-    struphy_cuda::matrix_vector(inv+9*i,v+3*i,out+3*i);
+    struphy_cuda::linalg_kernels::matrix_inv(a+9*i,inv+9*i);
+    struphy_cuda::linalg_kernels::matrix_vector(inv+9*i,v+3*i,out+3*i);
 }
 """
     inv, out = cp.empty(matrices.shape), cp.empty(vectors.shape)
@@ -86,7 +86,7 @@ def test_boundary_helpers(bc, newton):
 #include "struphy/pic/pushing/pusher_utilities_kernels.cuh"
 extern "C" __global__ void evaluate(MarkerArgs m, DomainArgs d, int newton) {
     int i=blockDim.x*blockIdx.x+threadIdx.x;
-    if(i<m.n_markers) struphy_cuda::apply_kinetic_bc_marker(i,m,d,newton);
+    if(i<m.n_markers) struphy_cuda::pusher_utilities_kernels::apply_kinetic_bc_marker(i,m,d,newton);
 }
 """
     results = []
@@ -114,7 +114,7 @@ def test_span_with_cunumpy_wrapper(degree):
     kernel = device_function_kernel(
         '#include "struphy/bsplines/bsplines_kernels.cuh"\n'
         "__device__ int span_at(const double* t, int nt, int p, double x) {"
-        "return struphy_cuda::find_span(t, nt, p, x);}",
+        "return struphy_cuda::bsplines_kernels::find_span(t, nt, p, x);}",
         "int span_at(const double* t, int nt, int p, double x)",
         **CUDA_OPTIONS,
     )
@@ -141,14 +141,14 @@ def test_cuboid_helpers():
     source = r"""
 #include "struphy/geometry/evaluation_kernels.cuh"
 extern "C" __global__ void evaluate(const double* eta, const double* params, double* out) {
-    struphy_cuda::cuboid(eta[0], eta[1], eta[2], params[0], params[1], params[2],
+    struphy_cuda::cuboid_kernels::cuboid(eta[0], eta[1], eta[2], params[0], params[1], params[2],
                         params[3], params[4], params[5], out);
-    struphy_cuda::cuboid_df(params[0], params[1], params[2],
+    struphy_cuda::cuboid_kernels::cuboid_df(params[0], params[1], params[2],
                            params[3], params[4], params[5], out + 3);
     DomainArgs args = {};
     args.kind_map = 10;
     args.params = const_cast<double*>(params);
-    struphy_cuda::df(eta[0], eta[1], eta[2], args, out + 12);
+    struphy_cuda::evaluation_kernels::df(eta[0], eta[1], eta[2], args, out + 12);
 }
 """
     eta = np.array([0.2, 0.4, 0.8])
@@ -177,8 +177,8 @@ def test_get_spans_helper():
 extern "C" __global__ void evaluate(const double* eta, DerhamArgs args_derham, double* out, int n) {
     int ip = blockDim.x * blockIdx.x + threadIdx.x;
     if (ip >= n) return;
-    struphy_cuda::SplineScratch scratch;
-    struphy_cuda::get_spans(eta[3*ip], eta[3*ip+1], eta[3*ip+2], args_derham, scratch);
+    struphy_cuda::evaluation_kernels_3d::SplineScratch scratch;
+    struphy_cuda::evaluation_kernels_3d::get_spans(eta[3*ip], eta[3*ip+1], eta[3*ip+2], args_derham, scratch);
     const int spans[3] = {scratch.span1, scratch.span2, scratch.span3};
     const double* bn[3] = {scratch.bn1, scratch.bn2, scratch.bn3};
     const double* bd[3] = {scratch.bd1, scratch.bd2, scratch.bd3};
@@ -257,20 +257,20 @@ def test_spline_mapping_helpers(name):
 __device__ double metric_entry(double eta1, double eta2, double eta3, int kind_coeff, int entry,
                                bool avoid_round_off, const DomainArgs& args) {
     double tmp1[9], tmp2[9], tmp3[9], out[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
-    if (kind_coeff == 0) struphy_cuda::f(eta1, eta2, eta3, args, out);
-    else if (kind_coeff == 1) struphy_cuda::df(eta1, eta2, eta3, args, out);
-    else if (kind_coeff == 2) out[0] = struphy_cuda::det_df(eta1, eta2, eta3, args, tmp1);
-    else if (kind_coeff == 3) struphy_cuda::df_inv(eta1, eta2, eta3, args, tmp1, avoid_round_off, out);
-    else if (kind_coeff == 4) struphy_cuda::g(eta1, eta2, eta3, args, tmp1, tmp2, avoid_round_off, out);
-    else struphy_cuda::g_inv(eta1, eta2, eta3, args, tmp1, tmp2, tmp3, avoid_round_off, out);
+    if (kind_coeff == 0) struphy_cuda::evaluation_kernels::f(eta1, eta2, eta3, args, out);
+    else if (kind_coeff == 1) struphy_cuda::evaluation_kernels::df(eta1, eta2, eta3, args, out);
+    else if (kind_coeff == 2) out[0] = struphy_cuda::evaluation_kernels::det_df(eta1, eta2, eta3, args, tmp1);
+    else if (kind_coeff == 3) struphy_cuda::evaluation_kernels::df_inv(eta1, eta2, eta3, args, tmp1, avoid_round_off, out);
+    else if (kind_coeff == 4) struphy_cuda::evaluation_kernels::g(eta1, eta2, eta3, args, tmp1, tmp2, avoid_round_off, out);
+    else struphy_cuda::evaluation_kernels::g_inv(eta1, eta2, eta3, args, tmp1, tmp2, tmp3, avoid_round_off, out);
     return out[entry];
 }
 __device__ double transform_entry(double a0, double a1, double a2, double eta1, double eta2, double eta3,
                                   int kind_transform, int kind_fun, int entry, const DomainArgs& args_domain) {
     double a[3] = {a0, a1, a2}, out[3] = {0, 0, 0};
-    if (kind_transform == 0) struphy_cuda::pull(a, eta1, eta2, eta3, kind_fun, args_domain, out);
-    else if (kind_transform == 1) struphy_cuda::push(a, eta1, eta2, eta3, kind_fun, args_domain, out);
-    else struphy_cuda::tran(a, eta1, eta2, eta3, kind_fun, args_domain, out);
+    if (kind_transform == 0) struphy_cuda::transform_kernels::pull(a, eta1, eta2, eta3, kind_fun, args_domain, out);
+    else if (kind_transform == 1) struphy_cuda::transform_kernels::push(a, eta1, eta2, eta3, kind_fun, args_domain, out);
+    else struphy_cuda::transform_kernels::tran(a, eta1, eta2, eta3, kind_fun, args_domain, out);
     return out[entry];
 }
 """
@@ -403,9 +403,10 @@ def test_transform_helpers(domain_index):
     int ip = blockDim.x * blockIdx.x + threadIdx.x;
     if (ip >= n) return;
     const double* fill = fills + 9 * ip;
-    struphy_cuda::m_v_fill_b_v1_symm(args_derham, eta[3 * ip], eta[3 * ip + 1], eta[3 * ip + 2], mat11, mat12, mat13,
-                                     mat22, mat23, mat33, fill[0], fill[1], fill[2], fill[3], fill[4], fill[5], vec1,
-                                     vec2, vec3, fill[6], fill[7], fill[8]);
+    struphy_cuda::particle_to_mat_kernels::m_v_fill_b_v1_symm(args_derham, eta[3 * ip], eta[3 * ip + 1],
+                                                              eta[3 * ip + 2], mat11, mat12, mat13, mat22, mat23,
+                                                              mat33, fill[0], fill[1], fill[2], fill[3], fill[4],
+                                                              fill[5], vec1, vec2, vec3, fill[6], fill[7], fill[8]);
 }
 """