cs.LGSep 29, 2026

Neural Constitutive Learning for Generalized Reaction-Diffusion Systems

Authors: Shang-Ke Chen, Yu-Peng Wang, Shih-Hsuan Hung, Wei-Fang Sun, Chao-Shun Zhan, Simon See, Min-Jhe Lu

Organizations: National Tsing Hua University · NVIDIA

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Generalized reaction-diffusion systems encompass diverse transport mechanisms and coupled reaction kinetics. A central question for neural PDE solvers is what should be learned so that a common interface can accommodate phase-field and degenerate transport, local reactions, and multispecies coupling. We propose the Neural Constitutive Laws--Mass-Compression-Transport (NCL-MCT) Solver, which learns PDE-specific constitutive responses while retaining temporal evolution in a shared MCT integrator. Transport is represented through mobility and thermodynamic driving force, and reaction through relative reaction rates. These constitutive responses depend on the current density rather than explicitly on the initial condition or elapsed time, motivating their reuse across different initial conditions and time horizons. The same interface supports velocity-data supervision and known-law supervision, neither of which requires time integration during training. When constitutive laws are known, supervision can be evaluated on independently sampled density fields, enabling trajectory-free constitutive learning without generating solution trajectories. Across seven systems, separately trained constitutive modules share the same interface and MCT integrator and achieve relative rollout L2L^2 errors of 10−410^{-4} to 10−210^{-2}. Tests with unseen initial-condition families and an extended time horizon assess reuse beyond training conditions, while separate experiments demonstrate trajectory-free constitutive learning. These results support constitutive responses as an effective learning target for a shared neural PDE framework.

Figures & tables

Figure 1: Overview of NCL-MCT Solver. (a) At inference, frozen constitutive modules are coupled to the shared MCT integrator. At each step, the transport velocity u=ξf\mathbf{u}=\xi\mathbf{f} and relative reaction rate rr are reevaluated from the current density to evolve II and MM and reconstruct ρ=MI\rho=MI . During training, transport is supervised either (b) from reference velocities on PDE snapshots, with curl regularization, or (c) directly from known mobility and force laws on density fields, including independently sampled density fields. For reactive systems, rr is supervised by the known reaction law in both settings.
SystemNTN_{T}MetricF-FNO ( Tran et al., 2023 )CFO ( Hou et al., 2026 )DOOL ( Chang et al., 2026 )NCL-MCT (Vel.) (Ours)NCL-MCT (Law) (Ours)
(a) Generalized diffusion
Linear Diffusion600600ErollE_{\mathrm{roll}}(6.52±1.65)×10−4(6.52\pm 1.65)\times 10^{-4}(9.90±4.91)×10−3(9.90\pm 4.91)\times 10^{-3}(5.65±1.48)×10−2(5.65\pm 1.48)\times 10^{-2}(9.87±2.61)×10−4(9.87\pm 2.61)\times 10^{-4}(1.57±0.46)×10−3(1.57\pm 0.46)\times 10^{-3}
EmaxE_{\mathrm{max}}(9.76±2.78)×10−4(9.76\pm 2.78)\times 10^{-4}(1.19±0.60)×10−2(1.19\pm 0.60)\times 10^{-2}(7.36±1.48)×10−2(7.36\pm 1.48)\times 10^{-2}(1.18±0.29)×10−3(1.18\pm 0.29)\times 10^{-3}(1.82±0.59)×10−3(1.82\pm 0.59)\times 10^{-3}
CH10,00010{,}000ErollE_{\mathrm{roll}}(1.45±0.09)×10−1(1.45\pm 0.09)\times 10^{-1}(6.16±1.47)×10−2(6.16\pm 1.47)\times 10^{-2}(3.32±0.53)×10−1(3.32\pm 0.53)\times 10^{-1}(7.50±1.92)×10−3(7.50\pm 1.92)\times 10^{-3}(7.57±1.28)×10−3(7.57\pm 1.28)\times 10^{-3}
EmaxE_{\mathrm{max}}(1.87±0.14)×10−1(1.87\pm 0.14)\times 10^{-1}(8.78±2.35)×10−2(8.78\pm 2.35)\times 10^{-2}(4.71±0.74)×10−1(4.71\pm 0.74)\times 10^{-1}(1.04±0.30)×10−2(1.04\pm 0.30)\times 10^{-2}(1.05±0.25)×10−2(1.05\pm 0.25)\times 10^{-2}
PM2,0002{,}000ErollE_{\mathrm{roll}}(7.87±4.00)×10−2(7.87\pm 4.00)\times 10^{-2}(2.68±0.61)×10−3(2.68\pm 0.61)\times 10^{-3}(7.30±2.61)×10−1(7.30\pm 2.61)\times 10^{-1}(1.67±0.17)×10−2(1.67\pm 0.17)\times 10^{-2}(1.66±0.25)×10−2(1.66\pm 0.25)\times 10^{-2}
Table 1: Applicability across generalized reaction-diffusion systems. ErollE_{\mathrm{roll}} and EmaxE_{\mathrm{max}} denote the space-time and maximum per-time relative L2L^{2} errors, reported as mean ±\pm standard deviation over ten test trajectories; NTN_{T} is the rollout length. Vel./Law denote velocity-data and known-law supervision, respectively. Red/orange mark the best/second-best mean per row. Gray ∗ denotes F-FNO results computed from the finite trajectories only; −- denotes not evaluated.
Figure 2: Schnakenberg comparison for both species. (a) Concentrations UU (top) and VV (bottom) at t=1t=1 for a test trajectory on which all methods remain finite. (b) Per-time relative L2L^{2} errors over the 20,00020{,}000 -step rollout.
Figure 3: Generalization beyond training conditions on reactive Cahn–Hilliard. (a) From an unseen two-nucleus initial condition, F-FNO distorts the large domains and loses secondary droplets, while CFO merges them into elongated structures; both NCL-MCT Solver configurations better preserve the reference morphology. (b) Models trained on the first 60%60\% of the interval are rolled out over the full horizon. Beyond the training window, CFO loses interior contrast and smooths narrow channels, while both NCL-MCT Solver configurations better preserve the interfaces.
 MetricF-FNO ( Tran et al., 2023 )CFO ( Hou et al., 2026 )NCL-MCT (Vel.) (Ours)NCL-MCT (Law) (Ours)
