PaperlayerView exact revision on arXiv

Source-linked AI summary

Machine-learned spectral PDE solvers for fluids and plasmas

Mellow Baixuan Chen, Alan A. Kaptanoglu

arXiv:2610.27317v1physics.comp-phcs.LGphysics.plasm-ph

TL;DR

PDE discovery from data is sensitive to observation noise, and derivative fitting does not ensure accurate forward integration. WARP trains explicit PDE coefficients through weak-form spectral Galerkin rollouts anchored at observed states. Across eight fluid and plasma systems, it improves coefficient recovery, rollout accuracy, and statistical measures over derivative fitting on the same weak formulation.

  • Problem

    PDE discovery from trajectories is confounded by observation noise, while derivative-fit methods do not ensure that fitted models reproduce trajectories under forward integration.

    Source: §1 · p. 1

  • Method

    WARP combines weak spatial projection in a spectral Galerkin solver with coefficient optimization through multistep rollouts anchored at observed states.

    Source: §1 · p. 1, §6 · p. 21

  • Results

    Across eight fluid and plasma systems, WARP improves coefficient recovery, rollout accuracy, and prediction of statistical measures over derivative fitting on the identical weak formulation.

    Source: §6 · p. 21

  • Takeaways & Limitations

    The comparison separates complementary robustness contributions: weak projection limits derivative amplification of observation noise, while multistep rollout improves integrated-dynamics accuracy.

    Source: §6 · p. 21

    Experiments assume periodic domains and candidate libraries containing the correct operators; experimental noise, missing fields, and shocks remain untested.

    Source: §5 · p. 21

Abstract

from arXiv · show

Sparse regression methods for discovering governing partial differential equations (PDEs) from data are sensitive to observation noise and are not constrained to reproduce the data under forward integration. We address both limitations by performing weak-form discovery in the spatial discretization in which the PDE is solved, learned by differentiating through multistep rollouts anchored at the observed states. The new method, Weak-form Anchored Rollout for PDEs (WARP), produces a machine-learned PDE solver. In this work, we focus on systems amenable to a spectral Galerkin solver. Across eight benchmark fluid and plasma systems, including high-dimensional turbulent simulations, WARP produces numerically stable rollouts from new random initial conditions and substantially improves coefficient recovery, rollout accuracy, and statistical measures compared with a single-step weak-form SINDy baseline.

1 Introduction

PDE discovery is confounded by observation noise and by methods that do not ensure accurate forward integration. WARP trains coefficients through anchored multistep rollouts of a weak-form spectral Galerkin model, improving results across eight fluid and plasma systems.

1 Introduction

  • Observation noise and inverse-problem ambiguities complicate PDE discovery, while small derivative-fit residuals do not ensure stable forward integration.Finite trajectories may visit only a small part of state space, where candidate operators can be correlated and distinct equations can fit local data.Source: §1 · p. 1
  • WARP combines weak spatial projection with solver-in-the-loop coefficient training, differentiating through multistep rollouts of a spectral Galerkin model initialized at observed states.It outputs explicit PDE coefficients and an autonomous solver in the same finite-dimensional weak formulation used during training.Source: §1 · p. 1
  • Across eight systems, WARP improves coefficient recovery and typically autonomous rollout accuracy over single-step weak-form fitting, isolating the benefit of multistep rollout.Both methods use the same weak spatial projection.Source: §1 · p. 2

2 Background and Motivation

The paper places WARP among methods that learn interpretable dynamics through integration, while distinguishing its specified physical-operator coefficients and weak spatial formulation from related approaches.

2 Background and Motivation

  • WARP estimates coefficients of specified differential operators, unlike neural representations that learn a vector field or spatial operators without naming a compact physical-operator set.Neural methods can be more expressive when relevant terms are unknown, while WARP’s target is an explicit PDE.Source: §2 · p. 2
  • Neural ODEs and related methods train by differentiating through numerical integration, evaluating the learned right-hand side on its evolving predicted state [19] [13] [24] [25] [26] [27].WARP uses the same self-consistency condition, rather than evaluating candidate operators only on measured states.Source: §2 · p. 2
  • Operator inference with rollouts already uses observation-anchored multistep integration; WARP instead ties shared PDE coefficients across retained modes through weak projections of named operators [16, 28].The distinction is parameterization and spatial weak formulation, not rollout training itself.Source: §2 · p. 2
  • Integration-aware sparse-identification methods include short-prediction SINDy, Runge–Kutta SINDy [8], and Unrolled-SINDy; WARP fixes analytic basis derivatives and weakly projects its PDE.Unrolled-SINDy evaluates candidate libraries on intermediate predicted states.Source: §2 · p. 2
  • Solver-coupled correction methods learn effective corrections or closures, whereas WARP targets a small vector of explicit PDE coefficients, which can absorb unresolved-scale effects in projected-DNS experiments.The distinction is less sharp for the projected-DNS experiment.Source: §2 · p. 2
  • Weak projection does not remove statistical bias: noisy responses and libraries remain coupled, so WARP’s noisy anchors and targets do not make it an unbiased estimator [49].WARP avoids a finite-difference time derivative, but noise-dependent biases from nonlinear library evaluation remain relevant.Source: §2 · pp. 2–3
  • WARP’s stated contribution is explicit PDE-coefficient identification using a shared Fourier representation for state and tests, not the novelty of combining weak forms with rollout training.On periodic domains, this choice removes boundary terms and avoids localized-test placement and shape choices.Source: §2 · p. 3
  • The motivating application scope includes plasma PDE discovery, MHD, kinetic-to-fluid modeling, boundary control, and reduced models of oscillations and turbulence [55] [57].Source: §2 · p. 3

