Particles on domains with a pole (markers through the axis, polar splines in particle propagators)

[spossann]

Stacked on #759 (polar-profiling). It retargets to devel once #758/#759 are merged.

Summary

Particle markers can now run on domains with a pole at η1 = 0 (no hole around the axis), with polar splines for the FEEC parts. Examples are IGAPolarCylinder and IGAPolarTorus (HollowCylinder/HollowTorus with a1=0 also work in kernels). Before this PR, Vlasov and GuidingCenter failed on such domains: several propagators did not handle PolarVector, no marker boundary condition handled the axis, and the kernels divided by det DF = 0.

This prepares the Vlasov/GuidingCenter profiling examples (next PR, stacked on this one).

Marker push through the axis

• New update_position_stage (pic/pushing/pusher_utilities_kernels.py): the stage update of the explicit Runge–Kutta pushers (PushEta, explicit GC bxE* and B*) is done in ξ = η1·(cos 2πη2, sin 2πη2), keeping η3.
  • In these coordinates the mapping and the marker velocity dξ/dt = (∂ξ/∂η)·dη/dt are regular at the pole, and markers move through the axis.
  • Straight lines through the axis of a cylinder are exact (tested to 1e-10).
  • Without a pole the update is unchanged.
• Axis map in apply_kinetic_bc_marker (Pyccel and CUDA) and Particles.apply_kinetic_bc:
  • A marker with η1 < 0 is mapped to (−η1, η2 + ½), with unchanged velocities.
  • η1 is kept ≥ 1e-10, also after loading (Sobol loading starts at η = 0).
  • The boundary condition of axis 0 then applies at η1 = 1 only.
• New field MarkerArguments.pole, also in CudaMarkerArguments and pusher_args.cuh.
• Not supported yet: pole domains on the CuPy backend, and poc/pol_period ≠ 1. Both raise NotImplementedError.

Polar splines in particle propagators

• PushVxB: the projected equilibrium b2 (already tensor-product) was extracted a second time, which failed with polar splines. Only the b2 variable is extracted now.
• PushGuidingCenterBxEstar/Parallel: φ, the electric field and b̃ are passed to the kernels as tensor-product coefficients (PolarVector has no _data). The gradients use derham.grad, boundary operators included, followed by E^T. Without polar splines nothing changes.

Other fixes found along the way

• Reflecting boundary: only Cartesian velocities (vdim = 3) are reflected. For Particles5D the columns 3:6 (v∥, μ, weight) were overwritten.
• CanonicalMaxwellian2D: removed an unused read of domain.params["a1"], which raised a KeyError on IGAPolarTorus.
• AdhocTorus, AdhocTorusQPsi: the ψ derivatives were 0/0 (NaN) on the magnetic axis. Off the axis they are unchanged (4e-16); on the axis they take the correct limits.
  • Note: AdhocTorusQPsi's interpolated ψ has ψ′(0) ≈ −6e-5 ≠ 0, so its field is not exactly regular on the axis. That is pre-existing and not changed here.
• GuidingCenter: en_fB summed over holes and ghosts instead of the valid markers. Each removed marker added +1 to en_tot.

Disc loading under MPI (pre-existing bug, not pole-related)

With spatial="disc" on more than one process, markers were lost in the first push or the run hung:

• Disc transform after the selection: pseudo_random loading selected the uniform draws on each process's domain and applied η1 → √η1 afterwards. That moved many markers to another process's domain.
  • Without sorting boxes they were not sent there, so the first push evaluated the fields outside the local domain.
  • Example: GuidingCenter on 4 processes removed 15% of the markers in the first step.
  • The transform is now applied before the domain check.
• Sort after loading: markers are sent to their process even without sorting boxes (not for restart loading). Sobol loading draws on the whole domain on each process, so it needed this too.
• Marker array size: the number of markers per process was estimated with the uniform measure. With the disc measure 2η1, processes at larger η1 overflowed the marker array: a ValueError on some processes and a hang on the others.
• New test: test_draw_disc_on_own_process (pseudo_random and sobol_standard) in test_draw_parallel.py.

Known limitations

• Discrete-gradient GC pushers: they only get the axis map. Their mid-points in η are not regular near the axis.
  • In the test (ε = 1), relative energy errors near the axis reach 1e-3 to 5e-2.
  • With the explicit RK4 pushers, markers that cross the axis conserve their energy to about 1e-6.
• Buffer size: with few markers per process, the discrete-gradient pushers need a larger bufsize on 4 ranks (load-imbalance abort in the box sort). This is also the case on a torus with a hole, so it is not pole-related.

Tests

New file src/struphy/pic/tests/test_pole.py, which passes serially and on 4 MPI ranks:

• the axis map in the kernel and in Python, including a marker with η1 > 1;
• reflection for vdim = 2 and 3;
• PushEta straight lines through the axis (forward Euler, Heun2, RK4) against the exact Cartesian end points;
• Vlasov on IGAPolarCylinder with ScrewPinch and polar splines: kinetic energy conserved to 1e-12;
• GuidingCenter on IGAPolarTorus with AdhocTorus and polar splines, explicit and discrete-gradient pushers: energies per marker, compared by ID, including markers that cross the axis.

Also run locally:

• test_pushers.py, test_reflect.py, test_kernel_setup.py, test_kernel_backends.py, test_accum_gc_gradB.py, test_maxwellians.py and fields_background/tests (serial);
• the default-params runs of Vlasov and GuidingCenter.

A separate review agent checked the branch before the push. Its findings are fixed in the second commit:

• the Python axis map caught η1 > 1;
• the clamp covered only one loading branch;
• the GC gradients had dropped the boundary operators.

🤖 Generated with Claude Code

Follow-up (stacked): #762 adds the Vlasov/GuidingCenter profiling examples that use this.