Total magnetic field in B* = B0 + B~ + eps v_par curl b0 (PushGuidingCenterParallel, DriftKineticMaxwellCoupling)

src/struphy/models/drift_kinetic_maxwell.py

@@ -185,11 +185,12 @@ def doc_pde(cls):
         .. math::
 
             \mathbf E^* = \mathbf E - \varepsilon \mu \nabla |\mathbf B|\,,\qquad
-            \mathbf B^* = \mathbf B_0 + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad
-            B^*_\parallel = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B}) + \varepsilon v_\parallel
-            (\nabla \times \mathbf b_0) \cdot \mathbf b_0\,,
+            \mathbf B^* = \mathbf B_0 + \tilde{\mathbf B} + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad
+            B^*_\parallel = \mathbf b_0 \cdot \mathbf B^*\,,
 
-        where :math:`|\mathbf B| = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B})`.
+        where :math:`|\mathbf B| = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B})`: the total magnetic field in
+        the first term of :math:`\mathbf B^*` (parallel streaming along the perturbed field lines), the curvature
+        term with the equilibrium unit vector :math:`\mathbf b_0`.
 
         Ampère's law, with the displacement current and the (linearized) polarization current of the guiding centers:
 

src/struphy/propagators/drift_kinetic_maxwell_coupling.py

@@ -53,9 +53,11 @@ class DriftKineticMaxwellCoupling(Propagator):
 
     .. math::
 
-        \mathbf B^* = \mathbf B_0 + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad
-        B^*_\parallel = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B}) + \varepsilon v_\parallel
-        (\nabla \times \mathbf b_0) \cdot \mathbf b_0\,.
+        \mathbf B^* = \mathbf B_0 + \tilde{\mathbf B} + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad
+        B^*_\parallel = \mathbf b_0 \cdot \mathbf B^*\,,
+
+    with the total magnetic field in the first term (parallel streaming along the perturbed field lines) and the
+    curvature term of the equilibrium.
 
     The current is the guiding-center current :math:`\iint f^\textnormal{vol} \dot{\mathbf X}` without the
     magnetic drifts (which move the positions): parallel streaming, including the curvature drift in
@@ -371,8 +373,10 @@ def allocate(self):
 
         self._info = self.options.solver_params.info
 
-        # magnetic background: B* (2-form, equilibrium) and B*_parallel with the perturbation b (frozen in a step)
-        self._b2 = self.projected_equil.b2
+        # magnetic field: B* (2-form) and B*_parallel with the total field B0 + b (b frozen in a step)
+        self._b2_eq = self.projected_equil.b2
+        self._b2 = self._b2_eq.space.zeros()
+        self._E2T = self.derham.extraction_ops["2"].transpose()
         self._curl_unit_b2 = self.projected_equil.curl_unit_b2
         self._unit_b1 = self.projected_equil.unit_b1
         self._absB0 = self.projected_equil.absB0
@@ -382,7 +386,7 @@ def allocate(self):
         self._PB = getattr(self.basis_ops, "PB")
         self._B_dot_b = self.derham.V0.zeros()
         self._B_dot_b_pol = self._PB.codomain.zeros()
-        self._update_b_dot_b()
+        self._update_b_star()
 
         args_fields = (
             self._b2[0]._data,
@@ -530,9 +534,15 @@ def allocate(self):
             pushes_eta=False,
         )
 
-    def _update_b_dot_b(self):
-        """Parallel magnetic field b0 . (B0 + b) in the tensor-product space V0 (used in B*_parallel)."""
-        B_dot_b_pol = self._PB.dot(self.variables.b.spline.vector, out=self._B_dot_b_pol)
+    def _update_b_star(self):
+        """Total magnetic field B0 + b (2-form, in B*) and parallel field b0 . (B0 + b) (0-form, in B*_parallel) in the
+        tensor-product spaces."""
+        b = self.variables.b.spline.vector
+        b2 = self._E2T.dot(b, out=self._b2)
+        b2 += self._b2_eq
+        b2.update_ghost_regions()
+
+        B_dot_b_pol = self._PB.dot(b, out=self._B_dot_b_pol)
         B_dot_b = self._E0T.dot(B_dot_b_pol, out=self._B_dot_b)
         B_dot_b += self._absB0
         B_dot_b.update_ghost_regions()
@@ -627,7 +637,7 @@ def _mid_point_field(self, dt: float):
         bn = self.variables.b.spline.vector
 
         # B*_parallel with b^n, then A_par, S (or their fillings per marker) and j
-        self._update_b_dot_b()
+        self._update_b_star()
         self._accumulate()
 
         # rhs = M_E e^n + dt/2 (c1 C^T M2 b^n - c2 j)

src/struphy/propagators/push_guiding_center_bx_estar.py