3 Methodology

The methodology projects candidate PDEs into a finite spectral Galerkin system and trains their coefficients through anchored rollouts, then evaluates learned solvers on training and held-out trajectories or turbulent statistics.

3.1 Spectral Galerkin discretization

  • Periodic states are projected onto a finite orthonormal Fourier band used for both representation and testing; band-closed data evolve this Galerkin system without discarded-mode feedback.Nonlinear products can still generate modes outside the retained band, which projection and dealiasing remove.Source: §3.1 · p. 3

3.2 Weak-form modal ODE

  • Integration by parts moves derivatives onto known basis functions, yielding a weak-form modal ODE whose library defines the synthetic-data, training, and deployment model.Source: §3.2 · p. 4

3.3 Differential-operator representations

  • Candidate differential operators form an interpretable library whose coefficients are estimated; experiments use the correct operator support, leaving term selection and larger or misspecified libraries for future work.The fixed support isolates the effect of the training criterion.Source: §3.3 · p. 4

3.4 Observation model and noise

  • Each field’s grid observations receive independent additive Gaussian noise, and modal observations are obtained by projecting the noisy fields onto the basis.The noise-to-signal ratio is defined per field.Source: §3.4 · p. 4

3.5 Data generation

  • Reference trajectories and model evaluations use adaptive Dormand–Prince integration and 2/3 dealiasing; one experiment also uses externally generated DNS data.Source: §3.5 · p. 4

3.6 Single-step derivative-fit comparison

  • The FD baseline regresses central finite differences against the same library, basis, and modal observations as WARP, isolating multistep rollout from single-step derivative fitting.Source: §3.6 · p. 4

3.7 WARP anchored multistep rollout training

  • WARP initializes rollouts at observed anchor states, integrates to selected offsets, differentiates the objective through the solve, and outputs an autonomous spectral Galerkin solver.Offsets may be bidirectional or forward-only; the Hasegawa–Wakatani setup uses a shorter backward horizon.Source: §3.7 · pp. 4–5

3.8 Evaluation metrics

  • Deterministic held-out tests assess trajectory error, whereas turbulent benchmarks use invariant statistics after spin-up and compare against the true model’s finite-sampling floor.The protocol uses fresh initial conditions from the trained regime; turbulent evaluations include independent or perturbed draws as specified for each system.Source: §3.8 · pp. 5–6

4 Benchmarks and Results

The benchmarks evaluate noisy coefficient recovery, rollout accuracy, and turbulent-flow statistics across nine fluid and plasma configurations, with FD and WARP trained at the same sampling rate and adaptive Dormand–Prince integration throughout. Across the modified Hasegawa–Wakatani and Orszag–Tang benchmarks, WARP improves coefficient recovery and deployed behavior over finite-difference discovery, while results show that coefficient error alone does not determine invariant-statistics accuracy and that Orszag–Tang event accuracy degrades at high noise.

4.2 Burgers’ equation

  • On Burgers, WARP’s coefficient advantage grows with noise to about an order of magnitude at the highest level, though FD generalizes better at low noise (0.01 ≤ η ≤ 0.05).At the highest noise, WARP recovers convection within a few percent with a seed-to-seed spread more than an order of magnitude smaller, while FD’s estimate is strongly seed-dependent.Source: §4.2 · p. 7, Figure 1 · p. 8

4.3 Korteweg–de Vries equation

  • For KdV, WARP improves coefficient error by one to two orders of magnitude and yields one to three orders of magnitude lower testing rollout error at every nonzero noise level.At the highest noise level, FD loses the dispersive coefficient almost entirely, while WARP recovers both coefficients.Source: Figure 3 · p. 8, §4.3 · p. 8

