Hybrid-PIC: operator-split implicit magnetic diffusion with exact curl-curl preconditioning

[tomzhu0225]

Summary

This PR adds an optional, default-off operator-split implicit magnetic diffusion step for hybrid-PIC.

In boundary-driven hybrid problems with a large vacuum (or vacuum-like) resistivity, leaving resistive diffusion of B entirely in the explicit Faraday / Ohm path makes the vacuum region stiff: substeps collapse and the driven B signature tends to lock to the wall instead of the plasma edge. The new split advances

dB/dt = -curl((eta / mu0) curl B)

after each hybrid B half-step with a theta-method linear solve. Ohm's-law resistivity used by the explicit Faraday update is capped by mag_diff_eta_explicit_max (default 0) and the implicit solve advances the residual max(eta - mag_diff_eta_explicit_max, 0), so resistivity is not applied twice.

Zero places all resistive diffusion in the implicit solve. In coupled production problems, a small positive, explicitly CFL-safe cap can improve robustness by retaining controlled resistive damping in the explicit substeps. The appropriate value is problem-, mesh-, and timestep-dependent; the Z-pinch validation below used 2.0e-2 Ohm m.

The split is inspired by the FLASH MagDiff approach; the implementation is native WarpX/AMReX. It does not make the Hall or electron-pressure terms implicit.

Motivation and application validation: RCYL dynamic Z-pinch

The motivating application is a 164 ns RCYL dynamic Z-pinch with an outer-radius B_theta drive (pec_insulator), a 15 mm deuterium column, and a high-resistivity vacuum gap. Without implicit magnetic diffusion, the hybrid field update is either unstable at the required vacuum resistivity or does not transmit the applied field from the wall to the moving plasma boundary correctly.

The current comparison uses the calibrated hybrid run rather than the earlier RZ demonstration:

	Hybrid with this PR	Fully kinetic reference
Geometry and duration	RCYL, 0-164 ns	RCYL, 0-164 ns
Model	Hybrid ions + QDSMC electron energy	Theta-implicit kinetic electrons and ions
Mesh and particles	128 radial cells; 400 ion particles/cell	224 radial cells; 400 particles/cell/species
Time step	7.826 ps	1 ps
Resistive field treatment	Variable-eta operator split with PETSc exact curl-curl Pmat; explicit eta cap 2.0e-2 Ohm m	Fully implicit EM response
Observed step cost	9-11 ms during run-in; 16.1 ms full-run average with D-D collisions	About 0.1 s/step

The timings above are application measurements from different hosts, not a controlled same-node microbenchmark. They nevertheless describe the practical runs: the hybrid also takes 7.8x fewer steps. The corrected full hybrid run completed in 338 s locally; extrapolating the measured implicit step cost gives roughly 1.64e4 s for the reference trajectory, about a 50x observed time-to-solution difference for this setup.

Magnetic diffusion is responsible for the close B_theta front and vacuum-field profile. Matching the density sheath additionally requires deuterium-deuterium Coulomb collisions; resistive electron-ion drag alone does not supply ion-ion shock broadening. With both pieces present:

• Over 8-80 ns, the normalized RMS errors are 0.114 for B_theta and 0.113 for density.
• Over 8-120 ns, including first stagnation, they are 0.174 and 0.192.
• Both models reach the axis at the 112 ns diagnostic. At 120 ns the peak density is 39.5 n0 in hybrid and 40.1 n0 in the kinetic reference.
• Over the complete 8-164 ns record, the normalized RMS errors are 0.221 and 0.205. The remaining error is concentrated in the post-stagnation rebound: the trajectories begin to drift out of phase after about 144 ns.

The figure shows the full trajectory rather than only the favorable run-in interval. The upper row is deuterium density and the lower row is B_theta.

Scope