Unseen initial conditionsErollE_{\mathrm{roll}}1.81×10−11.81\times 10^{-1} ∗1.69×10−11.69\times 10^{-1}2.99×10−22.99\times 10^{-2}2.86×10−22.86\times 10^{-2}
(5 systems)EmaxE_{\mathrm{max}}2.97×10−12.97\times 10^{-1} ∗2.38×10−12.38\times 10^{-1}4.39×10−24.39\times 10^{-2}4.44×10−24.44\times 10^{-2}
Extended time horizonErollE_{\mathrm{roll}}6.75×10−26.75\times 10^{-2} †6.89×10−26.89\times 10^{-2}5.85×10−35.85\times 10^{-3}5.49×10−35.49\times 10^{-3}
(4 systems)ΔE\Delta E1.25×10−21.25\times 10^{-2} †4.82×10−34.82\times 10^{-3}1.97×10−31.97\times 10^{-3}2.76×10−32.76\times 10^{-3}
∗ F-FNO is nonfinite on CH and Fisher-KPP; the reported average uses only the 3/5 finite systems and is not directly comparable.
† For reactive CH, F-FNO statistics use the 9/10 trajectories that remain finite.
Table 2: Generalization beyond training conditions, aggregated across systems: (a) unseen initial conditions (five systems, one test trajectory per system) and (b) extended time horizon (four systems, ten test trajectories per system). Entries are averages of the per-system ErollE_{\mathrm{roll}} , EmaxE_{\mathrm{max}} , or ΔE\Delta E values. Per-system results are reported in Appendix G . DOOL is omitted because it covers too few systems for aggregation. Red/orange indicate the best/second-best comparable result per row.
SystemNTN_{T}MetricDOOL ( Chang et al., 2026 )NCL-MCT (Law) (Ours)
Linear Diffusion4,0004{,}000ErollE_{\mathrm{roll}}(1.53±1.07)×10−3(1.53\pm 1.07)\times 10^{-3}(3.13±1.69)×10−4(3.13\pm 1.69)\times 10^{-4}
(1D, DOOL setting)EmaxE_{\mathrm{max}}(2.39±1.67)×10−3(2.39\pm 1.67)\times 10^{-3}(4.05±2.32)×10−4(4.05\pm 2.32)\times 10^{-4}
CH2,0002{,}000ErollE_{\mathrm{roll}}(7.46±2.90)×10−3(7.46\pm 2.90)\times 10^{-3}(6.40±6.66)×10−3(6.40\pm 6.66)\times 10^{-3}
(DOOL setting)EmaxE_{\mathrm{max}}(1.17±0.37)×10−2(1.17\pm 0.37)\times 10^{-2}(8.71±8.33)×10−3(8.71\pm 8.33)\times 10^{-3}
Table 3: Trajectory-free learning from independently sampled density fields. Both methods use the same unlabeled density fields and are evaluated on ten held-out initial conditions. ErollE_{\mathrm{roll}} and EmaxE_{\mathrm{max}} are relative L2L^{2} errors (mean ±\pm standard deviation), and NTN_{T} is the rollout length. Law denotes known-law supervision on independently sampled states. Red marks the lower mean per row.
Appendix figures & tables32 assets