4.4 Zakharov–Kuznetsov equation (2D)

  • For 2D ZK, WARP’s coefficient advantage grows monotonically under noise to roughly an order of magnitude, with improved integrated predictions from clean initial conditions.The testing curves compare both models from the same clean initial states, isolating the prediction improvement from initial observation noise.Source: Figure 5 · p. 9, §4.4 · p. 9

4.5 Decaying 2D incompressible Navier–Stokes vorticity

  • In decaying 2D Navier–Stokes, WARP’s coefficient advantage grows from one to nearly two orders of magnitude with noise; at the highest level, it estimates advection within a percent.FD underestimates the advection coefficient by about half; testing from identical clean initial conditions also shows improved predictions with WARP.Source: §4.5 · p. 9, Figure 7 · p. 10

4.6 Body-forced Kolmogorov turbulence (2D)

  • Under noise, FD coefficient error approaches unity as all three coefficients attenuate together, while WARP degrades gradually and preserves coefficient signs and magnitudes; through η = 0.1, WARP remains statistically indistinguishable from the true model, while FD drifts past laminar energy.Source: §4.6 · p. 10

4.7 Hasegawa–Wakatani coupled-field discovery in saturated turbulence

  • At c1 = 1, WARP improves saturated-HW coefficient error by one to two orders of magnitude under noise and matches the true-coefficient rollout through twenty training horizons.Its error curve coincides with the true-coefficient baseline deployed from the same noisy snapshot; FD’s error is several times larger during growth and saturates earlier.Source: §4.7 · pp. 12–13, Figure 16 · p. 13, Figure 17 · p. 13

4.8 Closure discovery from independent 5122 DNS data

  • On independent-DNS HW data, WARP has lower coefficient error than FD at c1 = 1 for every retained band, and its zero-noise offset decreases as M increases.Under observation noise, FD error rises while WARP remains near its zero-noise value; the offset reflects closure and discretization differences, not a proven lower bound on estimation error.Source: §4.8 · p. 15

4.9 Modified Hasegawa–Wakatani

  • The mHW solver scan locates an eddy-to-zonal transition and identifies four regimes for discovery: C ∈ {0.01, 0.1} (eddy), C = 0.3 (transition), and C = 1 (zonal).Doubling the band moves the transition close to the DNS transition and brings the zonal-regime flux exponent within a tenth of the DNS value.Source: §4.9 · p. 15
  • At 10% noise, WARP retains roughly one to two orders of magnitude lower mHW coefficient error than FD, and at 5% its deployments are statistically indistinguishable from the true model for C ≥ 0.1.At 10% noise, residual WARP coefficient error shifts the zonal fraction upward and can cross the threshold near the transition; invariant statistics are not a monotone function of coefficient error, and FD differs most decisively in the zonal regime.Source: §4.9 · pp. 15–16, Figure 23 · p. 16, Figure 24 · p. 17

4.10 Decaying MHD and the Orszag–Tang vortex (2D)

  • The Orszag–Tang generator’s M = 21 peak-current amplitude is 4% low relative to the 128² pseudo-spectral reference, with the first-peak time agreeing.The transient’s joint dissipation peaks at t∗ = 1.05.Source: §4.10 · p. 18, Figure 28 · p. 19
  • At η = 1%, WARP’s Orszag–Tang event errors are a few percent, while FD shifts the dissipation peak by an amount comparable to t∗ and overestimates peak current by 5–24× for η ≤ 0.1.At η = 20%, WARP’s event errors grow to order one and the event is no longer resolved; its held-out deployment reproduces the reference vorticity and current-density structure, including the current sheet, within the band truncation.Source: §4.10 · pp. 18, 20, Figure 30 · p. 19, Figure 31 · p. 20

4.11 Robustness checks

  • ZK coefficient recovery is statistically flat across tested basis sizes and flat in rollout horizon beyond H ≈ 25, placing the default H = 50 on the plateau.Because discovery and data generation share the discretization, no truncation floor arises; neither retained mode count nor horizon acts as a delicate tuning parameter in the tested cases.Source: §4.11 · p. 20

5 Limitations and Future Work

The comparison isolates multistep rollout training, not the separate effect of weak spatial projection. The experiments also assume periodic domains and correct candidate operators, while turbulent cases retain specific closure and resolution limits.