Capability	Status
Default off; existing hybrid decks unchanged	Yes
Geometries	Cartesian 1D_Z / XZ / 3D, RZ, and RCYLINDER; RSPHERE is rejected
Spatially varying parser eta	Yes, through the matrix-free operator in all supported geometries
Constant eta	Flat MultiFab through the same matrix-free path
Time integration	Theta method for 0 < theta <= 1; backward Euler is the default and Crank-Nicolson is tested
Physical field boundaries	Periodic, PEC, and pec_insulator in Cartesian 1D/XZ/3D; axis, PEC, and pec_insulator handling in RZ/RCYLINDER
Embedded boundaries	Stair-step EB in XZ, 3D, and RZ. AMReX does not provide EB in 1D builds, so Cartesian 1D and RCYLINDER cannot enable EB
Linear solvers	AMReX GMRES with scaled Jacobi; optional PETSc KSP
PETSc preconditioner matrix	Assembled frozen-eta exact discrete curl-curl operator in Cartesian 1D/XZ/3D, RZ, and RCYLINDER; covered EB degrees of freedom are omitted
Linear tolerances used in regression tests	rtol=1e-8, atol=0
PICMI	Yes - all magnetic-diffusion controls, including PETSc algebraic-PC controls, are exposed on HybridPICSolver

Main inputs

hybrid_pic_model.implicit_mag_diffusion = 0|1   # default 0
hybrid_pic_model.mag_diff_theta = 1.0            # (0, 1]; 1 = backward Euler
hybrid_pic_model.mag_diff_eta_explicit_max = 0.0 # explicit Ohm cap; 0 = fully implicit
hybrid_pic_model.mag_diff_use_variable_eta = 0|1 # parser / MultiFab eta
hybrid_pic_model.mag_diff_constant_eta = ...     # optional constant eta
hybrid_pic_model.mag_diff_linear_solver = amrex_gmres|petsc
hybrid_pic_model.mag_diff_rtol / atol / max_iter / verbose

# Optional PETSc algebraic preconditioner controls
hybrid_pic_model.mag_diff_petsc_pc_type = asm
hybrid_pic_model.mag_diff_petsc_asm_overlap = 0
hybrid_pic_model.mag_diff_petsc_sub_ksp_type = preonly
hybrid_pic_model.mag_diff_petsc_sub_pc_type = ilu
hybrid_pic_model.mag_diff_petsc_ilu_factor_levels = 2

The same controls are available through PICMI:

solver = picmi.HybridPICSolver(
    grid=grid,
    Te=Te,
    n0=n0,
    plasma_resistivity=eta,
    implicit_mag_diffusion=True,
    mag_diff_theta=1.0,
    mag_diff_use_variable_eta=True,
    mag_diff_eta_explicit_max=0.0,
    mag_diff_linear_solver="petsc",
    mag_diff_rtol=1.0e-8,
    mag_diff_petsc_pc_type="asm",
    mag_diff_petsc_asm_overlap=0,
    mag_diff_petsc_sub_ksp_type="preonly",
    mag_diff_petsc_sub_pc_type="ilu",
    mag_diff_petsc_ilu_factor_levels=2,
)

When variable eta is enabled, resistivity comes from the existing hybrid resistivity parsers and is frozen on Yee E/J faces for the diffusion solve. PETSc always receives the exact assembled curl-curl Pmat matching the matrix-free operator. The algebraic preconditioner applied to that matrix is configurable through the mag_diff_petsc_* inputs or PETSc runtime options.

Implementation

• HybridMagDiffusion applies the theta method to the matrix-free FDTD curl(eta curl B) / mu0 operator.
• Constant and parser-defined eta share one operator path; the earlier MLCurlCurl path has been removed.
• The solve runs after each hybrid BfieldEvolve half-step, with the explicit Ohm eta cap applied in HybridPICSolveE.
• AMReX GMRES uses a geometry-aware scaled Jacobi diagonal.
• PETSc KSP uses an assembled frozen-eta exact curl-curl Pmat, including Cartesian component coupling and the cylindrical metric/axis terms.
• Boundary values are eliminated consistently in the matrix-free operator and PETSc matrix. Stair-step EB masks keep covered fields fixed and remove covered PETSc unknowns.
• Under CUDA, EB masks and PETSc vectors use pinned host staging because the current PETSc path is host-side.