Supplementary material from the paper’s appendix.

Appendix

SystemDomainOrderReactiontt rangeΔt\Delta tNTN_{T}
Linear Diffusion[0,1)2[0,1)^{2}2None[0,0.03][0,0.03]5×10−55\times 10^{-5}600
CH[0,1)2[0,1)^{2}4None[0,0.1][0,0.1]10−510^{-5}10410^{4}
PM[−4,4)2[-4,4)^{2}2 †None[0.1,0.3][0.1,0.3]10−410^{-4}2000
Fisher-KPP[0,1)2[0,1)^{2}2Logistic[0,0.015][0,0.015]3×10−53\times 10^{-5}500
Reactive CH[0,1)2[0,1)^{2}4Double-well[0,0.1][0,0.1]10−510^{-5}10410^{4}
Reactive PM[−4,4)2[-4,4)^{2}2 †Logistic[0.1,0.6][0.1,0.6]5×10−45\times 10^{-4}1000
Appendix
Table 4: Per-system configurations. “Order” denotes the highest spatial derivative order in the transport term of the density equation, and “Reaction” the local reaction source type. Δt\Delta t and NTN_{T} denote the reference time-step size and number of time steps, respectively. † Degenerate transport, for which the flux vanishes as ρ→0\rho\to 0 at the free boundary.
SystemFamilyRange or random parameters
Linear DiffusionGFRF, Equation 37ρ0∈[0.1,1.0]\rho_{0}\in[0.1,1.0]
CHGFRF, Equation 37ρ0∈[0.7,1.3]\rho_{0}\in[0.7,1.3]
Fisher-KPPGFRF, Equation 37ρ0∈[0.2,0.8]\rho_{0}\in[0.2,0.8]
Reactive CHGFRF, Equation 37ρ0∈[0.7,1.3]\rho_{0}\in[0.7,1.3]
PMBarenblatt, Equation 38R∼U[0.8,1.4]R\sim\mathcal{U}[0.8,1.4]
Reactive PMRadial bump, Equation 39R∼U[0.8,2.0]R\sim\mathcal{U}[0.8,2.0] , A∼U[0.5,0.9]A\sim\mathcal{U}[0.5,0.9]
Appendix
Table 5: Initial-condition family of each system. GFRF denotes Equation 37 rescaled to the stated density range. Random parameters are sampled independently for each trajectory.
Figure 4: Linear diffusion, 600-step rollout to T=0.03T=0.03 .
Figure 5: Cahn-Hilliard, 10,000-step rollout to T=0.1T=0.1 .
Figure 6: Porous medium, 2,000-step rollout from t=0.1t=0.1 to T=0.3T=0.3 .
Figure 7: Fisher-KPP, 500-step rollout to T=0.015T=0.015 .
Figure 8: Reactive Cahn-Hilliard, 10,000-step rollout to T=0.1T=0.1 .
Figure 9: Reactive porous medium, 1,000-step rollout from t=0.1t=0.1 to T=0.6T=0.6 .
Figure 10: Schnakenberg, activator UU , 20,000-step rollout to T=1T=1 .
Figure 11: Schnakenberg, inhibitor VV , 20,000-step rollout to T=1T=1 .
SystemFamilyNTN_{T}
Linear DiffusionSix Gaussian blobs on a circle, Equation 41600
CHTwo Gaussian blobs, Equation 4210410^{4}
Fisher-KPPGaussian ring, Equation 43500
Reactive CHTwo Gaussian blobs, Equation 4210410^{4}
SchnakenbergSpot lattice, Equation 442×1042\times 10^{4}
Appendix
Table 6: Out-of-distribution initial-condition families. Distances are periodic on each system’s domain, and all fields are evaluated on the same 128×128128\times 128 grid and rolled out over the same interval as the corresponding in-distribution test.
SystemNTN_{T}MetricF-FNO ( Tran et al., 2023 )CFO ( Hou et al., 2026 )DOOL ( Chang et al., 2026 )NCL-MCT (Vel.) (Ours)NCL-MCT (Law) (Ours)
Linear Diffusion600600ErollE_{\mathrm{roll}}1.71×10−11.71\times 10^{-1}1.65×10−11.65\times 10^{-1}5.74×10−15.74\times 10^{-1}4.07×10−24.07\times 10^{-2}3.81×10−23.81\times 10^{-2}
EmaxE_{\mathrm{max}}2.87×10−12.87\times 10^{-1}2.13×10−12.13\times 10^{-1}7.01×10−17.01\times 10^{-1}4.67×10−24.67\times 10^{-2}4.44×10−24.44\times 10^{-2}
CH10,00010{,}000ErollE_{\mathrm{roll}}NaN2.31×10−12.31\times 10^{-1}4.28×10−14.28\times 10^{-1}1.96×10−21.96\times 10^{-2}1.83×10−21.83\times 10^{-2}
EmaxE_{\mathrm{max}}NaN3.30×10−13.30\times 10^{-1}5.33×10−15.33\times 10^{-1}3.63×10−23.63\times 10^{-2}3.32×10−23.32\times 10^{-2}
Fisher-KPP500500ErollE_{\mathrm{roll}}NaN1.35×10−11.35\times 10^{-1}−-5.71×10−25.71\times 10^{-2}5.36×10−25.36\times 10^{-2}
EmaxE_{\mathrm{max}}NaN1.71×10−11.71\times 10^{-1}−-7.02×10−27.02\times 10^{-2}6.57×10−26.57\times 10^{-2}
Appendix
Table 7: Rollout errors from unseen initial conditions. Each system uses one fixed out-of-distribution initial condition, so entries are single rollouts without standard deviations. ErollE_{\mathrm{roll}} and EmaxE_{\mathrm{max}} denote space-time and maximum per-time relative L2L^{2} errors. Red/orange mark the best/second-best result per row. NaN denotes a nonfinite rollout; −- denotes a method not evaluated under the adopted formulation. Schnakenberg UU and VV count as one system.
Figure 12: Linear diffusion from six Gaussian blobs on a circle, Equation 41 , 600-step rollout to T=0.03T=0.03 . F-FNO loses the six-fold symmetry early and replaces the ring by two smeared diagonal lobes, while CFO retains separated blobs after the reference has largely merged them. DOOL shows distinct blobs on a striped background where the reference is nearly uniform. Both NCL-MCT Solver configurations more closely follow the merging and flattening dynamics.
Figure 13: Cahn-Hilliard from two Gaussian blobs, Equation 42 , 10,000-step rollout to T=0.1T=0.1 . F-FNO develops oscillations near the domain boundary at t=0.01t=0.01 and becomes nonfinite thereafter. CFO recovers the two large domains but forms an irregular set of secondary droplets, while DOOL produces a different labyrinthine pattern. Both NCL-MCT Solver configurations reproduce the two domains and the regular array of smaller droplets between them.
Figure 14: Fisher-KPP from a Gaussian ring, Equation 43 , 500-step rollout to T=0.015T=0.015 . The blank panel denotes a nonfinite state. F-FNO breaks the ring into two bright lobes and develops diagonal striping before becoming nonfinite, while CFO retains a pronounced central depression after the reference has filled in. Both NCL-MCT Solver configurations more closely follow the filling of the ring.
Figure 15: Reactive Cahn-Hilliard from two Gaussian blobs, Equation 42 , 10,000-step rollout to T=0.1T=0.1 . F-FNO preserves the large domains but distorts their shape and misses the secondary-droplet array, while CFO merges the droplets into elongated structures. Both NCL-MCT Solver configurations recover the array more closely.
Figure 16: Schnakenberg from a spot lattice, Equation 44 , activator UU , 20,000-step rollout to T=1T=1 . All methods remain finite but differ in pattern selection: F-FNO retains smooth diagonal bands where the reference has formed a spot lattice, while CFO and both NCL-MCT Solver configurations resolve the spots more closely.
Figure 17: Schnakenberg from a spot lattice, inhibitor VV , with the same rollout and layout as Figure 16 . The F-FNO banding is more visible for this species.
SystemΔt\Delta tNTtrainN_{T}^{\mathrm{train}} ( TtrainT_{\mathrm{train}} )NTN_{T} ( TT )Extrapolated fraction
Linear diffusion5×10−55\times 10^{-5}360360 ( 0.0180.018 )600600 ( 0.030.03 )40%40\%
CH10−510^{-5}6,0006{,}000 ( 0.060.06 )10,00010{,}000 ( 0.10.1 )40%40\%
Fisher-KPP3×10−53\times 10^{-5}300300 ( 0.0090.009 )500500 ( 0.0150.015 )40%40\%
Reactive CH10−510^{-5}6,0006{,}000 ( 0.060.06 )10,00010{,}000 ( 0.10.1 )40%40\%
Appendix
Table 8: Training and evaluation horizons. NTtrainN_{T}^{\mathrm{train}} and NTN_{T} denote the rollout lengths of the truncated training trajectories and full evaluation, respectively. The step size is unchanged and equals the reference Δt\Delta t in Table 4 .
SystemMetricF-FNO ( Tran et al., 2023 )CFO ( Hou et al., 2026 )NCL-MCT (Vel.) (Ours)NCL-MCT (Law) (Ours)
Linear DiffusionErollE_{\mathrm{roll}}(4.01±1.26)×10−4(4.01\pm 1.26)\times 10^{-4}(5.13±1.04)×10−2(5.13\pm 1.04)\times 10^{-2}(1.44±0.27)×10−3(1.44\pm 0.27)\times 10^{-3}(2.26±0.37)×10−3(2.26\pm 0.37)\times 10^{-3}
EmaxE_{\mathrm{max}}(6.09±1.87)×10−4(6.09\pm 1.87)\times 10^{-4}(6.25±1.52)×10−2(6.25\pm 1.52)\times 10^{-2}(1.62±0.34)×10−3(1.62\pm 0.34)\times 10^{-3}(2.60±0.46)×10−3(2.60\pm 0.46)\times 10^{-3}
Eˉ(Ttrain)\bar{E}(T_{\mathrm{train}})3.78×10−43.78\times 10^{-4}5.82×10−25.82\times 10^{-2}1.58×10−31.58\times 10^{-3}2.53×10−32.53\times 10^{-3}
ΔE\Delta E1.69×10−41.69\times 10^{-4}2.09×10−32.09\times 10^{-3}−5.72×10−5-5.72\times 10^{-5}−1.75×10−4-1.75\times 10^{-4}
CHErollE_{\mathrm{roll}}(1.39±0.17)×10−1(1.39\pm 0.17)\times 10^{-1}(9.96±0.86)×10−2(9.96\pm 0.86)\times 10^{-2}(9.27±3.03)×10−3(9.27\pm 3.03)\times 10^{-3}(1.03±0.27)×10−2(1.03\pm 0.27)\times 10^{-2}
EmaxE_{\mathrm{max}}(1.85±0.24)×10−1(1.85\pm 0.24)\times 10^{-1}(1.26±0.13)×10−1(1.26\pm 0.13)\times 10^{-1}(1.44±0.89)×10−2(1.44\pm 0.89)\times 10^{-2}(1.49±0.37)×10−2(1.49\pm 0.37)\times 10^{-2}
Appendix
Table 9: Rollout errors with the horizon extended beyond training. Models are trained on the first 60%60\% of each reference interval and evaluated over the full interval (Table 8 ). ErollE_{\mathrm{roll}} and EmaxE_{\mathrm{max}} are mean ±\pm standard deviation over ten test trajectories. Eˉ(Ttrain)\bar{E}(T_{\mathrm{train}}) is the trajectory-averaged per-time error at the training horizon, and ΔE\Delta E is its change to the final time, Equation 45 . Red/orange mark the best/second-best mean per row.
Figure 18: Linear diffusion, trained to Ttrain=0.018T_{\mathrm{train}}=0.018 and rolled out to T=0.03T=0.03 . Only the first row lies within the training horizon.
Figure 19: Cahn-Hilliard, trained to Ttrain=0.06T_{\mathrm{train}}=0.06 and rolled out to T=0.1T=0.1 .
Figure 20: Fisher-KPP, trained to Ttrain=0.009T_{\mathrm{train}}=0.009 and rolled out to T=0.015T=0.015 .
Figure 21: Reactive Cahn-Hilliard, trained to Ttrain=0.06T_{\mathrm{train}}=0.06 and rolled out to T=0.1T=0.1 .
Figure 22: One-dimensional linear diffusion in the DOOL setting after trajectory-free constitutive learning. Space-time density for the visualized initial condition ρ0(x)=2+sin⁡x\rho_{0}(x)=2+\sin x , rolled out to T=1T=1 in 4,0004{,}000 steps. Columns show the analytic solution, NCL-MCT Solver (Law), and DOOL using a shared color scale.
Figure 23: Two-dimensional Cahn-Hilliard in the DOOL setting after trajectory-free constitutive learning, shown using the shifted density ρ=u+2\rho=u+2 for u0=sin⁡xsin⁡yu_{0}=\sin x\sin y . The rollout extends to T=0.2T=0.2 over 2,0002{,}000 steps. Rows show t=0.01t=0.01 , 0.100.10 , and 0.200.20 , with an independent color scale for each row. Both methods preserve the checkerboard structure, while at the final time DOOL exhibits grid-aligned texture and a distorted interface; NCL-MCT Solver (Law) remains smooth with slightly weaker interior contrast.
SystemMethodNTN_{T}ParamsVRAM (GiB)Train (s)Infer (s)ms/step
Linear DiffusionF-FNO600600733,506733{,}5068.3718.37111,11711{,}1170.80±0.000.80\pm 0.001.341.34
CFO8,095,4898{,}095{,}4896.9706.9705,7125{,}7124.91±0.014.91\pm 0.018.188.18
DOOL64,62064{,}6200.0820.0821321320.45±0.040.45\pm 0.040.750.75
NCL-MCT Solver (Vel.)97,41997{,}4191.0611.0611,9121{,}9121.22±0.001.22\pm 0.002.032.03
NCL-MCT Solver (Law)97,41997{,}4191.0491.0491,8661{,}8661.22±0.001.22\pm 0.002.032.03
CHF-FNO10,00010{,}000733,506733{,}5068.3718.37111,11711{,}11713.15±0.0113.15\pm 0.011.321.32
Appendix
Table 10: Computational cost of all methods. Params counts trainable parameters; peak VRAM is the framework-reported maximum live-tensor memory during training; Train is the wall-clock time for 100,000100{,}000 optimizer steps; Infer is the mean wall-clock time of one NTN_{T} -step rollout over the ten test trajectories; and ms/step is Infer /NT/N_{T} . −- denotes unevaluated cases.
Figure 24: Relative mass error for porous medium, measured against the total mass at the initial time. The dashed curve shows the corresponding drift of the reference evaluation.
SystemMetricHC-PINN ( Hao et al., 2026 )NCL-MCT (Vel.) (Ours)NCL-MCT (Law) (Ours)
(a) Generalized diffusion
Linear DiffusionErollE_{\mathrm{roll}}6.43×10−46.43\times 10^{-4}2.17×10−32.17\times 10^{-3}2.23×10−32.23\times 10^{-3}
EmaxE_{\mathrm{max}}1.55×10−31.55\times 10^{-3}3.78×10−33.78\times 10^{-3}4.18×10−34.18\times 10^{-3}
CHErollE_{\mathrm{roll}}2.44×10−12.44\times 10^{-1}3.75×10−23.75\times 10^{-2}6.89×10−36.89\times 10^{-3}
EmaxE_{\mathrm{max}}3.02×10−13.02\times 10^{-1}4.93×10−24.93\times 10^{-2}1.11×10−21.11\times 10^{-2}
PMErollE_{\mathrm{roll}}1.04×10−11.04\times 10^{-1}1.37×10−21.37\times 10^{-2}1.38×10−21.38\times 10^{-2}
Appendix
Table 11: Single-trajectory comparison with HC-PINN. HC-PINN and both NCL-MCT Solver configurations are trained and evaluated on the same initial-value problem; the NCL-MCT Solver results therefore differ from the ten-trajectory averages in Table 1 . Red/orange mark the best/second-best result per row.
  Schnakenberg ( UU )Schnakenberg ( VV )
