[nicolominichelli]

Follows #900 (merged).

Summary

Second-order and stable time stepping of the perturbative (delta-f) variational MHD model VariationalPerturbativeMHD. The changes came out of a comparison with PerturbativeMHDWithQ (exact nonlinear slab solutions, screw pinch). Each change keeps the discrete equations of the variational principle and the exact energy conservation; only the grouping of terms into substeps and the nonlinear solvers change.

1. Sound step with the total pressure as transport weight (VariationalLinearSoundEvolve(total_weight=True), set by the perturbative model). The pressure perturbation is transported with p0 + p~^{n+1/2} in the implicit sound step, instead of p0 there plus a separate transport with weight p~ in the nonlinear substep. The pressure enters the energy only linearly (int p_hat / (gamma-1)), so the separate substep had no force: the nonlinear part of the sound coupling was split across two substeps, i.e. explicit. It was unstable for the grid-scale sound waves (omega dt >~ 2) as soon as p~ is a finite fraction of p0 (growth x1.1 per step at 5 %, x1.25 at 10 %, and the blow-up with a constant offset p0 -> p0 + pc, p~ -> p~ - pc). Now the step depends only on the total pressure. Picard iteration on the weight, one Schur solve per iteration (4-6 iterations).
2. VariationalNonlinearPBEvolve keeps only the magnetic terms and is solved by a Picard iteration on the weight b^{n+1/2} with one SPD solve of (M_rho + dt^2/4 K) du = ... per iteration. The previous Picard iteration on all of (u, b) was an explicit fixed point on the wave operator and diverged for dt k_max |b~| >~ 2 (measured contraction (dt |b~| k_max / 2)^2); it also stopped after one iteration when b~ = 0 (explicit Euler for p~).
3. VariationalDensityEvolve: midpoint density weight (Options(weights=...), default "midpoint" for model="perturbative", "frozen" otherwise): the frozen weight made the step non-symmetric and the model first order in time. New option jacobian_solver ("pbicgstab" default, "rpbicgstab", "gmres"); a breakdown of the Newton linear solve (BiCGStab at flow CFL numbers dt |u| k_max ~ 10) now raises an error instead of continuing with NaN.

Evidence

Scripts and data: delta_f_docs/var_vs_q/data (not in the repo).

• Slab, exact nonlinear solutions (circularly polarized Alfven wave 10 %/20 %, simple sound wave 10 %/30 %): second order on all dt ladders (dt = 0.2 ... 0.006), energy exact to <= 3e-10. Before: sound wave first order; Alfven wave failed for dt >= 0.1 (10 %) and dt >= 0.025 (20 %).
• Seeded stability (grid noise 1e-8): uniform p~ of +-10 %, +20 %, and offsets pc = 0.5, 1.0 up to dt = 0.2: stable (before: blow-up after 30-60 steps).
• Screw pinch (R0 = 20, 12x12x4): stage pc = 0 / 0.1 / 1.0 gives the same velocity to 3e-11 (before: 0.9 % at pc = 0.1, blow-up at pc = 1); agreement with the previous scheme where it was stable (5e-4); finite b~ (10 %, 20 %): same evolution as PerturbativeMHDWithQ; long run at b~ 20 %, dt = 0.2 (4x the q time step): bounded |u| and b~ as with q and energy to 3e-12 for 28 steps, then the density Newton solve broke down at |u| ~ 1.7 (flow CFL ~10, see Open; q at dt = 0.05 itself lost energy at the 1e-4 level and failed at t = 39).
• Vanishing edge pressure (p~ proportional to p0, p0 = 1e-8 at the edge, 10 %/20 %, positive total pressure): same evolution as PerturbativeMHDWithQ (max |u| agrees to 3 digits), energy exact to 1e-12.
• Cost (pinch): +1 to +5 % per time step (the sound step is 2-4x more expensive, the nonlinear substep cheaper; both are 11-15 % of the step).

Tests

• New src/struphy/propagators/tests/test_variational_sound_total_weight.py: the sound step with total weight is symmetric in time and conserves the energy; the Strang group at a large time step gives the same velocity for p0 + pc, p~ - pc as for pc = 0 (< 1e-7; 1.39 without the total weight).
• test_variational_nonlinear_pb_evolve.py updated (no pressure in the nonlinear substep).
• All variational propagator tests pass on 1 and 2 ranks, as do the variational model regression tests.
• Baselines VariationalPerturbativeMHD-u_space=Hdiv and VariationalPerturbativeMHD (default, u in H1vec) updated (1 and 2 ranks); in both the merged sound step changes the second-order diagnostic bt2 by 7e-3 and div_u by 7e-4 (relative), the splitting difference between two second-order schemes for the same equations (verified: switching the total weight off removes the nonlinear pressure transport and changes it by 10 %).

Open

• VariationalDensityEvolve: at flow CFL numbers of order ten the BiCGStab solve of the Newton system breaks down (now an error); unpreconditioned GMRES works but needs many iterations in 3D.
• Inner-wall mode of the linear model (Pi^1, coarse radial grid), unchanged.

🤖 Generated with Claude Code

[spossann]

Looks good - please update with current devel.