Neural operators applied to transient-dynamics PDEs with strong discontinuities exhibit autoregressive instability: in concrete-penetration stress-field prediction, the wavelet neural operator (WNO) diverges in autoregressive rollout, while MeshGraphNets collapse to zero predictions. WNO's instability stems from the lack of a structural constraint on the spectral radius of its propagation operator; the Fourier neural operator (FNO) is stable in these measurements but only emergently, not by construction. We propose a constitutive Markov physics-informed neural operator (MPNO) modeling one-step evolution as a Markov (row-stochastic) propagation operator. Physics-coupled edge weights (acoustic-impedance harmonic mean, contact area, and traction amplitude) encode material-interface constitutive information into a nonnegative symmetric adjacency matrix W; after normalizing the graph Laplacian L = D - W by lambda_max, the propagator P = I - alpha*L~ is constructively constrained to spectral radius rho(P) <= 1, suppressing exponential amplification of autoregressive errors. Stability is thus a designable architectural property, not an optimized loss objective. On three PDEs (Burgers and two-dimensional transverse-section concrete penetration), MPNO rolls out stably with bounded error on all test seeds at 100/135/165 m/s; the single-step relative L2 error is 0.7304 +/- 0.0008, better than WNO and comparable to FNO at about one quarter of FNO's parameters. The edge-weight formula transfers across scenarios by replacing material-property variables. With about 20K parameters, MPNO delivers roughly 10^5x inference speedup over LS-DYNA.
Figures & tables
Figure 1: Autoregressive rollout of neural operators on concrete-penetration stress-field prediction (29 steps, frames 11–39 of the 40-frame window; schematic, anchored to measured single-step errors and end-of-rollout error bounds). MPNO (solid blue) stays bounded throughout; WNO (red dashed) diverges to infinity on all seeds; MeshGraphNets (gray dotted) fails to learn (error stays near 1); FNO (orange band) remains stable on the clean test seeds but, unlike MPNO, carries no spectral guarantee.
Parameter
Value
Target
ϕ500×200 mm (cylinder)
Aggregate fraction / size
42% / 5–25 mm
Aggregate / mortar constitutive
JH-2 (granite) / CSCM (C30)
Projectile / velocity
ϕ45 mm ogive steel projectile / 100, 135, 165, 200 m/s (100 each)
Grid resolution / time frames
64×64 / 40 frames
ITZ
Homogenized dispersion ( ∼ 50 μ m)
Table 1: Parameters of the concrete-penetration dataset
Figure 2: Problem geometry and data. (a) LS-DYNA setup: a ϕ 500 × 200 mm cylindrical concrete target (Fuller-graded aggregates 5–25 mm, 42% volume fraction, JH-2/CSCM) impacted by a ϕ 45 mm ogive-nose steel projectile (0.8 kg, 100–200 m/s). (b) An LS-DYNA σzz stress-field snapshot of the impact zone (seed 5, 100 m/s, frame 20; normalized by 30 MPa, range ±8.5 ; peak σzz 5.039). (c) An example aggregate volume-fraction field Vf .
Figure 4: Spectral-radius taxonomy for autoregressive stability. A one-step propagator with spectral radius ρ<1 contracts ( ∥et∥2→0 ), ρ>1 diverges ( ∥et∥2→∞ ), and ρ=1 is the critical boundary (error bounded, ∥et∥2≤∥e0∥2 ). MPNO’s constructed propagator P=I−αL~ places ρ(P) exactly at the critical boundary (Property 1), independent of training. Measured full-model Jacobian spectral radii (gray markers, representative operating point seed 55 frame 10, § 4.4 Fig. 8 ): ρ(J)=0.83 (MPNO) and ρ(J)=0.90 (no_spec). Re-measured over six clean rollout seeds, ρ(J) is operating-point-dependent (MPNO mean ≈1.05 vs. no_spec ≈1.13 ); the ordering MPNO < no_spec persists, but a single operating point does not certify ρ(J)≤1 . POD-DMD ρ≈0.88 – 0.97 (shaded band) is emergent without an explicit constraint. FNO (no guarantee) relies on emergent stability; WNO diverges (observed).
Method
Parameters
Relative to MPNO
LS-DYNA
0 (physics-based, no training)
N/A
FNO
∼ 79,000
2D implementation for evaluation
WNO
∼ 11,000
2D implementation for evaluation
MeshGraphNets
∼ 1,000,000
typical configuration
MPNO (hidden=32)
∼ 20,000
the proposed model
Table 2: Trainable parameter scale of each method
Model
Parameters
Energy retention after 30 steps
1st-step Rel-L2
Final-step Rel-L2
MPNO
4,482
Stable (energy@final=0.21)
0.37
1.49
FNO
57,825
0.32
0.22
1.28
WNO
17,121
1.64 (seed-dependent)
0.43
4.25
Table 3: Autoregressive-stability comparison on the Burgers equation
Figure 5: Cross-PDE autoregressive stability. (a) Burgers equation (wave/shock, 30-step autoregression): WNO diverges (0.43 → 4.25, red dashed) while MPNO (blue) and FNO (orange) remain bounded (0.37 → 1.49 and 0.22 → 1.28). (b) The same MPNO architecture is applied across PDEs by replacing only the material-property variables in the edge-weight formula (Burgers: uniform weight; Darcy: permeability; penetration: acoustic impedance Zm=1.0 / Za=1.495 ), with autoregressive stability maintained in every scenario. Single-step accuracies are reported in Table 3 (Burgers) and in the text (Darcy).
Figure 6: Spectral validation of the propagation operator on the cropped impact-zone graph (10 random aggregate distributions, K=20 power-iteration steps). (a) Relative error of the Rayleigh-quotient power-iteration estimate of λmax decays with K ; average underestimation 0.0005% at K=20 (max 0.003%). (b) Eigenvalue spectrum of P=I−αL~ ( α=0.577 ) lies in [0.423,1.000]⊂[0,1] (10/10 PSD), confirming ρ(P)≤1 by construction.
Note: MPNO–F13 is an improved variant of MPNO that recomputes the edge weights with an edge-MLP at every graph-propagation round (§ 4.3.1 ), letting the physics coupling run through the propagation process rather than being frozen once; with about 37K parameters (about half of FNO) it raises the single-step accuracy to 0.7278, narrowing the gap to FNO to ±0.007 , while remaining stable long-horizon.
Table 4: Single-step accuracy and autoregressive rollout comparison on the penetration data (single step with N=3 ; 29-step rollout = the entire active horizon from data frame 10)
Figure 7: (a) LS-DYNA σzz stress fields of the impact zone at frames t=10,20,30,39 (seed 205, 165 m/s), showing the penetration-induced stress wave propagating from the impact point; amplitude decays with geometric spreading. (b) Peak ∣σzz∣ in the impact zone over 40 frames: the stress peaks at t≈8 (10.82, normalized by 30 MPa) upon impact, then decays with oscillations due to wave propagation and reflection.
Figure 8: Pseudospectral transient growth of the full-model Jacobian (seed 55, frame 10, N=224 ). ∥Jt∥ (largest singular value of the t -step Jacobian power) is shown for MPNO (solid blue) and no_spec (dashed red); the horizontal dashed line marks the stability threshold ∥Jt∥=1 . MPNO’s transient peaks at t=1 (1.58, equal to σmax(J)=1.58 ) then collapses below 1 at t=2 and decays monotonically to 0.263, reconciling σmax(J)>1 with asymptotic contraction at this operating point ( ρ(J)=0.83≤1 ). no_spec retains growth above 1 through t=6 ( ρ(J)=0.90 ). Model spectral radii ρ(J) are positioned on the ρ -axis taxonomy in Fig. 4 .
Figure 9: Controlled counterfactual: inflating the learned λmax . The estimated λmax is multiplied by 1×/6×/8×/10× before re-normalizing the propagation operator (all else unchanged). MPNO (blue) preserves its stability plateau at every amplification factor (0/6 seeds diverge); no_spec (red) loses emergent stability monotonically, diverging in 2/6, 3/6, and 4/6 seeds at 6× , 8× , and 10× . The divergence of no_spec under spectrum inflation, despite identical stability on nominal data, shows that spectral normalization, not incidental training, underlies MPNO’s stability.
Configuration
Removed component
Three-case single-step Rel-L2
All-case single-step Rel-L2
Conclusion
baseline
(full model)
0.7304±0.0008
0.7108
Baseline ( N=3 )
no_spec
Spectral normalization
0.7366
0.7135
Consistent degradation across conventions
no_phys
Physics-coupled edge weight
0.7291
0.7256
Three-case neutral; benefits high-velocity cases
no_round_mlp
round_mlps cyclic post-processing
0.7314
0.7201
Three-case neutral; useful at high velocity
no_speed
Velocity coupling
–
0.7175
Useful
no_appnp
APPNP residual connection
–
0.7191
Useful
Table 5: Ablation study (single-step TEST Rel-L2; main baseline is the N=3 mean, ablations are single runs with SEED=42)
Figure 10: Efficiency and stability relative to baselines. Horizontal axis: inference time per 29-step rollout (log scale); vertical axis: stability type (constructed / emergent / diverging). MPNO is the only model combining a constructed stability guarantee with lightweight inference ( ∼ 0.1 s on CPU vs. 4–8 h for the LS-DYNA reference, a ≈ 10 5× speed-up); FNO is stable but without a guarantee; WNO diverges; POD-DMD is emergent. Parameter counts are listed in Table 2 .
Model
#Params
Single-step rel.-L2 ( ×1 )
Long-horizon completed
Diverged
Final-step error
MPNO (lightweight)
19,973
8.82±3.05
91/91
No
17.96
FNO
1,187,556
6.67±2.73
57/91
Yes
>50
WNO
73,956
6.68±4.29
4/91
Yes
>50
MeshGraphNet
63,780
11.79±3.49
91/91
No
22.74
Table 6: Cross-domain out-of-distribution generalization on the PDEBench 2D compressible-Euler Sod shock-tube (single trajectory, 101 frames; all models share one protocol). rel.-L2 = dimensionless relative L2 ratio ∥u^−u∥F/∥u∥F ; “diverged” = the error exceeds 50 or NaN before the rollout ends.
Predicting interfacial degradation and microstructural evolution through phase-field modelling is computationally intensive, particularly for parametric studies. While neural operators such as the Fourier neural operator (FNO) offer a promising route to acceleration, the absence of physical constraints compromises their generalisation and long-term stability. Here, we present PF-PINO, a physics-informed neural operator framework that embeds phase-field governing equation residuals directly into the training loss, together with gradient normalisation loss balancing and a staggered training scheme for coupled multi-field systems. PF-PINO is validated on four benchmarks including pencil-electrode corrosion, electro-polishing corrosion, dendritic solidification, and spinodal decomposition, collectively probing parametric generalisation across orders-of-magnitude kinetic regimes, non-periodic boundaries with stochastic morphologies, coupled multi-field dynamics with out-of-distribution extrapolation, and long-term autoregressive stability under complex pattern formation. Across all benchmarks, PF-PINO reduces relative L2 errors by up to an order of magnitude over FNO and achieves sub-mesh-size interface accuracy, establishing a data-efficient surrogate for high-throughput parametric phase-field assessment in computational materials science and engineering applications.
Nanxi Chen, Airong Chen, Rujin Ma
College of Civil Engineering, Tongji University, Shanghai 200092, China
Training neural PDE solvers is often bottlenecked by expensive data generation or unstable physics-informed neural network (PINN) involving challenging optimization landscapes due to higher-order derivatives. To tackle this issue, we propose an alternative approach using Monte Carlo approaches to estimate the solution to the PDE as a stochastic process for weak supervision during training. Leveraging the Walk-on-Spheres method, we introduce a learning scheme called Walk-on-Spheres Neural Operator (WoS-NO) which uses weak supervision from WoS to train any given neural operator. We propose to amortize the cost of Monte Carlo walks across the distribution of PDE instances using stochastic representations from the WoS algorithm to generate cheap, noisy, estimates of the PDE solution during training. This is formulated into a data-free physics-informed objective where a neural operator is trained to regress against these weak supervisions, allowing the operator to learn a generalized solution map for an entire family of PDEs. This strategy does not require expensive pre-computed datasets, avoids computing higher-order derivatives for loss functions that are memory-intensive and unstable, and demonstrates zero-shot generalization to novel PDE parameters and domains. Experiments show that for the same number of training steps, our method exhibits up to 8.75× improvement in L2-error compared to standard physics-informed training schemes, up to 6.31× improvement in training speed, and reductions of up to 2.97× in GPU memory consumption. We present the code at https://github.com/neuraloperator/WoS-NO
Hrishikesh Viswanath, Hong Chul Nam, Xi Deng +3
Purdue University, USA · ETH Zurich, Switzerland · California Institute of Technology, USA +1
Neural operators excel as deterministic surrogates, but inevitably collapse to the conditional mean when applied to stochastic PDEs, discarding the variance and tail structure upon which uncertainty quantification depends. Recovering this structure typically requires Monte Carlo rollouts or grafted generative models, both of which surrender the one-shot efficiency and resolution invariance that define the operator paradigm. To resolve this, we draw on the Doob-Meyer theorem, which establishes that any semimartingale fundamentally decomposes into a predictable drift and an unpredictable, zero-mean martingale. Translating this theorem into an architectural prior, we introduce the Martingale Neural Operator (MNO). MNO maps an initial condition directly to the conditional mean and covariance of the terminal law, parameterized by a drift-like mean and a low-rank factor Bφ with Bφ⊤Bφ positive semi-definite by construction. For our experiments, we use a Gaussian residual instantiation. Across 1D SPDEs, rough volatility, and 2D operator tasks, MNO reduces Wasserstein distance by up to 120× on φ4 field theory and 68× on stochastic Burgers, evaluating ∼3× faster than a conditional diffusion baseline at matched wall-clock training budgets. On 2D tasks, MNO is comparable to FNO on zero-shot resolution transfer and turbulent flow, while quasi-deterministic systems such as Gray-Scott remain a failure mode.
Kai Hidajat
Department of Applied Mathematics University of Washington Seattle, WA, USA