ParameterizationObjectiveErollE_{\mathrm{roll}}EmaxE_{\mathrm{max}}ErollE_{\mathrm{roll}}EmaxE_{\mathrm{max}}
u,r\mathbf{u},\mathbf{r}Ludata+Lrlaw\mathcal{L}^{\mathrm{data}}_{u}+\mathcal{L}^{\mathrm{law}}_{r}(3.40±1.70)×10−2(3.40\pm 1.70)\times 10^{-2}(1.11±0.66)×10−1(1.11\pm 0.66)\times 10^{-1}(7.84±4.27)×10−3(7.84\pm 4.27)\times 10^{-3}(2.96±2.18)×10−2(2.96\pm 2.18)\times 10^{-2}
(ξ,f),r(\xi,\mathbf{f}),\mathbf{r}Lvel\mathcal{L}_{\mathrm{vel}}(2.52±1.77)×10−2(2.52\pm 1.77)\times 10^{-2}(8.28±8.55)×10−2(8.28\pm 8.55)\times 10^{-2}(6.74±3.37)×10−3(6.74\pm 3.37)\times 10^{-3}(2.19±1.93)×10−2(2.19\pm 1.93)\times 10^{-2}
(ξ,f),r(\xi,\mathbf{f}),\mathbf{r}Llaw\mathcal{L}_{\mathrm{law}}(2.38±1.26)×10−2(2.38\pm 1.26)\times 10^{-2}(7.36±6.15)×10−2(7.36\pm 6.15)\times 10^{-2}(6.57±2.27)×10−3(6.57\pm 2.27)\times 10^{-3}(1.97±1.30)×10−2(1.97\pm 1.30)\times 10^{-2}
Appendix
Table 12: Ablation of constitutive parameterization and supervision on Schnakenberg. Errors are mean ±\pm standard deviation. Red/orange indicate the best/second-best unrounded mean in each column.
TransportReactionSpeciesErollE_{\mathrm{roll}}EmaxE_{\mathrm{max}}
NumericalNumericalUU(6.57±2.43)×10−3(6.57\pm 2.43)\times 10^{-3}(1.634±0.480)×10−2(1.634\pm 0.480)\times 10^{-2}
VV(2.11±0.80)×10−3(2.11\pm 0.80)\times 10^{-3}(6.25±1.95)×10−3(6.25\pm 1.95)\times 10^{-3}
NeuralNumericalUU(1.059±0.259)×10−2(1.059\pm 0.259)\times 10^{-2}(2.361±0.412)×10−2(2.361\pm 0.412)\times 10^{-2}
VV(3.07±0.58)×10−3(3.07\pm 0.58)\times 10^{-3}(7.52±1.35)×10−3(7.52\pm 1.35)\times 10^{-3}
NumericalNeuralUU(2.191±1.507)×10−2(2.191\pm 1.507)\times 10^{-2}(7.182±6.851)×10−2(7.182\pm 6.851)\times 10^{-2}
VV(6.63±3.40)×10−3(6.63\pm 3.40)\times 10^{-3}(2.178±1.809)×10−2(2.178\pm 1.809)\times 10^{-2}
Appendix
Table 13: Ablation of constitutive modularity on Schnakenberg. Numerical denotes constitutive responses computed from known laws using numerical discretization; Neural denotes responses predicted by trained constitutive modules. Errors are mean ±\pm standard deviation. Red/orange mark the best/second-best means for each species. Numerical-Numerical is the numerical constitutive reference within the same MCT integrator.
λcurl\lambda_{\mathrm{curl}}ErollE_{\mathrm{roll}}EmaxE_{\mathrm{max}}
00(1.03±0.23)×10−3(1.03\pm 0.23)\times 10^{-3}(1.32±0.36)×10−3(1.32\pm 0.36)\times 10^{-3}
0.010.01(5.01±1.39)×10−4(5.01\pm 1.39)\times 10^{-4}(5.96±1.96)×10−4(5.96\pm 1.96)\times 10^{-4}
0.10.1(8.61±1.59)×10−4(8.61\pm 1.59)\times 10^{-4}(1.17±0.26)×10−3(1.17\pm 0.26)\times 10^{-3}
11(1.13±0.19)×10−3(1.13\pm 0.19)\times 10^{-3}(1.72±0.53)×10−3(1.72\pm 0.53)\times 10^{-3}
Appendix
Table 14: Ablation of the curl-penalty weight λcurl\lambda_{\mathrm{curl}} on Fisher-KPP under velocity-data supervision. λcurl=0\lambda_{\mathrm{curl}}=0 removes the penalty. Errors are mean ±\pm standard deviation over ten test trajectories with 500500 rollout steps. Red/orange mark the best/second-best means. λcurl=0.01\lambda_{\mathrm{curl}}=0.01 is used throughout the paper.
More in this topic