5 Limitations and Future Work

  • The shared Galerkin representation, preprocessing, and library isolate the effect of multistep rollout; a strong-form control is needed to quantify weak spatial projection separately.Comparisons with other weak-form implementations would also mix training criteria with different test functions and discretizations.Source: §5 · p. 21
  • The experiments assume periodic domains and libraries containing the correct operators; noise, missing fields, and shocks remain untested.Larger libraries require sparsity promotion or stability selection, and nonperiodic domains require an appropriate Galerkin basis.Source: §5 · p. 21
  • The turbulent studies lack closure terms for independent DNS projection, and the Orszag–Tang current-sheet amplitude is limited by the M = 21 band.The truncated, resolution-matched system captures the qualitative mHW transition with a quantitative offset.Source: §5 · p. 21

6 Conclusion

WARP learns explicit PDE coefficients by differentiating through anchored multistep trajectories of the weak Galerkin solver, yielding an autonomous spectral solver. Across eight fluid and plasma systems, it improves coefficient recovery, rollout accuracy, and statistical prediction over derivative fitting, while broader operator libraries and unknown dynamics remain to be validated.

6 Conclusion

  • WARP differentiates through anchored multistep trajectories of the weak Galerkin model to identify interpretable PDE coefficients.The resulting coefficients define an autonomous spectral solver for direct forward integration.Source: §6 · p. 21
  • Across eight fluid and plasma systems, WARP improves coefficient recovery, rollout accuracy, and statistical-measure prediction over derivative fitting on the identical weak formulation.The advantage persists under strong noise, in coupled turbulence, a current-sheet-forming MHD transient, and coefficient error for unclosed 5122 DNS projections.Source: §6 · p. 21
  • The approach combines spatial weak projection, which limits derivative-noise amplification, with multistep rollout, which constrains and improves integrated-dynamics accuracy.Next steps include automated term selection, large operator libraries, nonperiodic and local basis functions, and a profiled FFT implementation for larger modal systems.Source: §6 · p. 21

Declaration of generative AI and AI-assisted technologies in the writing process

The authors report using large language models as coding assistants and for manuscript proofreading and derivation checks, with all code and text reviewed and edited by the authors.

Declaration of generative AI and AI-assisted technologies in the writing process

Appendix A Additional Hasegawa–Wakatani Diagnostics

The appendix supplements the main-text Hasegawa–Wakatani results with diagnostics and deployment visuals for additional settings. It also explains how multistep rollout training couples offsets through one trajectory and penalizes models that drift during integration.

Appendix A Additional Hasegawa–Wakatani Diagnostics

  • The appendix retains Hasegawa–Wakatani diagnostics for c1 = 0.1 and c1 = 5, alongside spatial snapshots of the c1 = 1 deployment.The main text shows identification and validation at c1 = 1.Source: Appendix A · p. 21, Figure 32 · p. 22, Figure 34 · p. 23
  • Figure 32 plots particle flux, resistive energy dissipation, shell energy, density fluctuation, and vorticity for band-closed HW data generation.Its top and bottom correspond to c1 = 0.1 and c1 = 5, respectively, using resolution-matched hyperdiffusion.Source: Figure 32 · p. 22
  • Figure 33 compares clean reference and WARP and FD deployments from a noisy snapshot at t = 250, at 10 and 20 time units after launch.The columns represent 20 and 40 training horizons; color scales are set by the reference at each time and shared down each column.Source: Figure 33 · p. 22

B.1 Why multistep rollout training helps

  • Finite temporal offsets increase signal relative to snapshot noise, while the weak library is evaluated at predicted states during every right-hand-side call.Coupling offsets through one model trajectory penalizes coefficients that match local increments but drift under integration, rejecting many unstable or rapidly diverging forward solutions.Source: §B.1 · p. 23

B.2 Coefficient recovery and noise sensitivity

WARP's coefficient-recovery advantage varies across benchmarks and noise regimes, with particularly large gains reported for KdV. Sampling rate can bias finite-difference fits when important dynamical rates are unresolved, whereas WARP integrates candidate models adaptively between observations.

B.2 Coefficient recovery and noise sensitivity

  • 1–2 orders of magnitude on KdV, versus a modest Burgers advantage, is an empirical trend rather than a scaling law.The problems also differ in resolution and time scale.Source: §B.2 · p. 23
  • FD degrades gradually on Burgers, ZK, Navier–Stokes, and HW, but becomes noise-dominated on KdV, Kolmogorov flow, and the Orszag–Tang vortex.After FD saturates, more noise mainly increases WARP's error and narrows the methods' ratio, although WARP remains substantially more accurate; central differences also carry sampling-rate bias in Kolmogorov flow and mHW.Source: §B.2 · pp. 23–24