Validation

Examples/Tests/hybrid_mag_diffusion/ covers:

• Cartesian constant-eta Fourier decay and variable-eta Crank-Nicolson smoke tests. Python-enabled builds run the constant-eta case through HybridPICSolver; non-Python builds use the equivalent native deck.
• AMReX/PETSc parity and exact-Pmat checks in Cartesian 1D, XZ, and 3D.
• Cartesian PEC / pec_insulator boundaries in 1D, XZ, and 3D, with XZ and 3D EB coverage.
• RZ variable eta, Crank-Nicolson, PETSc, EB, split parity, radial feed, axial PEC, and axial feed cases.
• RCYLINDER mild, harsh, PETSc parity, harsh-PETSc, and variable-eta pec_insulator feed cases.

After synchronizing with current development, all 29 local non-PETSc magnetic-diffusion run/analysis tests and all 10 PETSc-enabled RCYLINDER tests passed. The full GitHub Actions matrix remains the merge gate.

Follow-up work

1. Benchmark and document scalable PETSc algebraic preconditioners, including ASM/BJacobi/Hypre/GAMG, and investigate reuse of the assembled preconditioning matrix when eta changes slowly. User-selectable PETSc preconditioner configuration is already implemented.
2. Quantify the accuracy, stiffness-damping, and runtime trade-off between backward Euler (theta=1) and Crank-Nicolson (theta=0.5). Add old/new boundary-value blending for time-dependent pec_insulator feeds; driven CN boundaries are currently first-order in time.
3. Reduce host/device staging with device-aware PETSc assembly and solve paths.

[tomzhu0225]

The application validation in the PR description has been refreshed from the completed RCYL study. The obsolete 0-50 ns RZ figure and timing estimate are replaced by the 0-164 ns RCYL density/Btheta comparison, current 9-16 ms hybrid timing measurements, and the full-trajectory error metrics. The text now also separates the roles of implicit magnetic diffusion (vacuum-field evolution) and D-D collisions (ion-shock broadening), and explicitly documents the remaining post-stagnation phase drift after about 144 ns. The timing comparison is labeled as cross-host application performance rather than a same-node microbenchmark.

[tomzhu0225]

PICMI support is now included. HybridPICSolver exposes implicit_mag_diffusion and all current mag_diff_* controls while leaving unspecified values to the C++ defaults. The existing 2D constant-eta analytic regression selects a PICMI deck in Python-enabled builds and the equivalent native deck otherwise, so API coverage was added without adding another CI test. Local checks: all ten PICMI values emit correctly; the generated deck passes the analytic decay oracle (discrete relative error 4.1e-16); native fallback CTest passes 2/2; Ruff and check_inputs.py pass.

[ax3l]

ping @prkkumar for coordination with existing implicit implementation

[tomzhu0225]

Thanks @ax3l. Tagging @prkkumar-he and @RemiLehe for coordination. We are aware of the fully implicit work in #7208.

This PR implements the standard operator-split/IMEX middle ground:

• A fully explicit advance is inexpensive per step, but resistive magnetic diffusion imposes a parabolic timestep limit, $\Delta t \propto \Delta x^2/\eta$, which becomes prohibitive for large resistivity or vacuum-like regions.
• A fully implicit advance can remove broader stability restrictions, but requires solving a larger coupled particle-field nonlinear system.
• Here, only the stiff magnetic-diffusion operator is implicit; the established explicit hybrid particle-field advance remains unchanged. This removes the diffusion timestep bottleneck through a smaller linear solve without requiring a new coupled evolution scheme.

Treating the stiff diffusive term implicitly while retaining non-stiff terms explicitly is a standard IMEX/operator-splitting strategy. The exact curl-curl preconditioner makes that targeted solve practical for spatially varying resistivity. Since the feature is opt-in, preserves existing defaults, and has analytic, regression, PICMI, and application validation, it can be merged as an independently useful capability.