High-order polar splines via homotopy operators (PolarFEEC.jl port)

[max-models]

High-order polar splines of arbitrary regularity 𝔯 at the pole, built with Poincaré homotopy operators: a port of R. Remmerswaal's PolarFEEC.jl and of the methods in the NMPP seminar talk (2026). They are available as

DerhamOptions(polar_splines=True, polar_method="homotopy", polar_regularity=r)

C¹ (polar_method="c1") remains the default and is unchanged.

Solves the following issue(s):

None

Serial core: struphy.polar.homotopy

A pure numpy/scipy port of the algorithms of PolarFEEC.jl, with its conventions.

• spaces1d: clamped and periodic B-splines and trigonometric spaces, each with a derivative-scaled companion and exact derivative matrices.
• forms2d: tensor-product and physical 2D form spaces, with commuting and L² projections.
• mappings: disk-like mappings.
• polynomials: the Zernike basis and the centre-shifted polynomial basis.
• complex2d: PolarComplex2D, which builds:
  • the extraction matrices E and the duals F (FᵀE = I);
  • the enrichments cochain_chain, cochain and none;
  • the homotopies, the Dirichlet boundary homotopies and the mean projection 𝔅;
  • the cohomology projections and the polar masses and derivatives.
• complex3d: PolarComplex3D, the total complex with the 1D complex in η₃. It has the cohomology of the solid torus: [1,1,0,0] free and [0,0,1,1] Dirichlet.
• solvers: Poisson and curl–curl in the homotopy gauge. The 2D curl–curl "free lunch" is one M̄₂ solve followed by h₂.

Core changes

• DerhamOptions: new fields polar_method, polar_regularity, polar_enrichment, polar_polynomial_projection and polar_cross_section.
• PolarExtractionBlocksHomotopy (polar/homotopy/derham_blocks.py) provides the blocks of struphy.polar.linear_operators for the 3D total complex. In PolarDerhamSpace:
  • the ring layout comes from these blocks;
  • the C¹ checks (Nel2 % 3 == 0, p ≥ 2) apply only to C¹;
  • the domain-decomposition check is collective and suggests mpi_dims_mask=(False, True, True).
• CommutingProjector: the homotopy path is the tensor-product projection followed by F̃ᵀ. The polar inter-/histopolation matrix is the identity, and no polar_shift is needed. The C¹ path is unchanged.
• Any domain with a pole works: HollowCylinder(a1=0), HollowTorus(a1=0), IGAPolarTorus, and the new domains:
  • PerturbedDiskCylinder (kind_map 23) and RoundedSquareCylinder (24), with pyccel and CUDA kernels;
  • HelicalTube, a spline mapping with the presets "stellarator" and "islands".
• New equilibrium HelicalVacuumField: a single helical harmonic, B = B₀∇(z + δ I_m(kr) sin(mθ − kz)/k). It evaluates the Bessel functions with SciPy, so on GPU it is evaluated on the host, like GVEC and DESC.
• Fixes:
  • The Cartesian equilibria (CartesianFluidEquilibriumWithB, CartesianMHDequilibrium) returned NaN on a polar axis. They now evaluate b, j and ∇|B| directly in Cartesian coordinates; away from the axis the values are unchanged (≤ 5e-15).
  • grad_1d_matrix now works for n = 1.
  • The polar source of ImplicitDiffusion is fixed.

Model-specific changes

• New toy model ScalarWave: the 0-form wave equation, with the propagator ScalarWaveLeapfrog.
• New MaxwellWeakAmpere(algo="leapfrog"): velocity Verlet, the symplectic scheme of the talk.

Examples

Each example has a physics description, --test-run and a pproc script with the standard snapshots.

• examples/ScalarWave/gaussian_perturbed_disk: the 0-form wave across the axis. It matches the core leapfrog to 1e-11.
• examples/Maxwell/gaussian_perturbed_disk: the 1-form wave. It matches the core to about 1e-5.
• examples/Maxwell/gaussian_helical_tube: the 3D 1-form wave in HelicalTube("stellarator").
• examples/Poisson/polar_rounded_square and examples/CurlCurl/polar_rounded_square: PolarFEEC.jl example 05. They match the core.
• examples/GuidingCenter/vacuum_near_axis: orbits near the axis of HelicalVacuumField and the energy error max_t |E(t) − E(0)| / (μ|B|). convergence_vacuum_near_axis.py runs the scan over 𝔯 = 0…4 and Δt = 0.04…0.0025.

