Fuse the pusher preamble into one CUDA kernel

src/struphy/kernel_arguments/pusher_args.cuh

@@ -12,7 +12,7 @@ struct MarkerArgs {
     int vdim;
     int weight_idx;
     int first_diagnostics_idx;
-    int first_init_idx;
+    int first_pusher_idx;
     int first_shift_idx;
     int residual_idx;
     int first_free_idx;

src/struphy/kernel_arguments/pusher_args_cuda.py

@@ -39,7 +39,7 @@ class CudaMarkerArguments(CudaStructArguments):
         ("vdim", "int"),
         ("weight_idx", "int"),
         ("first_diagnostics_idx", "int"),
-        ("first_init_idx", "int"),
+        ("first_pusher_idx", "int"),
         ("first_shift_idx", "int"),
         ("residual_idx", "int"),
         ("first_free_idx", "int"),
@@ -71,7 +71,7 @@ def __init__(
         self.weight_idx = weight_idx
         self.n_markers = markers.shape[0]
         self.first_diagnostics_idx = first_diagnostics_idx
-        self.first_init_idx = first_pusher_idx
+        self.first_pusher_idx = first_pusher_idx
         self.first_shift_idx = first_shift_idx
         self.residual_idx = residual_idx
         self.first_free_idx = first_free_idx

src/struphy/kernel_arguments/pusher_args_kernels.py

@@ -71,7 +71,7 @@ def __init__(
 
         # useful indices
         self.first_diagnostics_idx = first_diagnostics_idx
-        self.first_init_idx = first_pusher_idx
+        self.first_pusher_idx = first_pusher_idx
         self.first_shift_idx = first_shift_idx  # starting idx for eta-shifts due to boundary conditions
         self.residual_idx = residual_idx  # residual in iterative solvers
         self.first_free_idx = first_free_idx  # index after which auxiliary saving is possible

src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB/cc_lin_mhd_5d_gradB_kernels.py

@@ -99,7 +99,7 @@ def cc_lin_mhd_5d_gradB(
 
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     mu_idx = args_markers.mu_idx
 
     # allocate for magnetic field evaluation
@@ -130,7 +130,7 @@ def cc_lin_mhd_5d_gradB(
             continue
 
         # if particle is refilled
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         # marker positions

src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB_dg/cc_lin_mhd_5d_gradB_dg_kernels.py

@@ -71,7 +71,7 @@ def cc_lin_mhd_5d_gradB_dg(
 
     markers = args_markers.markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
 
     # allocate for magnetic field evaluation
     eta_diff = empty(3, dtype=float)
@@ -108,10 +108,10 @@ def cc_lin_mhd_5d_gradB_dg(
             continue
 
         # marker positions, mid point
-        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_init_idx : first_init_idx + 3]) / 2.0
+        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_pusher_idx : first_pusher_idx + 3]) / 2.0
         eta_mid[:] = mod(eta_mid[:], 1.0)
 
-        eta_diff[:] = markers[ip, 0:3] - markers[ip, first_init_idx : first_init_idx + 3]
+        eta_diff[:] = markers[ip, 0:3] - markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         # marker weight and velocity
         weight = markers[ip, 5]

src/struphy/pic/diagnostics/kernels/eval_gradB_ediff/eval_gradB_ediff_kernels.py

@@ -42,7 +42,7 @@ def eval_gradB_ediff(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
@@ -51,10 +51,10 @@ def eval_gradB_ediff(
             continue
 
         # marker positions, mid point
-        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_init_idx : first_init_idx + 3]) / 2.0
+        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_pusher_idx : first_pusher_idx + 3]) / 2.0
         eta_mid[:] = mod(eta_mid[:], 1.0)
 
-        eta_diff = markers[ip, 0:3] - markers[ip, first_init_idx : first_init_idx + 3]
+        eta_diff = markers[ip, 0:3] - markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         # marker weight and velocity
         weight = markers[ip, 5]

src/struphy/pic/pushing/kernels/bstar_2form/bstar_2form_kernels.py

@@ -50,7 +50,7 @@ def bstar_2form(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -59,13 +59,13 @@ def bstar_2form(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v = alpha[3] * v_k + (1.0 - alpha[3]) * v_n
 
         # spline evaluation

src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_kernels.py

@@ -45,7 +45,7 @@ def bstar_parallel_3form(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -54,13 +54,13 @@ def bstar_parallel_3form(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v = alpha[3] * v_k + (1.0 - alpha[3]) * v_n
 
         # evaluate Jacobian, result in dfm

src/struphy/pic/pushing/kernels/driftkinetic_hamiltonian/driftkinetic_hamiltonian_kernels.py

@@ -44,7 +44,7 @@ def driftkinetic_hamiltonian(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -53,13 +53,13 @@ def driftkinetic_hamiltonian(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v = alpha[3] * v_k + (1.0 - alpha[3]) * v_n
 
         mu = markers[ip, mu_idx]

src/struphy/pic/pushing/kernels/grad_driftkinetic_hamiltonian/grad_driftkinetic_hamiltonian_kernels.py

@@ -58,7 +58,7 @@ def grad_driftkinetic_hamiltonian(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -67,7 +67,7 @@ def grad_driftkinetic_hamiltonian(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)

src/struphy/pic/pushing/kernels/push_deterministic_diffusion_stage/push_deterministic_diffusion_stage_kernels.py

@@ -49,7 +49,7 @@ def push_deterministic_diffusion_stage(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -121,7 +121,7 @@ def push_deterministic_diffusion_stage(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * a[stage] * k
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )

src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu

@@ -14,7 +14,7 @@
  * @param c Butcher nodes c (n_stages entries; unused, as in the pyccel kernel).
  * @param n_stages Number of stages; a raw device pointer carries no length.
  *
- * Holes (first_init_idx == -1) and boundary particles (last column == -2) are skipped.
+ * Holes (first_pusher_idx == -1) and boundary particles (last column == -2) are skipped.
  * The kinetic boundary conditions are applied per marker after the update.
  */
 extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_markers, DomainArgs args_domain,
@@ -24,12 +24,12 @@ extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_
     int n_markers = args_markers.n_markers;
     if (ip >= n_markers) return;
 
-    int first_init_idx = args_markers.first_init_idx;
+    int first_pusher_idx = args_markers.first_pusher_idx;
     int first_free_idx = args_markers.first_free_idx;
     double last = stage == n_stages - 1 ? 1. : 0.;
 
     // check if marker is a hole or a boundary particle
-    if (args_markers.markers(ip, first_init_idx) == -1. || args_markers.markers(ip, args_markers.markers.shape[1] - 1) == -2.)
+    if (args_markers.markers(ip, first_pusher_idx) == -1. || args_markers.markers(ip, args_markers.markers.shape[1] - 1) == -2.)
         return;
 
     double e1 = args_markers.markers(ip, 0);
@@ -51,7 +51,7 @@ extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_
         // accumulation for last stage
         args_markers.markers(ip, first_free_idx + j) += dt * b[stage] * k[j];
         // update positions for intermediate stages or last stage
-        args_markers.markers(ip, j) = args_markers.markers(ip, first_init_idx + j) + dt * a[stage] * k[j] +
+        args_markers.markers(ip, j) = args_markers.markers(ip, first_pusher_idx + j) + dt * a[stage] * k[j] +
                                       last * args_markers.markers(ip, first_free_idx + j);
     }
 

src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_kernels.py

@@ -43,7 +43,7 @@ def push_eta_stage(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     if stage == n_stages - 1:
@@ -57,7 +57,7 @@ def push_eta_stage(
     # fmt: on
     for ip in range(n_markers):
         # check if marker is a hole or a boundary particle
-        if markers[ip, first_init_idx] == -1.0 or markers[ip, -1] == -2.0:
+        if markers[ip, first_pusher_idx] == -1.0 or markers[ip, -1] == -2.0:
             continue
 
         e1 = markers[ip, 0]
@@ -85,7 +85,7 @@ def push_eta_stage(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * a[stage] * k
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )

src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order/push_gc_Bstar_discrete_gradient_1st_order_kernels.py

@@ -82,25 +82,25 @@ def push_gc_Bstar_discrete_gradient_1st_order(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta_mid[:] = (eta_k + eta_n) / 2.0
         eta_mid[:] = mod(eta_mid, 1.0)
         eta_diff[:] = eta_k - eta_n
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v_mid = (v_k + v_n) / 2.0
         v_diff = v_k - v_n
 

src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order_newton/push_gc_Bstar_discrete_gradient_1st_order_newton_kernels.py

@@ -143,23 +143,23 @@ def push_gc_Bstar_discrete_gradient_1st_order_newton(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
         eta_k_shifted[:] = eta_k + markers[ip, first_shift_idx : first_shift_idx + 3]
         eta_diff[:] = eta_k_shifted - eta_n
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v_diff = v_k - v_n
 
         mu = markers[ip, mu_idx]

src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_2nd_order/push_gc_Bstar_discrete_gradient_2nd_order_kernels.py

@@ -108,25 +108,25 @@ def push_gc_Bstar_discrete_gradient_2nd_order(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta_mid[:] = (eta_k + eta_n) / 2.0
         eta_mid[:] = mod(eta_mid, 1.0)
         eta_diff[:] = eta_k - eta_n
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v_mid = (v_k + v_n) / 2.0
         v_diff = v_k - v_n
 

src/struphy/pic/pushing/kernels/push_gc_Bstar_explicit_multistage/push_gc_Bstar_explicit_multistage_kernels.py

@@ -80,7 +80,7 @@ def push_gc_Bstar_explicit_multistage(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -93,7 +93,7 @@ def push_gc_Bstar_explicit_multistage(
 
     for ip in range(n_markers):
         # check if marker is a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         # if stage == 0.:
@@ -207,11 +207,13 @@ def push_gc_Bstar_explicit_multistage(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * a[stage] * k
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
-        markers[ip, 3] = markers[ip, first_init_idx + 3] + dt * a[stage] * k_v + last * markers[ip, first_free_idx + 3]
+        markers[ip, 3] = (
+            markers[ip, first_pusher_idx + 3] + dt * a[stage] * k_v + last * markers[ip, first_free_idx + 3]
+        )
 
         # apply kinetic boundary conditions
         pusher_utilities_kernels.apply_kinetic_bc_marker(ip, args_markers, args_domain, False)