Reaction-Diffusion Systems

See all 13 papers

Explore similar work

CardsList
Apr 2, 2026

Flow Learners for PDEs: Toward a Physics-to-Physics Paradigm for Scientific Computing

Partial differential equations (PDEs) govern nearly every physical process in science and engineering, but solving them at scale remains prohibitively expensive. Generative AI has transformed language, vision, and protein science, but learned PDE solvers have not undergone a comparable shift. Existing paradigms each capture part of the problem. Physics-informed neural networks embed residual structure, although they are often difficult to optimize in stiff, multiscale, or large-domain regimes. Neural operators amortize across instances, although they commonly inherit a snapshot-prediction view of solving and can degrade over long rollouts. Diffusion-based solvers model uncertainty, although they are often built on a solver template that still centers on state regression. We argue that the core issue is the abstraction used to train learned solvers. Many models are asked to predict states, while many scientific settings require modeling how uncertainty moves through constrained dynamics. The relevant object is transport over physically admissible futures. This motivates flow learners: models that parameterize transport vector fields and generate trajectories through integration, echoing the continuous dynamics that define PDE evolution. This physics-to-physics alignment supports continuous-time prediction, native uncertainty quantification, and new opportunities for physics-aware solver design. We explain why transport-based learning offers a stronger organizing principle for learned PDE solving and outline the research agenda that follows from this shift.