Appendix C Sampling-Rate and Resolution Checks

  • A 128^2, M = 42 Kolmogorov-flow companion changes mean energy by 0.34% and leaves mean dissipation within block-sampling error.The 64^2, M = 21 case reaches energy near 0.44Elam and dissipation and injection near 0.09Dlam, without band-edge pile-up.Source: Appendix C · p. 24
  • Within 2%, a small-amplitude unstable eigenmode reproduces the mHW dispersion-relation growth rate.Validation uses boxes containing ten fastest-growing wavelengths and resolution-matched dissipation in both fields.Source: Appendix C · p. 24
  • Unresolved dynamically important rates bias central-difference fits even without observation noise, while WARP's adaptive integration decouples integration accuracy from sampling rate.In mHW, the resistive relaxation rate is largest at the smallest nonzero wavenumbers; in Kolmogorov flow, the fastest viscous decay is at high wavenumbers where the measured spectrum is weak.Source: Appendix C · p. 24

Appendix D Hyperparameter Sensitivity

The retained sensitivity runs show little change in ZK coefficient error across the band sweep and in Burgers error once H ≥25.

Appendix D Hyperparameter Sensitivity

  • Burgers coefficient error changes little once H ≥25.The sensitivity comparison varies rollout horizon H.Source: Appendix D · p. 24, Figure 36 · p. 24
  • ZK coefficient error is statistically flat across the retained-band sweep.Figure 36 shows error bars for one standard deviation over three noise seeds.Source: Figure 36 · p. 24, Appendix D · p. 24
  • Figure 36 compares WARP coefficient error against retained band M for ZK and rollout horizon H for Burgers.Error bars represent one standard deviation over three noise seeds.Source: Figure 36 · p. 24

Appendix E Long-Horizon Deployment

In long-horizon deployment from unseen initial conditions, WARP remains stable and follows reference behavior across three systems, while FD’s bounded trajectories can still show incorrect dynamics.

Appendix E Long-Horizon Deployment

  • The comparison uses η = 10% models launched from a new initial condition for twice the original testing interval, against a reference integrated from that same condition.Figure 37 shows the deployment; HW uses c1 = 1, and |q|max denotes the largest absolute modal coefficient.Source: Appendix E · p. 24, Figure 37 · p. 25
  • WARP stays stable and follows the reference in all three systems during deployment over twice the original held-out horizon.The models were trained at η = 10% and deployed from an initial condition not used during training.Source: Appendix E · pp. 24–25, Figure 37 · p. 25
  • WARP reproduces the Orszag–Tang energy decay too quickly because viscosity and resistivity are overestimated, despite accurate advection and Lorentz coefficients.In HW, saturated fluxes match the reference; in Kolmogorov flow, the energy budget and intermittent dissipation bursts are reproduced.Source: Appendix E · pp. 24–25
  • FD remains bounded in every case but either develops spurious energy growth or settles on the wrong attractor over long horizons.In Orszag–Tang, total energy grows by an order of magnitude; in Kolmogorov flow, the dynamics run about an order of magnitude too slowly and miss the transition.Source: Appendix E · p. 25

Appendix F Identification with an Overcomplete Library

With an overcomplete 33-column library for Kolmogorov flow, WARP retains accurate coefficient recovery under noise and reproduces the energy budget at deployment, unlike FD’s attenuated and spurious coefficients.

Appendix F Identification with an Overcomplete Library

  • The benchmark uses 33 zero-mean-projected columns up to quadratic field products and second-order spatial derivatives, with five columns belonging to the true model.The true advection, Laplacian, and forcing terms are represented by the antisymmetric derivative pair, two second-derivative columns, and forcing, respectively.Source: Appendix F · p. 25
  • With noise, FD attenuates true coefficients and assigns weight to spurious columns at every noise level; above one percent noise, its spurious weight exceeds its viscosity coefficient.The spurious terms are mostly an opposite-sign pair of linear terms in ω and ψ, amounting to high-wavenumber linear damping in place of attenuated viscosity.Source: Appendix F · pp. 26–27, Figure 38 · p. 26
  • WARP recovers true coefficients to within a few percent at all noise levels and keeps spurious coefficients an order of magnitude smaller than FD’s.At ten percent noise, WARP’s largest spurious coefficients are as large as its recovered viscosity.Source: Figure 38 · p. 26, Appendix F · p. 27
  • At η = 10%, WARP reproduces the true energy budget, including intermittent dissipation bursts, while FD evolves several times too slowly and settles at the wrong energy.Both models remain bounded from a new initial condition; boundedness is observed for these fitted models, not guaranteed for arbitrary coefficients in the overcomplete library.Source: Appendix F · p. 27, Figure 39 · p. 27

Related papers

Loading 2610.27317v1…