@@ -92,7 +92,7 @@ def __init__(
             Electric potential in ``"H1"`` space contributing to :math:`E^*`.
             If ``None``, an empty ``FEECVariable(space="H1")`` is created internally.
         b_tilde : FEECVariable, default=None
-            Magnetic perturbation in ``"Hcurl"`` space contributing to :math:`B^*`.
+            Magnetic perturbation in ``"Hdiv"`` space contributing to :math:`B^*`.
             If ``None``, the perturbation is ignored.
         e_field : FEECVariable, default=None
             Electric field in ``"Hcurl"`` space, used in :math:`E^* = E - \varepsilon \mu \nabla |B|` instead of

src/struphy/propagators/push_guiding_center_parallel.py

@@ -108,7 +108,7 @@ def __init__(
             Electric potential in ``"H1"`` space contributing to :math:`E^*`.
             If ``None``, an empty ``FEECVariable(space="H1")`` is created internally.
         b_tilde : FEECVariable, default=None
-            Magnetic perturbation in ``"Hcurl"`` space contributing to :math:`B^*`.
+            Magnetic perturbation in ``"Hdiv"`` space contributing to :math:`B^*`.
             If ``None``, the perturbation is ignored.
         e_field : FEECVariable, default=None
             Electric field in ``"Hcurl"`` space, used in :math:`E^* = E - \varepsilon \mu \nabla |B|` instead of
@@ -209,11 +209,15 @@ def allocate(self):
         # PolarVector --> Stencil-/BlockVector (identity maps in case of no polar splines)
         self._E0T = self.derham.extraction_ops["0"].transpose()
         self._E1T = self.derham.extraction_ops["1"].transpose()
+        self._E2T = self.derham.extraction_ops["2"].transpose()
 
-        # magnetic perturbation
+        # magnetic perturbation (in B* = B0 + b_tilde + epsilon v_par curl b0 and in B*_parallel and |B|)
         if self.b_tilde is not None:
             self._B_dot_b = self.derham.V0.zeros()
             self._grad_b_full = self.derham.V1.zeros()
+            self._b2_eq = b2
+            self._b2_full = b2.space.zeros()
+            b2 = self._b2_full
 
             self._PB = getattr(self.basis_ops, "PB")
             self._B_dot_b_pol = self._PB.codomain.zeros()
@@ -536,7 +540,8 @@ def _update_phi(self):
         phi.update_ghost_regions()
 
     def _update_b_tilde(self):
-        """Parallel magnetic perturbation and its gradient (plus equilibrium) in the tensor-product spaces."""
+        """Total magnetic field (2-form), parallel magnetic field and its gradient (equilibrium plus perturbation) in the
+        tensor-product spaces."""
         B_dot_b_pol = self._PB.dot(self.b_tilde.spline.vector, out=self._B_dot_b_pol)
         B_dot_b = self._E0T.dot(B_dot_b_pol, out=self._B_dot_b)
         B_dot_b.update_ghost_regions()
@@ -548,6 +553,11 @@ def _update_b_tilde(self):
         grad_b_full += self._gradB1
         B_dot_b += self._absB0
 
+        # total magnetic field B0 + b_tilde (2-form) in B*
+        b2_full = self._E2T.dot(self.b_tilde.spline.vector, out=self._b2_full)
+        b2_full += self._b2_eq
+        b2_full.update_ghost_regions()
+
     @profile
     def __call__(self, dt):
         # electric potential

src/struphy/propagators/tests/test_drift_kinetic_maxwell_coupling.py

@@ -286,6 +286,9 @@ def test_parallel_force_applied_once(algo, evaluate_e_field):
     model = _cold_model()
     prop = model.propagators.push_para
     prop.options = prop.Options(algo=algo, evaluate_e_field=evaluate_e_field)
+    # re-allocate with the new options: drop the kernels registered by the first allocation
+    prop._init_kernels = ()
+    prop._eval_kernels = ()
     prop.allocate()
     particles = model.kinetic_ions.var.particles
     v_0 = _all_vpar(particles)
@@ -295,6 +298,9 @@ def test_parallel_force_applied_once(algo, evaluate_e_field):
     # the same push with the electric field changes v_par (by about dt / epsilon * b* . E)
     model.propagators.push_para.e_field = model.em_fields.e_field
     prop.options = prop.Options(algo="explicit")
+    # re-allocate with the new options: drop the kernels registered by the first allocation
+    prop._init_kernels = ()
+    prop._eval_kernels = ()
     prop.allocate()
     prop(0.1)
     assert xp.max(xp.abs(_all_vpar(particles) - v_0)) > 1e-2

src/struphy/propagators/tests/test_push_guiding_center_full_b.py

@@ -0,0 +1,106 @@
+import logging
+
+import cunumpy as xp
+import pytest
+
+from struphy import set_logging_level
+
+set_logging_level(logging.WARNING)
+
+LX, LY, LZ = 2.0, 2.0, 4.0
+DB = 0.1  # uniform perturbation B~_x on B0 = e_z
+
+
+def _model(**coupling_options):
+    """DriftKineticMaxwell in a box with B0 = e_z and the uniform perturbation B~ = DB e_x, cold guiding centers
+    (mu ~ 0) streaming with v_par ~ 1, no electric field."""
+    from struphy import (
+        DerhamOptions,
+        LoadingParameters,
+        Simulation,
+        domains,
+        equils,
+        grids,
+        maxwellians,
+        models,
+        perturbations,
+    )
+
+    vth = 1e-6
+    model = models.DriftKineticMaxwell(epsilon=0.5, lam=2.0, kappa=2.0)
+    model.kinetic_ions.set_markers(
+        loading_params=LoadingParameters(ppc=5, seed=1234, moments=(1.0, 0.0, vth, vth), B0=1.0)
+    )
+    model.kinetic_ions.var.add_background(
+        maxwellians.GyroMaxwellian2D(
+            n=(1.0, None), u_para=(1.0, None), vth_para=(vth, None), vth_perp=(vth, None), B0=1.0
+        )
+    )
+    model.em_fields.b_field.add_perturbation(
+        perturbations.ModesCos(ls=(0,), ms=(0,), ns=(0,), amps=(DB,), given_in_basis="physical", comp=0)
+    )
+    prop = model.propagators.coupling
+    prop.options = prop.Options(**coupling_options)
+    Simulation(
+        model=model,
+        domain=domains.Cuboid(r1=LX, r2=LY, r3=LZ),
+        equil=equils.HomogenSlab(B0z=1.0),
+        grid=grids.TensorProductGrid(num_elements=(4, 4, 4)),
+        derham_opts=DerhamOptions(degree=(2, 2, 2)),
+    ).allocate()
+    return model
+
+
+def _ids_and_positions(particles):
+    """Marker ids and logical positions of all valid markers on all processes, sorted by id."""
+    from maybempi import MPI
+
+    local = xp.vstack((particles.marker_ids, particles.markers[:3, particles.valid_mks]))
+    comm = MPI.COMM_WORLD
+    if comm.Get_size() > 1:
+        local = xp.concatenate(comm.allgather(local), axis=1)
+    return local[:, xp.argsort(local[0])]
+
+
+@pytest.mark.parametrize("algo", ["discrete_gradient_1st_order", "explicit"])
+def test_push_parallel_along_full_b(algo):
+    """PushGuidingCenterParallel with b_tilde streams along B* = B0 + B~ (+ curvature term): with B~ = DB e_x on
+    B0 = e_z the guiding centers move in x by DB times their displacement in z (before, B* had B0 only)."""
+    model = _model()
+    prop = model.propagators.push_para
+    prop.options = prop.Options(algo=algo)
+    # re-allocate with the new options: drop the kernels registered by the first allocation
+    prop._init_kernels = ()
+    prop._eval_kernels = ()
+    prop.allocate()
+
+    particles = model.kinetic_ions.var.particles
+    before = _ids_and_positions(particles)
+    assert len(xp.unique(before[0])) == before.shape[1]
+    prop(0.05)
+    after = _ids_and_positions(particles)
+    assert xp.array_equal(before[0], after[0])
+
+    # displacements (periodic, unwrapped) in physical units
+    d_eta = (after[1:] - before[1:] + 0.5) % 1.0 - 0.5
+    dx, dz = LX * d_eta[0], LZ * d_eta[2]
+    assert xp.min(xp.abs(dz)) > 1e-3
+    assert xp.max(xp.abs(dx - DB * dz)) < 1e-10 * xp.max(xp.abs(dz))
+    assert xp.max(xp.abs(LY * d_eta[1])) < 1e-12
+
+
+def test_coupling_parallel_direction_full_b():
+    """The parallel response of DriftKineticMaxwellCoupling uses the same B* = B0 + B~ (+ curvature term): the
+    stored a_p = B*2 / B*3_par (logical components, matrix-free) points along the total field."""
+    model = _model(matrix_free=True)
+    prop = model.propagators.coupling
+    prop._update_b_star()
+    prop._accumulate()
+
+    particles = model.kinetic_ions.var.particles
+    column = particles.args_markers.first_free_idx
+    a = particles.markers[column : column + 3, particles.valid_mks]
+    # physical B* = (DB, 0, 1), |b0 . B*| = 1: logical vector components DF^{-1} B* = (DB / LX, 0, 1 / LZ)
+    assert xp.allclose(a[0], DB / LX, rtol=1e-10)
+    assert xp.allclose(a[1], 0.0, atol=1e-12)
+    assert xp.allclose(a[2], 1.0 / LZ, rtol=1e-10)