Yilong Dai, Shengyu Chen, Xiaowei Jia +1

The University of Alabama, Tuscaloosa, Alabama, USA · University of Pittsburgh, Pittsburgh, Pennsylvania, USA

Neural PDE SolversPDE SolvingFlow-Based Generative Modeling
AskPDF
May 16, 2026

Prediction Is Not Physics: Learning and Evaluating Conserved Quantities in Neural Simulators

A diffusion model trained on Hamiltonian trajectories can achieve rollout MSE near 10−310^{-3}, but the standard deviation of its energy over time is between 7500 and 36000 times larger than the ground-truth energy standard deviation, indicating a failure to preserve conservation laws. This gap motivates our central question of whether neural networks can learn or select globally conserved quantities from physical trajectories. We investigate this across three Hamiltonian systems: projectile motion, pendulum, and spring-mass. We use a structured T(v)+V(q)T(v)+V(q) energy model, a black-box Conservation Discovery Network (CDN), a polynomial CDN, and a conditional diffusion baseline. The structured network reaches R2≥0.9999R^2 \geq 0.9999 against analytical energy on clean data, while the black-box CDN reaches R2≥0.996R^2 \geq 0.996 when trained with temporal consistency plus a small alignment loss to analytical energy at t=0t=0 (λalign=0.2λ_{\mathrm{align}}=0.2). With λalign=0λ_{\mathrm{align}}=0, CDN Pearson R2R^2 collapses on pendulum and spring-mass (<10−3< 10^{-3}), showing that temporal consistency alone is not enough to reliably identify the true energy. Under 1%1\% additive Gaussian noise, the CDN outperforms the structured model on the projectile and spring-mass systems, suggesting that the CDN may be more robust to noisy inputs in this setting. However, the polynomial CDN is sensitive to training configuration: it achieves R2=0.78R^2=0.78 under a short training schedule on the pendulum system, but reaches R2=0.9998R^2=0.9998 with more training time and data, regardless of whether noise is added.