Demos, tutorials and figures

struphy.polar.homotopy.demos reproduces the figures of the talk and of PolarFEEC.jl:

Demo	Figures
spaces	jl 01–03
basis, regularity	the bases and regularity at the centre (talk; jl 04)
projections	convergence and commutativity of the projections (talk)
pde	Poisson and curl–curl, jl 05
maxwell	the spectrum on the disk, jl 06
fieldlines	Poincaré sections and ι(r) in the twisted tube, jl 07
performance	jl 08
waves	the 0- and 1-form wave studies of the talk

The tutorials tutorial_polar_splines.ipynb and dev_tutorial_polar_homotopy.ipynb use these demos. By default they run at a coarse resolution for CI; STRUPHY_TUTORIAL_FULL=1 switches to the resolution of the talk.

Tests

In src/struphy/polar/tests/:

• test_homotopy_core: a port of the PolarFEEC.jl tests (paper theorems, polar complex, boundary homotopies, the 2D curl–curl case).
• test_homotopy_convergence (marker convergence): rates 4, 3, 3 on the perturbed disk.
• test_homotopy_3d: the total complex, the torus volume and the 3D Poisson problem.
• test_homotopy_derham: the Derham blocks equal the core matrices; commuting projectors, cohomology, masses, the preconditioner and Dirichlet, for 𝔯 = 1…4. Runs serially and on 4 MPI ranks.
• test_homotopy_demos.
• test_homotopy_applications:
  • energy conservation of the leapfrog;
  • ScalarWave equals the core leapfrog;
  • HelicalTube;
  • HelicalVacuumField is curl- and divergence-free and finite on the axis;
  • the Cartesian equilibria are unchanged off the axis;
  • field lines: div b ≈ 1e-15 and ι matches.

Every new example passes test_examples. The existing C¹ polar tests, fields_background tests (including a new HelicalVacuumField CuPy case), Maxwell and CurlCurl verification tests and toy-model tests all pass. Locally, the 2 DESCunit domain tests fail only because desc-opt is not installed.

The core was cross-checked against PolarFEEC.jl (Julia) on small cases: E, F, the homotopies, the dimensions and the Maxwell eigenvalues.

Documentation changes

• Polar splines (numerics-geomFE.rst): a new section on the homotopy construction, a paper ↔ code table and a study ↔ code table.
• Model and propagator lists: ScalarWave and ScalarWaveLeapfrog are added.

Reconstructed parts

The seminar's wave and guiding-centre scripts were not available. Those studies follow the methods stated in the talk and in the wave script from the PolarFEEC.jl history:

• the staggered leapfrog at 0.9 · 2/√λ_max;
• the Gaussian pulse next to the axis;
• the maximum relative L∞ error over the final 10%;
• RK4 orbits.

The following are choices made here and are listed in the docs:

• the reference solution (the finest mesh of the highest regularity);
• the L∞ sampling;
• T = 1;
• the helical-tube geometries;
• the vacuum field.

Known gaps (follow-ups)

• Not in this PR:
  • a solver="homotopy" option in the CurlCurl propagator (the free lunch in 2D; PCG on D̄ᵀM̄₂D̄ followed by h₂ in 3D), plus a homotopy case in test_verif_CurlCurl.py;
  • a homotopy variant of examples/CurlCurl/gaussian_on_axis_torus;
  • entries for the new domains in tutorial_domains.ipynb;
  • fake-CuPy tests for the homotopy blocks;
  • the trigonometric panel of the guiding-centre study, which needs a small core guiding-centre integrator.
• Not run: the long runs (the guiding-centre convergence scan, the wave studies and tutorials at full resolution, the animations) and a full Sphinx build.
• For plasma-plots: the curvilinear-grid plotting helper and the log-log convergence helper in demos/plotting.py would fit there better.

🤖 Generated with Claude Code

[max-models]

[ronaldremmerswaal]

It seems the code was ported entirely, without making use of struphy's existing B-spline functionality.

Similarly, it seems the existing polar construction still exists, which should (this should be verified) coincide exactly with setting regularity == 1 in the general framework.

It seems only the mappings that come with the PolarFEEC package are suppoted; instead the mappings from struphy's existing geometry should be supported.