Andrew Bukowski, Aditya Kothari, Simba Shi +1

Yale University

Neural Surrogate ModelingDynamical SystemsScientific DiscoveryHamiltonian DynamicsDiffusion Models
AskPDF
Aug 31, 2026

Learning PDE Time-Stepping with Neural Cellular Automata

Classical numerical solvers for partial differential equations (PDEs) are computationally expensive to solve repeatedly across varying initial conditions, motivating the need for learned surrogates. In this paper, we propose a trainable Neural Cellular Automata (NCA) based surrogate model for learning long time PDE dynamics. Rather than mapping an entire initial field to a full trajectory in one shot, our proposed model learns a small, local, homogeneous update rule that is applied identically and repeatedly at every grid cell, mirroring the locality of differential operators. We benchmark this framework against three baselines: PDE - Net, a modified physics-informed neural network (PINN), and a Fourier Neural Operator (FNO), on five canonical PDEs (heat, advection, Burgers, Allen - Cahn, and Fisher - KPP), evaluated at temporal domain two times beyond the training temporal domain. The proposed model achieves the lowest long-horizon relative errors on the majority of the experiments.

Esha Saha, Hao Wang

1Interdisciplinary Lab for Mathematical Ecology and Epidemiology (ILMEE) & The Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton (AB), T6G 2J5, Canada

Neural Cellular AutomataPDE SolvingPDE Surrogate Modeling
AskPDF