HotPOT: Operator Transformer with Heterogeneous Physics Experts for Large-Scale PDE Pre-training

Introduction

Partial differential equations (PDEs) are fundamental to modeling a wide range of physical systems. The high cost of repeated numerical simulations has motivated learning-based PDE solvers, particularly neural operators ([11, 8]) that learn solution mappings directly from data. Recent methods have shown that a single model can be trained jointly on data from multiple PDEs, rather than trained separately for each task, and have achieved strong performance across multiple PDE tasks.

However, adding data from heterogeneous PDEs with different governing equations, physical variables, and variable interactions does not always improve performance on existing tasks. For example, the MoE-POT model proposed in Wang et al. [27] reports increased error on some original tasks when expanding the training mixture of a dense pre-trained model from six to twelve datasets. To mitigate cross-task interference, MoE-POT uses learned routing to select specialized experts. These experts nevertheless use homogeneous architectures, which may be difficult to directly learn inherent physical distinctions in state evolution across equations.

On the other hand, Neural–numerical methods offer a basis for constructing heterogeneous experts by combining learned predictions with explicit numerical operations ((Karlbauer et al., 2022; [3]. Different PDEs can require different update operations and prediction quantities. An expert may predict a rate for time integration, an increment for a constrained update, or state components for direct assignment. Pairing a predictor with the corresponding update allows the expert’s computation to reflect the requirements of its PDE. This motivates a neural–numerical pre-training paradigm in which each PDE pairs learnable predictions with explicit numerical operations according to its update requirements, while shared history representations integrate the resulting experts into a unified pre-training model.

Specifically, we introduce HotPOT, an operator transformer that combines a shared encoder with heterogeneous neural–numerical experts for large-scale PDE pre-training. Figure 1(a) summarizes the workflow from heterogeneous PDE pre-training to prediction and transfer. The shared encoder processes scalar-field histories across PDEs, while each expert handles system-specific variable coupling and state evolution. Each expert pairs a learnable predictor with a PDE-specific fixed numerical update, which contains no trainable parameters. The predictor outputs the quantities required by its paired numerical update, which then computes the next state. The shared encoder and expert predictors are jointly trained using next-state supervision, bringing the different neural–numerical experts together in a single pre-training model. The Navier–Stokes (NS) example in Figure 1(b) illustrates the resulting prediction: at the tenth forecast step, HotPOT follows the reference vorticity contours more closely than MoE-POT in the displayed central region.

The main contributions are:

• We propose HotPOT, an operator transformer that leverages PDE-specific experts containing learnable predictors with fixed numerical updates, allowing different equations to use different prediction quantities and update operations.

• We integrate these heterogeneous experts through a shared history encoder into one multi-PDE model. The encoder and predictors are trained end-to-end using only next-state supervision, while variable coupling and state evolution remain system-specific.

• We evaluate HotPOT across five PDE families and six settings. Relative to MoE-POT, HotPOT reduces rollout error by 50.7% on average across 24 matched-scale comparisons before and after fine-tuning, and further improves zero-shot prediction under unseen physical conditions.

Related Work

Multi-PDE pre-training. Previous studies including MPP [14], DPOT [5], Poseidon [7], and Walrus [15] learn reusable representations through joint training across PDE systems. Some related approaches propose to better organize representations around physical variables or scalar fields, as in CoDA-NO [18] and Tadpole [12], or condition prediction on equation and system information, as in PDEformer [29], PROSE-PDE [22], and Unisolver [32]. These studies provide the foundation for learning across heterogeneous PDE data. We focus on incorporating PDE-specific prediction and update operations into this shared pre-training framework.

Mixture-of-experts for PDEs. MoE-POT [27] combines shared experts with input-dependent routed experts, while NESTOR [21] routes computation at both image and token levels. These methods allocate specialized neural capacity to different inputs through learned routing, while retaining a common expert structure. HotPOT extends expert specialization to the quantities predicted and the operations used to evolve the state. Each PDE family uses a specific neural predictor paired with its corresponding fixed numerical update. A shared history encoder connects these heterogeneous experts within one jointly trained model.

Neural–numerical methods. Neural–numerical methods place learned components within structured updates. Universal differential equations (Rackauckas et al., 2021) and APHYNITY [30] augment known dynamics, while learned discretizations and solver-in-the-loop methods learn numerical operators or corrections [1, 10, 26]. FINN [9] embeds trainable components in PDE-specific finite-volume structures, and OpsSplit [3] combines learned nonlinear operators with fixed finite-difference operators and time integration. CompNO [8] combines pretrained operators with boundary enforcement, and ComPhy [25] couples equation-specific modules. Origo [23] decomposes evolution into global spectral and local constitutive operators, while Serrano et al. (2026) compose learned operators at test time. DISCO [16] generates a trajectory-conditioned operator for time integration. Building on this line of work, HotPOT organizes PDE-specific predictor–update pairs as heterogeneous experts for multi-PDE pre-training. Their prediction quantities and update operations differ across PDEs, while the shared encoder and expert predictors are jointly trained using next-state supervision.

Problem Setting

We consider joint pre-training over a collection of time-dependent PDE systems. Let $\mathcal{R}$ denote the set of prediction settings, and let $r \in\mathcal{R}$ index one setting. A setting specifies a physical system together with its benchmark configuration, such as the observed variables, equation parameters, geometry, or boundary conditions. Different settings may therefore have different state variables, variable couplings, and state-transition mechanisms.

Each setting $r$ has a trajectory dataset $\mathcal{D}_r$ with $C_r$ observed physical fields. Flattening the spatial grid in a fixed order gives $\mathbf{q}_{r,c}^{n} \in\mathbb{R}^{N_r}$ for field $c$ at saved step $n$, where $N_r$ is the number of grid cells and $q_{r,c,i}^{n}$ is its scalar value at cell $i$. The operator concat concatenates its entries in the specified order. The full state is the vector $\mathbf{q}_r^n = \operatorname{concat}(\{\mathbf{q}_{r,c}\}_{c=1}^{C_r}) \in\mathbb{R}^{N_r C_r}$. Geometry and support channels describe the spatial domain and valid grid locations. They are provided as auxiliary inputs rather than physical fields to be predicted, so they are excluded from $C_r$ and the state vector $\mathbf{q}_r^n$. The pre-training corpus is $\mathcal{D} = \bigcup_{r \in\mathcal{R}} \mathcal{D}_r$.

Let $g_{\mathbf{w}}$ denote the next-state predictor with trainable parameter vector $\mathbf{w}$. The operator stack stacks its indexed row vectors in order. Let $\widehat{\mathbf{q}}_r^{n+1}$ denote the predicted next state. The most recent $T = 10$ saved states form the history matrix $\mathbf{H}_r^n$:

Thus $\mathbf{H}_r^n \in\mathbb{R}^{T \times N_r C_r}$ contains states in chronological order. The matrix $\mathbf{H}_{r,c}^n \in\mathbb{R}^{T \times N_r}$ is its field-$c$ block. The argument $r$ identifies the PDE setting. The predictor also uses the geometry and support information, which we omit from the arguments of $g_{\mathbf{w}}$ to keep the notation compact. Each prediction is appended to the history for autoregressive rollout. We use bold lowercase symbols for vectors, bold uppercase symbols for matrices, and calligraphic uppercase symbols for sets or spaces. Indexed braces gather entries in index order. Parentheses group arguments or pairs.

Method

HotPOT builds its shared encoder and neural predictors using patch encoding, temporal aggregation, and grouped Fourier/local blocks from DPOT and MoE-POT ([5, 27]). As shown in Figure 2, HotPOT realizes $g_w$ by connecting PDE-family predictor–update pairs through a shared history encoder.

Let $\tau:\mathcal{R}\to\mathcal{T}$ map settings to the set of expert families. The family $\tau(r)$ is assigned to setting $r$. Settings from the same PDE family use the corresponding predictor and numerical update program. In particular, the two NS viscosity settings share one NS expert family. The shared encoder has parameter vector $\theta$. For family $e\in\mathcal{T}$, $\phi_e$ parameterizes field coupling and evolution, and $\psi_e$ parameterizes decoding and readout. Their aggregation is $w=\operatorname{concat}(\theta,\{(\phi_e,\psi_e)\}_{e\in\mathcal{T}})$, where parameter pairs are flattened and concatenated in a fixed family order.

Shared Field-Wise History Representation

The deterministic function $f_{\mathrm{feat},r,c}$ normalizes field $c$’s history and constructs six input features: normalized field value, two coordinates, boundary distance, geometry support, and field validity. This mapping does not mix physical fields. Each scalar history then passes through the same encoder $f_{\mathrm{enc},\theta}$, with no learned cross-field mixing. At saved time $n$, let $Z_{r,c}\in\mathbb{R}^{P\times d_{\mathrm{enc}}}$ denote its representation, where $P$ is the number of spatial patches and $d_{\mathrm{enc}}$ the shared width. Then

A shared patch encoder and positional embedding produce tokens for each history position. Temporal aggregation then combines these tokens within each field. Shared Fourier/local and channel feed-forward network (FFN) blocks produce this feature matrix. The saved-time index $n$ is suppressed. Next-state supervision jointly trains these features across datasets. Appendix § 9.§ 1 gives the aggregation formula, dimensions, and block equations.

System-Specific Prediction

Each PDE family’s predictor supplies the quantities required by its numerical update: transport rates for NS, conservative-variable actions for CNS, height increments for SWE, concentration-increment proposals for DR, or next-state velocities for CFDB. We denote this prediction space by $\mathcal{U}_r$.

Field coupling and evolution. Let $f_{\mathrm{mix},\phi_{\tau(r)}}$ denote field coupling followed by system-specific evolution. For field $c$, let $\mu_{r,c}$ and $\sigma_{r,c}$ be the history mean and floored standard deviation over space and the $T$ input frames, and collect them as $\mathbf{s}_{r,c} = (\mu_{r,c}, \log\sigma_{r,c})^\top\in\mathbb{R}^2$. The coupled context $\mathbf{Z}^{\mathrm{mix}}_r \in\mathbb{R}^{P \times d_{\mathrm{priv}}}$, with private width $d_{\mathrm{priv}}$, is

Coupling concatenates field representations, applies a $1 \times1$ projection, and adds history-statistic and positional embeddings. System-specific Fourier/local/FFN blocks then produce the context.

Spatial decoding and prediction. The decoder $f_{\mathrm{dec},\psi_{\tau(r)}}$ expands $\mathbf{Z}^{\mathrm{mix}}_r$ to full resolution using a transposed convolution, Gaussian error linear unit (GELU) activation, and pointwise projection. Its readout combines this cached context with a query state $\boldsymbol{\xi}^{n,k}_r \in\mathbb{R}^{d_{\xi,r}}$, representing the current or an intermediate state. Here $k = 1,\ldots,K_r$ indexes queries within the saved interval $\Delta t_r$, the time between consecutive saved states. Let $\widehat{\mathbf{u}}^{n,k}_r \in\mathcal{U}_r \subseteq\mathbb{R}^{d_{u,r}}$ denote the neural prediction at query $k$; $d_{\xi,r}$ and $d_{u,r}$ count the spatial and channel entries of the flattened query and prediction, respectively. Then

The same symbol $\widehat{\mathbf{u}}^{n,k}_r$ represents the different physical quantities listed above. The following section specifies their interpretation, reshaping, and conversion into the next state.

System-Specific Numerical Execution

Let $f^{\Delta t_r}_{\mathrm{upd},\tau(r)}$ be the PDE family’s fixed numerical update over $\Delta t_r$. It combines the current state $\mathbf{q}^n_r$ with the neural predictions from Eq. (4) to produce the next-state estimate $\widehat{\mathbf{q}}^{n+1}_r$:

The program has no trainable parameters. For multi-stage updates, it constructs successive query states from earlier predictions and reuses the history context and readout weights. Fixed input/output mappings connect its grid to the observed state layout; the PDE’s physical coefficients are implicit in its arguments. Below we give the core updates, omitting setting index $r$. Let $N$ be the program’s grid-cell count, $\mathbf{1}_N$ the all-ones vector, and $\langle\cdot\rangle$ spatial mean. A superscript $+$ denotes a program’s updated state. Operator details and update properties are in Appendices § 9.§ 2–§ 9.§ 3.

Incompressible Navier–Stokes (NS). Here the neural output $\widehat{\mathbf{u}}^{n,k}$ is an effective vorticity transport rate. Let $\boldsymbol{\omega} \in\mathbb{R}^N$ be the current substep vorticity, $\mathbf{f} \in\mathbb{R}^N$ the forcing, $\nu$ the viscosity, and $\delta t$ the substep duration. Write $f^s_{\mathrm{diff},\nu}$ for numerical diffusion over duration $s$. Let $\boldsymbol{\omega}_a$ be the diffused state and $\boldsymbol{\omega}_m$ the midpoint state; $\widehat{\mathbf{r}}_1$ and $\widehat{\mathbf{r}}_2$ are the corresponding rate outputs of Eq. (4) at these two states. The updated vorticity $\boldsymbol{\omega}^{+}$ is computed in this order:

Composing these substeps yields the next saved vorticity state. Transport rates are learned; diffusion and midpoint integration are fixed operations.

Compressible Navier–Stokes (CNS). Reshape the neural output $\widehat{\mathbf{u}}^{n,1}$ as a normalized action matrix $\mathbf{A} \in\mathbb{R}^{N\times4}$ for density, two momentum components, and total energy. Let $\mathbf{Q} \in\mathbb{R}^{N\times4}$ be the current conservative state, $\bar{\mathbf{a} } \in\mathbb{R}^4$ the action’s channel means, and $f_{\mathrm{scale}}(\mathbf{Q}) \in\mathbb{R}^4$ its conservative-field channel scales; $\operatorname{diag}$ forms a diagonal matrix. Denote the increment by $\Delta\mathbf{Q}$, the candidate by $\mathbf{Q}^{\star}$, and the updated state by $\mathbf{Q}^{+}$. A global coefficient $\alpha_{\mathrm{CNS}} \in[0,1]$ enforces positive density and pressure:

The program restores candidate conservative means before the segment update and converts $\mathbf{Q}^{+}$ back to primitive variables for the next state. For admissible inputs, global mass, momentum, and energy are preserved in exact arithmetic (Proposition 2).

Shallow-water equation (SWE). The neural output $\hat{\mathbf{u}}^{n,1} = \Delta\mathbf{h} \in\mathbb{R}^{N}$ is a height increment. Let $\mathbf{h}^{n}$ and $\hat{\mathbf{h}}^{n+1}$ be the current and predicted height vectors, $\mathbf{d}$ the centered increment, and $\alpha\in[0,1]$ a global coefficient enforcing a positive height floor. The update is

Centering and uniform scaling preserve mass on the uniform grid in exact arithmetic while respecting the height floor for admissible inputs (Proposition 1).

Diffusion–reaction (DR). The neural outputs $\delta_{1} = \hat{\mathbf{u}}^{n,1}$ and $\delta_{2} = \hat{\mathbf{u}}^{n,2}$ are concentration-increment proposals in $\mathbb{R}^{2N}$ for the two fields. The first is queried at the current state $\mathbf{q}^{n}$; the second at the intermediate state $\boldsymbol{\xi}^{n,2}$. The program constructs this state and then adds the averaged proposals to obtain $\hat{\mathbf{q}}^{n+1}$:

CFDBench (CFDB). For the velocity-prediction tasks in CFDB, reshape the neural output $\hat{\mathbf{u}}^{n,1}$ as $\hat{\mathbf{V}} \in\mathbb{R}^{N \times2}$, containing the two predicted velocity fields. The next velocity state $\mathbf{V}^{n+1}$ is assigned directly:

The program retains auxiliary support channels; output mapping restores the observed layout.

Joint Training and Autoregressive Prediction

Let $\boldsymbol{\pi}$ be the nonnegative setting-weight vector, with $\sum_{r\in\mathcal{R}}\pi_{r}=1$, and let $\mathcal{O}_{r}\subseteq\{1,\ldots,C_{r}\}$ be the set of scored physical-channel indices. We sample settings from the categorical distribution $\operatorname{cat}(\boldsymbol{\pi})$. Sampling from $\mathcal{D}_{r}$ below means drawing a history–target window from its trajectories. With denominator stabilizer $\epsilon=10^{-8}$, the expected next-state loss $\ell_{\mathrm{train}}(\mathbf{w})$ supervises the complete predictor–program pair:

The norms run over valid spatial cells of each field vector, excluding padded cells. Scored channels are weighted equally within each setting. Gradients through executed numerical arithmetic train the selected expert and shared encoder without intermediate labels. During rollout, each prediction enters the history window and replaces its oldest frame.

New-family adaptation retains this objective, reusing the shared encoder with a newly initialized system-specific predictor and a configured fixed numerical update program, as detailed in Appendix A.7.

Experiments

We follow MoE-POT’s benchmark Wang et al. [27] for joint pre-training and task-specific fine-tuning, complemented by ablations, within-family generalization, and transfer to a new PDE family.

Datasets and evaluation. Six settings cover five PDE families: NS at viscosities $\nu= 10^{-5}$ and $\nu= 10^{-3}$ ([11]), PDEBench CNS, SWE, and DR ([24]), and CFDBench velocity prediction ([13]). We retain the benchmark’s ten-frame inputs, padding, and masking, and evaluate full-horizon autoregressive predictions using mean relative $L_2$ error (L2RE). Joint pre-training and task-specific fine-tuning run for 1000 and 200 epochs, respectively, with final checkpoints used for evaluation.

Baselines and model scales. We compare eight baseline methods across fifteen configurations: six separately trained models FNO, UNet, FFNO, GK-T, Oformer, and GNOT and three scales each of pre-trained FNO, DPOT, and MoE-POT. We also report fine-tuned DPOT and MoE-POT results. Baseline results and activated-parameter counts are taken from MoE-POT Table 1. HotPOT-T/S use approximately 17M/90M activated parameters, matching MoE-POT-T/S. Appendices § 9.§ 4 and § 9.§ 5 provide detailed evaluation, training, and resource settings.

Model	Activated Params	NS benchmarks		PDEBench			CFDBench
		NS ($10^{-5}$)	NS ($10^{-3}$)	CNS(0.1, 0.01)	SWE	DR
Small models
FNO	0.5M	0.156	0.0128	0.170	0.00440	0.120	0.00761
UNet	25M	0.198	0.0245	0.357	0.0521	0.0971	0.0209
FFNO	1.2M	0.161	0.0256	0.183	0.00458	0.161	0.0990
GK-T	1.1M	0.260	0.0148	0.919	0.0453	0.0120	0.419
Oformer	1.8M	0.289	0.00319	0.161	0.00474	0.991	0.00444
GNOT	2.2M	0.590	0.316	0.533	0.00199	0.930	0.0216
Pre-trained
FNO-T	10M	0.191	0.0245	0.0859	11.0	0.530	0.00601
FNO-S	30M	0.157	0.0225	0.357	0.184	0.385	0.00460
FNO-M	150M	0.141	0.00730	–	0.0104	0.0112	–
DPOT-T	7.5M	0.107	0.0155	0.0168	0.00631	0.0577	0.00673
DPOT-S	30.8M	0.0688	0.00781	0.0244	0.00392	0.0367	0.00870
DPOT-M	122M	0.0569	0.00708	0.0224	0.00247	0.0288	0.0113
MoE-POT-T	17M	0.0682	0.00768	0.0105	0.00640	0.0411	0.00529
MoE-POT-S	90M	0.0552	0.00583	0.00959	0.00289	0.0342	0.00448
MoE-POT-M	288M	0.0528	0.00570	0.00914	0.00299	0.0300	0.00513
HotPOT-T	17M	0.0315	0.00325	0.00671	0.00177	0.0150	0.00219
HotPOT-S	90M	0.0299	0.00323	0.00623	0.00166	0.0326	0.00225
Fine-tuned
DPOT-T	7.5M	0.0700	0.00725	0.0168	0.00313	0.0289	0.00391
DPOT-S	30.8M	0.0502	0.00635	0.0238	0.00315	0.0215	0.00586
DPOT-M	122M	0.0424	0.00593	0.0221	0.00260	0.0175	0.00653
MoE-POT-T	17M	0.0456	0.00493	0.00746	0.00305	0.0188	0.00330
MoE-POT-S	90M	0.0361	0.00376	0.00742	0.00251	0.0182	0.00313
MoE-POT-M	288M	0.0351	0.00388	0.00744	0.00193	0.0140	0.00398
HotPOT-T	17M	0.0215	0.000481	0.00583	0.00162	0.00163	0.00208
HotPOT-S	90M	0.0193	0.000483	0.00536	0.00154	0.00436	0.00202

Main Experiments

Pre-training and fine-tuning. Across six settings, two activated-parameter scales, and both training stages, HotPOT outperforms MoE-POT in all 24 matched-scale comparisons (Table 1). Mean relative error reductions are 47.1% after pre-training and 54.3% after fine-tuning, or 50.7% across both stages. For example, after fine-tuning, HotPOT-T reduces L2RE relative to MoE-POT-T by 90.3% on NS at $\nu= 10^{-3}$ and 91.3% on DR. Appendices § 9.9–§ 9.10 provide paired field visualizations and space–time error plots.

Accuracy with fewer activated parameters. The 17M HotPOT-T model also outperforms the 90M and 288M MoE-POT models on every setting in both stages. Across the twelve setting–stage comparisons, its mean relative error reductions are 47.7% and 46.5%, respectively, while using approximately 81.1% and 94.1% fewer activated parameters.

Parameter accounting. Each HotPOT-T active path uses 16.77–17.18M parameters, comprising the 2.012M shared encoder and the selected system-specific predictor. Detailed parameter accounting and numerical query counts are provided in Appendix § 9.§ 5.

Ablation Studies

Full denotes HotPOT-T. All models use 17M configurations and fixed E1000 checkpoints; Mean equally weights six full-horizon L2REs.

Ablations of cross-PDE sharing. We test whether to share history encoding and where sharing belongs relative to private coupling. Learned private history gives each PDE family a separate trainable encoder, retaining expert predictors and numerical updates, to test history sharing. Late-shared evolution moves shared processing after private field coupling to test this sharing boundary.

In Table 2, Full lowers mean error by 3.2% relative to private history while reducing stored parameters by 7.6%. Private history performs better on SWE and slightly better on DR. Late sharing raises mean error from 0.0101 to 0.0286, worsening five settings. The placement controls in Appendix § 9.§ 6 favor history processing before private coupling over late evolution processing; sharing at the history stage further improves mean accuracy and reduces storage.

Model	NS $10^{-5}$	NS $10^{-3}$	CNS	SWE	DR	CFDB	Mean
Full HotPOT	0.0315	0.00325	0.00671	0.00177	0.0150	0.00219	0.0101
Learned private history	0.0322	0.00420	0.00790	0.000962	0.0149	0.00222	0.0104
Late-shared evolution	0.0331	0.00413	0.00941	0.00421	0.119	0.00217	0.0286

Ablations of numerical experts. We test the predictor–update interface across PDE families and, specifically for NS, whether its numerical update provides a benefit beyond repeated neural refinement. Direct-state predictor (Direct) replaces predicted quantities and numerical execution with residual state prediction, retaining the neural architecture, inputs, masks, and nominal readout-call counts. This changes the interface and stage-state queries together. Generic iteration targets the NS expert, replacing its numerical update with repeated, state-dependent residual predictions using shared readout weights. Other expert paths retain their updates.

In Table 3, Full lowers mean error from 0.0324 to 0.0101 relative to direct prediction, improving five settings; Direct is slightly better on CFDB. Relative to Generic iteration, Full lowers NS error by 1.5% at $\nu= 10^{-5}$ and 19.7% at $\nu= 10^{-3}$. These comparisons support the paired predictor–update interface for aggregate accuracy and NS numerical execution over generic learned iteration on both viscosity settings. Appendix § 9.§ 6 details the controls and temporal analyses; Appendix § 9.8 reports physical diagnostics.

Model	NS $10^{-5}$	NS $10^{-3}$	CNS	SWE	DR	CFDB	Mean
Full HotPOT	0.0315	0.00325	0.00671	0.00177	0.0150	0.00219	0.0101
Direct-state predictor	0.0525	0.00625	0.0104	0.01900	0.122	0.00211	0.0324
Generic iteration	0.0320	0.00404	0.00720	0.00177	0.0150	0.00234	0.0104

Within-Family Generalization

We evaluate pre-trained 17M models at the original resolutions under unseen conditions within supported PDE families.

Zero-shot prediction. For CNS, $M$ is the Mach number and $\eta$ the shear viscosity. The bulk viscosity is $\zeta= \eta$ in the reported datasets Takamoto et al. [24]. Fixed pre-training checkpoints predict CNS at $(M, \eta) = (1, 10^{-2})$ and CFDB boundary-condition, geometry, and physical-property subsets within the cavity and cylinder families, without target training or checkpoint selection. We report L2RE over ten forecast steps, denoted H10, with case-uniform masked aggregation for CFDB. This differs from the full-horizon main benchmark. Table 4 shows that HotPOT lowers error relative to MoE-POT in all seven groups, including a 55.5% CNS reduction.

Model	$(M, \eta) = (1, 10^{-2})$	BC	Geo†	Prop	BC	Geo†	Prop
HotPOT-T	0.189	0.00568	0.0169	0.00467	0.00119	0.00180	0.000577
MoE-POT-T	0.423	0.0276	0.0439	0.0116	0.00802	0.0117	0.00223

BC, Geo, and Prop denote boundary-condition, geometry, and physical-property subsets. ^†Each geometry subset contains one held-out case within its respective cavity or cylinder family.

Few-shot adaptation. Pre-trained models adapt to NS ($\nu= 10^{-4}$) and CNS ($(M, \eta) = (1, 10^{-2})$) using 16, 64, or 256 training trajectories. On 200 test trajectories per target, HotPOT achieves lower full-rollout L2RE in all six comparisons, with error decreasing across budgets (Table 5). Rollout validation selects checkpoints; Appendix § 9.§ 4.1 details the protocol.

Model / shots	16	64	256	16	64	256
HotPOT-T	0.124	0.0964	0.0819	0.0983	0.0599	0.0365
MoE-POT-T	0.304	0.206	0.108	0.365	0.153	0.0735

Transfer to a New PDE Family

Wave-Layer and Wave-Gauss Herde et al. [7] belong to one wave-equation family absent from both pre-training mixtures. HotPOT jointly trains its transferred encoder and a new scalar predictor. Scratch uses the same architecture, initial predictor weights, residual updates, and training protocol, but a random encoder. This comparison tests whether shared encoder pre-training benefits adaptation to a new PDE family.

At 16, 32, 64, and 128 shots, rollout validation selects checkpoints evaluated on 200 held-out trajectories per dataset using mean per-step L2RE. Encoder transfer outperforms both baselines in all eight comparisons. Error reductions are 11.9–39.8% relative to MoE-POT full fine-tuning; relative to Scratch, they are 6.5–8.4% on Wave-Layer and 0.6–16.4% on Wave-Gauss, supporting transfer of the pre-trained history representation. Appendix § 9.7 gives results (Table 11), protocol, and splits.

Target	Model	16	32	64	128
Wave-Layer	MoE-POT full fine-tuning	0.576	0.517	0.435	0.365
Wave-Layer	HotPOT scratch	0.553	0.428	0.331	0.235
Wave-Layer	HotPOT encoder transfer	0.506	0.395	0.304	0.220
Wave-Gauss	MoE-POT full fine-tuning	0.514	0.452	0.394	0.325
Wave-Gauss	HotPOT scratch	0.542	0.372	0.284	0.225
Wave-Gauss	HotPOT encoder transfer	0.453	0.370	0.279	0.218

Discussion and Conclusion

HotPOT combines shared field-history encoding with PDE-specific neural–numerical experts, allowing different equations to use different prediction quantities and fixed update operations. Across six pre-training settings on regular grids, HotPOT outperforms MoE-POT at matched activated-parameter scales before and after fine-tuning. Ablations support the sharing and prediction–execution designs, while transfer to an unseen wave-equation family supports encoder reuse with limited target data. System-specific experts increase storage as the model covers more PDEs, and adapting to a new family requires training a predictor and configuring its numerical update. Multi-stage updates can also increase inference cost through repeated neural readouts. Future work will explore reusable experts, automatic numerical-program composition, and efficient execution for broader PDE collections.

AI Use Statement

Generative AI tools assisted with language editing.

Reproducibility Statement

The appendix specifies the data interface, model dimensions, system-specific numerical updates, optimization schedule, parameter accounting, and checkpoint selection. Code, configurations, and fixed evaluation scripts are available in the public repository linked in the abstract.

Implementation and Additional Results

Backbone implementation

The patch encoding, temporal aggregation, and grouped Fourier/local blocks below build on DPOT and MoE-POT ([5, 27]).

All inputs use ten saved frames at a common $128^2$ analysis resolution. The deterministic feature map constructs six channels per scalar field: normalized value, two coordinates, boundary distance, geometry support, and field validity. Both scales use patch size 8, four adaptive Fourier neural operator (AFNO) channel groups, and 16 retained modes. The scale-dependent dimensions are listed in Table 6. The context decoder expands to full resolution before two $3 \times3$ convolutions and a $1 \times1$ typed readout. Per-field history means and standard deviations are computed over the ten frames and spatial grid, with standard deviations floored at $10^{-6}$. The system-specific mixers receive these means and log-scales in addition to shared features and task embeddings; the statistics have no learned parameters. Coordinates span $[-1,1]$; boundary distance is constant on periodic grids and the minimum coordinate distance to an edge otherwise. These are input-feature conventions, separate from the expert’s numerical update program.

Component	HotPOT-T (17M active)	HotPOT-S (90M active)
Shared width / depth / FFN width	256 / 2 / 512	512 / 3 / 2048
System-specific width / depth / FFN width	512 / 5 / 2048	1024 / 7 / 3968
Decoded context channels	32	64
Typed-head hidden width	128	256
Task-feature channels / hidden width	16 / 64	32 / 128
Task-embedding hidden width	128	128

Input encoding and positional embedding. For a grid of $N_x \times N_y$ cells, flattening only the spatial axes gives the feature matrix $\mathbf{X}^j \in\mathbb{R}^{N_xN_y \times d_{\mathrm{in}}}$ at history position $j \in\{1,\ldots,T\}$, with $d_{\mathrm{in}}=6$. The implementation embeds these six features in an 11-slot tensor with five fixed zero slots before the patch convolution. This fixed embedding is included in $f_{\mathrm{patch}}$; parameter counts include the stored convolution weights. A convolutional patch encoder $f_{\mathrm{patch}}$ maps non-overlapping $p \times p$ patches to $d$-dimensional tokens, where $p=8$ and $d=d_{\mathrm{enc}}$ in this encoder. With $P=(N_x/p)(N_y/p)$ and learned positional matrix $\mathbf{E}_{\mathrm{pos}} \in\mathbb{R}^{P \times d}$,

The patch encoder consists of a strided convolution, GELU (Hendrycks & Gimpel, [6]), and a pointwise projection. Its parameters, including $\mathbf{E}_{\mathrm{pos}}$, are part of $\theta$ and are shared across history frames, scalar fields, and datasets. Spatial convolutions and FFTs reshape these matrices to the two-dimensional grid and flatten them back in the same order.

Temporal aggregation. With field, setting, and saved-time indices suppressed, let $\mathbf{z}_i^j=(\mathbf{Z}_{i,:}^j)^\top\in\mathbb{R}^d$ be the token vector at patch $i$. The matrix $\mathbf{Z}^{\mathrm{hist}} \in\mathbb{R}^{T \times d}$ in Figure 2 stacks $\{(\mathbf{z}_i^j)^\top\}_{j=1}^{T}$. A real-valued adaptation of Fourier-feature aggregation uses learned channel-mixing matrices $\mathbf{W}_j \in\mathbb{R}^{d \times d}$ and frequency vector $\boldsymbol{\kappa} \in\mathbb{R}_{>0}^{d}$. Let $\lambda_j=(j-T)/(T-1)$ denote the normalized temporal lag, $\odot$ elementwise multiplication, and $\mathbf{v}_i^0 \in\mathbb{R}^d$ the aggregated token. The cosine acts componentwise:

Let $\mathbf{V}^0 \in\mathbb{R}^{P \times d}$ denote the matrix formed by stacking $\{(\mathbf{v}_i^0)^\top\}_{i=1}^{P}$.

Fourier and local mixing blocks. Spatial processing uses the grouped spectral mixer inherited from the MoE-POT implementation, following the AFNO/Fourier operator-transformer lineage ([4]. The mixer applies a two-dimensional fast Fourier transform (FFT), splits feature channels into groups, transforms each group with a two-layer frequency-domain multilayer perceptron (MLP), and returns to the spatial grid by inverse FFT. The MLP parameters are shared over retained frequency locations; GELU acts separately on real and imaginary components, with no soft-thresholding. Let $f_{\mathrm{fft}}$ and $f_{\mathrm{ifft}}$ denote the spatial FFT and its inverse, and let $f_{\mathrm{freq},b}$ be the grouped MLP with unretained modes set to zero. For block index $b$ and input $\mathbf{U} \in\mathbb{R}^{P \times d}$, including the inherited internal skip connection gives

Each block combines this spectral path with a local convolutional function $f_{\mathrm{local},b}$ and a channel feed-forward function $f_{\mathrm{ffn},b}$. The functions $f_{\mathrm{gn},b,1}$ and $f_{\mathrm{gn},b,2}$ are the block’s two group normalizations ([28]). Let $\mathbf{V}^{b}$, $\mathbf{U}^{b}$, $\widetilde{\mathbf{V}}^{b}$, and $\mathbf{V}^{b+1} \in\mathbb{R}^{P \times d}$ denote the block input, normalized input, mixed features, and output, respectively:

All spatial feature matrices have shape $P \times d$, where $d$ is the width of the shared or system-specific block in question. The local path uses a depthwise $3 \times3$ convolution, GELU, and a $1 \times1$ convolution; the FFN uses two $1 \times1$ convolutions with GELU. After the shared blocks and a final normalization, the encoder produces $Z_{r,c}$; their parameters belong to $\theta$. The corresponding system-specific blocks use $\phi_{\tau(r)}$.

System-specific numerical update programs

Each expert instantiates Eq. (5) by pairing a system-specific neural predictor with a numerical update program. Field-wise history encoding, system-specific state evolution, and spatial context decoding run once per saved interval. The program then constructs query states and reuses the predictor’s readout. All predicted quantities are learned through next-state supervision; they require no separate rate, action, or increment labels. The programs below detail the five alternatives illustrated in Figure 2. We suppress the setting index $r$ below; $N$ denotes the grid-cell count and a superscript $+$ the updated state. Vectors collect spatial entries; multi-channel arrays are written as matrices when channels form columns. The all-ones vector is $\mathbf{1}_{N} \in\mathbb{R}^{N}$, and $\langle a\rangle= N^{-1}\mathbf{1}_{N}\mathbf{1}_{N}^{\top}a$ is the spatial mean of a field vector $a$.

NS: transport-rate integration

The neural transport-rate query in Eq. (4) uses stage vorticity, history context, and deterministic Biot–Savart transport and forcing features. Let $f_{\mathrm{rate}} : \mathbb{R}^{N} \to\mathbb{R}^{N}$ denote this query with the decoded history context and physical coefficients held fixed. Using the diffused and midpoint states $\omega_{d}$, $\omega_{m}$ from Eq. (6), its two outputs are

The program adds forcing to each rate and uses the explicit midpoint update between two diffusion half-steps, forming a Strang split (Strang, [20]). The same decoded context and readout weights are reused at both queries and throughout the saved interval.

For the diffusion operator $f_{\mathrm{diff},\nu}^{s}$, let $\delta t_{\mathrm{CN}} > 0$ be a Crank–Nicolson time increment, $s$ the diffusion duration, and $m = s/\delta t_{\mathrm{CN}} \in\mathbb{N}$ its number of increments. Let $f_{\mathrm{fft}}$ and $f_{\mathrm{ifft}}$ be the discrete Fourier transform on the periodic grid and its inverse. For Fourier mode $\ell$, let $\lambda_{\ell} \ge0$ be the eigenvalue of the negative Laplacian. With viscosity $\nu$, define the mode multiplier $b_{\ell}$ and the diffusion update of a vorticity vector $\omega\in\mathbb{R}^{N}$ by

This composes $m$ Crank–Nicolson diffusion increments ([2]). The implementation uses a positive multiplier, with $\delta t_{\mathrm{CN}}\nu\lambda_{\ell}/2 < 1$, and evaluates its power in logarithmic form. Each half-step takes $s = \delta t/2$; composing the split substeps covers the saved interval $\Delta t_{r}$.

CNS: conservative action and admissible update

At cell $i$, let $u_{x,i}, u_{y,i}$ denote velocity components, $\rho_i$ density, and $p_i$ pressure. Let $\gamma> 1$ be the specific-heat ratio; the total energy density is $E_i = p_i/(\gamma- 1) + \rho_i(u_{x,i}^{2} + u_{y,i}^{2})/2$. The conservative vector is $\mathbf{c}_i = (\rho_i, \rho_i u_{x,i}, \rho_i u_{y,i}, E_i)^{\mathsf{T}} \in\mathbb{R}^{4}$. The full conservative state is the matrix $\mathbf{Q} \in\mathbb{R}^{N \times4}$ with row $\mathbf{c}_i^{\mathsf{T}}$. One predictor query supplies a normalized action, reshaped as $\mathbf{A} \in\mathbb{R}^{N \times4}$. Its channel means are $\bar{\mathbf{a}} = N^{-1}\mathbf{A}^{\mathsf{T}}\mathbf{1}_{N} \in\mathbb{R}^{4}$. Let $\epsilon_s > 0$ be a numerical scale floor, and let $f_{\mathrm{scale}} : \mathbb{R}^{N \times4} \to\mathbb{R}^{4}$ return the per-channel spatial root mean square (RMS) of centered conservative fields, floored at $\epsilon_s$; $\operatorname{diag}$ places these scales on a matrix diagonal. Denote the conservative increment by $\Delta\mathbf{Q}$ and its candidate state by $\mathbf{Q}^{\star}$. The program constructs

Multiplication by the diagonal matrix scales each channel uniformly over the grid. The program restores the candidate conservative means to those of $\mathbf{Q}$. Let $\rho_{\min}, p_{\min} > 0$ be density and pressure floors, and let $\alpha_{\mathrm{CNS}} \in[0,1]$ be one global coefficient ensuring these strict bounds in the updated state $\mathbf{Q}^{+}$. The segment update is

For an admissible starting state, an admissible candidate uses $\alpha_{\mathrm{CNS}} = 1$; otherwise a prescribed number $B_{\mathrm{bis}}$ of bisections determines an admissible lower endpoint of the segment. Boolean tests select the segment, and autograd differentiates its executed arithmetic. The result is converted back to primitive variables. The preserved quantities are the global totals of the four conservative channels on the uniform grid. Proposition 2 establishes the admissibility and conservation properties.

DR: two-query proposal refinement

The sum of two predictor output groups gives a finite-interval concentration proposal $\delta_1 \in\mathbb{R}^{2N}$ at the current state $\mathbf{q}^{n}$; the two concentration fields are concatenated in the same order as $\mathbf{q}^{n}$. A second query at an intermediate state gives $\delta_2$; the program averages the two proposals as in Eq. (9). Both queries use the same context and readout weights. These outputs are concentration proposals learned through the resulting next-state error; the fixed operation is the query-and-average refinement.

SWE: mass-preserving height increment

Only height is observed, with current field vector $\mathbf{h}^{n} \in\mathbb{R}^{N}$ and components $h_i^{n}$. One predictor query directly predicts a height increment $\Delta\mathbf{h} \in\mathbb{R}^{N}$ in the observed field units. The program centers it and applies one coefficient $\alpha\in[0,1]$ to the entire grid, as specified in Eq. (8). Let $h_{\min} > 0$ be a prescribed height floor and $\epsilon_d > 0$ a numerical denominator floor. For an initial height satisfying $h_i^{n} \ge h_{\min}$, the coefficient in Eq. (20) ensures $h_i^{n+1} \ge h_{\min}$ at every cell $i$. Since the increment has zero spatial mean and $\alpha$ is uniform, the update preserves mass in exact arithmetic. The role of the limiter is to shorten the complete increment, rather than clip individual cells. Proposition 1 gives the formula and proof.

CFDB: next-velocity assignment

One predictor query estimates the next two velocity fields, concatenated as a vector in $\mathbb{R}^{2N}$. The program replaces the current velocity channels and copies auxiliary support channels unchanged, before the output mapping restores the evaluation resolution and channel layout. Within CFDB, spectral/local layers are shared across geometry variants, while pointwise FFNs and velocity heads are geometry-specific.

Properties of the numerical updates

We establish the update properties used by the SWE and CNS numerical programs on a uniform grid. All arrays have finitely many finite-valued entries, arithmetic is exact, and each admissibility factor is a scalar shared by all spatial cells of one sample.

denominator floor $\epsilon_d > 0$, denote the cellwise step bounds by $r_i$, their clipped minimum by $\alpha$, and the updated height by $\mathbf{h}^{+}$. Define

Then $0\le\alpha\le1$, $h_i^{+}\ge h_{\min}$ for every cell, and $\sum_i h_i^{+}=\sum_i h_i$.

Proof. Every ratio is nonnegative. For $d_i\ge0$, the update cannot decrease $h_i$. For $d_i<0$, the denominator in Eq. (20) is at least $-d_i$, so $\alpha(-d_i)\le h_i-h_{\min}$, which gives the height bound. Since the same $\alpha$ multiplies every cell, $\sum_i h_i^{+}=\sum_i h_i+\alpha\sum_i d_i=\sum_i h_i$. If $\mathbf{d}=0$, then $\alpha=1$ and the state is unchanged. A cell at the floor with a negative increment instead gives $\alpha=0$.

Centering the directly predicted height increment gives $\mathbf{d}$. The denominator floor can shorten the permitted step for very small negative increments, but leaves both guarantees intact. The common scale is essential to this construction: cellwise clipping would generally alter the zero-sum increment and hence its mass balance.

Optimization and evaluation

Joint pre-training. Both scales train for 1000 epochs, with 545 updates per epoch, using two devices and local batch 20 (global batch 40). Each trajectory contributes one ten-frame-to-one-frame training window per epoch. Over a ten-epoch cycle, sampling visits the initial window once and early, middle, and late non-initial windows three times each. Batches contain compatible tasks; the epoch plan retains the benchmark’s per-dataset trajectory exposure. Adam uses $\beta_1 = \beta_2 = 0.9$ and weight decay $10^{-6}$. OneCycleLR reaches $10^{-3}$ after 20% of training, with initial and final division factors both $10^4$. NS (on its native grid) and DR receive norm-scaled Gaussian input noise with factor $5 \times10^{-4}$; CNS, SWE, CFDB, and all targets remain clean. The loss is the one-step physical-channel relative error in Eq. (11). The 9000 CFDB training rows comprise 581 cavity, 8231 cylinder, and 188 tube examples; all three geometry paths therefore participate in pre-training.

Fine-tuning on the pre-training tasks. Starting from the fixed E1000 checkpoint, the shared encoder and selected system-specific predictor train for 200 complete epochs; other predictors remain frozen. Both scales use global batch 40, the same Adam settings and 20% warmup, and target-only phase-aware sampling. Maximum learning rates are specified per task in Table 7; input noise is disabled. Fine-tuning retains full per-epoch target exposure, giving 5,000 updates for each NS setting, 45,000 for CNS and CFDB, and 4,600 for SWE and DR. Learning rates are fixed from the T-scale held-out E200 rollout-validation diagnostics and reused at both scales. Final refits use the complete original training populations and the fixed E200 checkpoint.

Table 7: Main evaluation populations and fine-tuning schedule for the pre-training tasks at both model scales. Train is the trajectory count used per epoch; test is the fixed evaluated row count. $C_{\mathrm{pred},r} = \lvert O_r \rvert$ counts scored physical channels, and $L_{\mathrm{roll}}$ counts model-feedback prediction steps after ten input frames. The learning rate is the maximum during E200 fine-tuning.

Setting	Native grid	Train	Test	$C_{\mathrm{pred},r}$	$L_{\mathrm{roll}}$	Fine-tune LR
NS $10^{-5}$	$64^2$	1,000	200	1	10	$10^{-5}$
NS $10^{-3}$	$64^2$	1,000	200	1	20	$10^{-3}$
CNS	$128^2$	9,000	200	4	11	$10^{-4}$
SWE	$128^2$	900	60	1	91	$10^{-3}$
DR	$128^2$	900	60	2	91	$10^{-3}$
CFDB	$64^2$	9,000	1,000	2	10	$10^{-5}$

Rollout metric and population. After the first ten frames, every predicted frame is fed back to the model. For each trajectory and physical channel, the released relative-$L_2$ metric flattens the spatial grid and the complete reported prediction window together. Channel errors are then averaged over scored channels, and trajectory errors equally over evaluated rows; the denominator includes $10^{-8}$. Masks exclude padded cells and unscored channels. Thus an H10 entry measures the whole ten-frame prediction window. CFDB scores velocity channels, with geometry/support channels used only as inputs and state assembly. The fixed main population uses the first 1000 rows (zero-based indices 0–999) of the 1182-row packed CFDB test file, comprising 110 cavity and 890 cylinder rows. It is row-weighted; the separate condition study uses case-uniform aggregation. Main and ablation results use fixed final checkpoints. Reported runs use one training seed; few-shot and wave-transfer runs use seed 0.

Few-shot selection protocol. NS at $\nu= 10^{-4}$ supplies 9800 target-training trajectories, and CNS at $M = 1$ supplies 9000. The last 256 trajectories of each training source form a fixed validation set. The reported training subsets contain 16, 64, or 256 trajectories sampled from the remaining rows. Paired methods use identical training and validation indices. Training runs for 500 epochs with batch min(20, shots), Adam with $\beta_1 = \beta_2 = 0.9$, weight decay $10^{-6}$, and OneCycleLR with 20% warmup and division factors $10^4$. NS HotPOT uses maximum learning rates $10^{-4}$, $10^{-4}$ and $10^{-5}$ at the three budgets, respectively, selected from the evaluated $\{10^{-4}, 10^{-5}\}$ candidates by rollout validation. MoE-POT uses $10^{-4}$ on both targets; CNS HotPOT also uses $10^{-4}$. Each training trajectory supplies one phase-aware ten-frame-to-next-frame window per epoch. Every ten epochs, validation evaluates autoregressive rollout on all 256 validation trajectories over the target test horizon. The checkpoint with lowest finite validation L2RE is used for testing. Both models update their shared and selected task parameters during adaptation. MoE-POT starts from its local E1000 pre-training checkpoint.

Few-shot rollout evaluation

Each target–budget pair uses its validation-selected checkpoint, including normalization statistics. Evaluation starts from ten input frames and uses model feedback for 20 steps on NS and 11 on CNS, over 200 test trajectories per target. Masked per-channel space–time L2RE is averaged equally over trajectories.

Parameter storage and numerical query counts

Table 8 gives exact parameter counts at both scales; Table 9 compares activated and stored parameters across models.

Quantity	HotPOT-T	HotPOT-S
All-PDE shared history encoder	2,011,648	11,901,952
Active NS/SWE path	16,770,577	89,584,673
Active DR path	16,903,188	90,112,036
Active CFDB geometry path	16,905,234	90,116,130
Active CNS path	17,178,004	91,185,956
Stored total	97,902,720	518,480,368

Model	Activated	Stored total
MoE-POT-T	17	30
HotPOT-T	16.77–17.18	97.90
MoE-POT-S	90	166
HotPOT-S	89.58–91.19	518.48
MoE-POT-M	288	489

Repeated use of the same predictor counts its weights once in the activated parameter total and contributes to runtime according to the number of calls.

Each setting updates the shared encoder and its selected family-specific predictor. The two NS viscosity settings use the same NS predictor. For CFDB, the encoder and geometry-shared layers train together with the selected geometry’s pointwise FFNs and velocity head. The family selector is deterministic and has no learned router.

Numerical query counts. The history encoder, system-specific context blocks, and spatial decoder run once per saved interval. NS uses two readout calls per numerical substep; DR uses two per saved interval, and CNS, SWE, and each CFDB geometry path use one. Repeated calls reuse the same readout weights, so they increase computation without increasing the number of activated parameters. The present configuration has no matched wall-clock benchmark; parameter counts and query counts describe different aspects of its computational cost.

Additional ablation details and analyses

History sharing and placement. Full uses one learned history encoder across PDE families, followed by private coupling and evolution. Learned private history copies the same encoder architecture into a separately trainable encoder for each family. The two NS settings retain the same family assignment. Predictor architectures, numerical updates, task exposure, and active-path dimensions are retained; storage increases because the history weights are no longer shared. Late-shared evolution uses field-history tokenization before private coupling and moves the shared processing to the evolution stage. Late-private evolution uses the corresponding placement with family-private weights. Table 10 summarizes the four combinations.

Model	Placement	Weights	Stored (M)	Mean L2RE
Full	History	Shared	97.90	0.0101
Learned private history	History	Private	105.95	0.0104
Late-shared evolution	Evolution	Shared	97.90	0.0286
Late-private evolution	Evolution	Private	105.95	0.0287

Direct-state and iterative controls. Direct abbreviates the Direct-state predictor control. Full and Direct retain the same encoder, system-specific mixers, context decoders, and readout-head architectures, with 97,902,720 stored parameters and identical activated counts. Their history features, statistics, task inputs, and masks match. Direct replaces the prediction–update interface with residual state prediction, retaining nominal readout-call counts. It therefore changes prediction quantities, numerical execution, and stage-state queries together. Generic iteration uses the same cached history context and shared readout weights to refine the NS state through state-dependent residual queries. Its intervention tests the NS numerical update against a learned iterative alternative. Other paths retain their configured update structures. All rows are trained under the same E1000 multi-PDE protocol.

Space–time visualization of two-dimensional tasks. The SWE and DR tasks have two spatial dimensions, $x$ and $y$. The two DR concentration fields are denoted by $u$ and $v$. Figure 3 stacks their two-dimensional error fields over forecast steps to form $(x, y, n)$ space–time volumes, where $n$ indexes the forecast step. The third rendering axis therefore represents time, not a third spatial dimension. We compare HotPOT with Late-shared evolution using fixed example 0 and all 91 forecast steps. Two complementary two-dimensional cuts accompany each volume pair: an $x$–$y$ map at step 91 locates the final spatial error differences, and an $x$–time map at zero-based grid index $y = 64$ shows how the local differences evolve. These cuts expose local differences without occlusion by the volume rendering.

In these illustrated trajectories, the DR spatial cuts contain extended regions of lower HotPOT error, and the time cuts show differences persisting across many forecast steps. For SWE, the time cut localizes differences along the moving wave pattern. These local views complement the six-setting aggregate comparison in Table 2.

Full-population temporal errors. Figure 4 aggregates every trajectory in the six main-test populations. It connects the final aggregate scores in Table 2 and Table 3 to the evolution of prediction error through time. The per-step metric differs from full-horizon L2RE, which forms one norm over the complete space–time window.

Figure 5 shows where Full and Generic iteration differ on four matched NS trajectories.

Few-shot transfer to the wave family

Data and initialization. Wave-Layer and Wave-Gauss are two datasets from one wave-equation family in the Poseidon benchmark [7]. Neither belongs to either source checkpoint’s six-dataset pre-training mixture. For each target, we use a fixed custom split of mutually disjoint original trajectories: 1,800 in the training pool, 100 for validation, and 200 for testing. Each original trajectory belongs to exactly one split; sampled temporal windows remain within that split. All methods use these same splits. The two targets are evaluated separately at 128 × 128 resolution with ten input frames and 11/6 forecast frames for Layer/Gauss. Nested subsets contain 16, 32, 64, or 128 training trajectories.

HotPOT transfer initializes the 2.012M-parameter field-wise history encoder from the 17M E1000 source checkpoint. A scalar system-specific predictor is randomly initialized, and all 16.771M target-model parameters are trained. The scratch control has the same architecture and identical initial predictor weights, differing only in encoder initialization. Both use a common residual-prediction architecture that adds the predicted increment to the current displacement. Each model retains its native input normalization. MoE-POT loads its full local 1000-epoch pre-training checkpoint, updates all its trainable parameters, including the feature extractor and gating parameters.

Training and testing. All methods use the same target splits, sampled windows, training budget, optimizer, learning-rate schedule, and validation-selection rule. All runs use 500 epochs, batch $\min(20,\mathrm{shots})$, Adam with $\beta_{1} = \beta_{2} = 0.9$ and weight decay $10^{-6}$, and OneCycleLR with maximum $10^{-4}$, 20% warmup, and both division factors $10^{4}$. Each epoch samples one ten-frame-to-nextframe window per trajectory; the sampling indices and order are matched across arms. Validation evaluates the full rollout every ten epochs. The lowest validation error selects the checkpoint, with the earliest checkpoint breaking ties. The selected model is evaluated on the 200 held-out test trajectories.

The metric is the trajectory-uniform mean of per-step spatial L2RE, with denominator floor $10^{-8}$ and each model’s own autoregressive feedback.

Results. Encoder transfer lowers error relative to MoE-POT full fine-tuning by 12.1–39.8% on Wave-Layer and 11.9–32.9% on Wave-Gauss. It also improves on scratch training in all eight target–budget pairs, with reductions of 6.5–8.4% on Layer and 0.6–16.4% on Gauss. The comparison with scratch measures the effect of pre-trained encoder initialization under this matched target architecture and training protocol.

Physical diagnostics

For NS, relative enstrophy error compares the spatial mean squared vorticity with the reference at each forecast step and divides by the reference value, floored at $10^{-8}$. For SWE and CNS, conservation drift measures the absolute change of a predicted spatial mean from the last observed state. SWE normalizes by the initial absolute mean height; CNS normalizes each conservative component by its initial spatial mean absolute value. Metrics average all forecast steps and then all test trajectories.

Full has lower mean enstrophy error on both NS settings and lower conservation drift on SWE and CNS than Direct in Table 12. Propositions 1 and 2 establish constraint preservation in exact arithmetic; the table measures the implemented floating-point rollouts. The CFDB program copies auxiliary support channels exactly.

System	Diagnostic	Full	Direct
NS $10^{-5}$	Relative enstrophy error	$1.78 \times10^{-3}$	$3.63 \times10^{-3}$
NS $10^{-3}$	Relative enstrophy error	$2.16 \times10^{-3}$	$2.51 \times10^{-3}$
CNS	Normalized mass drift	$7.63 \times10^{-8}$	$2.04 \times10^{-4}$
	Normalized energy drift	$8.15 \times10^{-8}$	$2.06 \times10^{-4}$
SWE	Relative mass drift	$1.19 \times10^{-9}$	$9.31 \times10^{-5}$

Observed numerical execution. To distinguish mean correction from admissibility limiting, we observe the actual E1000 Full program on all 60 SWE trajectories over 91 steps and all 200 CNS trajectories over 11 steps. At each step, we record the uncentered neural proposal, the mean-corrected candidate, and the executed state along the model’s own feedback trajectory.

Figure 6 shows the drift before and after mean correction. The limiter coefficient equals one at all 7,660 evaluated updates, so each mean-corrected candidate is executed unchanged. The minimum executed SWE height is 0.2519; the minimum CNS density and pressure are 0.05715 and 0.3181, respectively.

SWE reference compatibility. Across all 60 reference trajectories and their 91 forecast frames, float64 spatial reductions give the following compatibility measures. Relative to the last observed frame, the maximum absolute change in mean height is $1.20 \times10^{-9}$. After shifting each forecast reference to that initial mean, both the initial and mean-shifted reference fields have minimum height above 0.2559, well above the $10^{-6}$ floor. The shifted references are compatible with the mass-preserving update in Appendix § 9.§ 3. Let $\mathbf{h}^{\mathrm{ref},n} \in\mathbb{R}^{N}$ be the forecast reference field, $\bar{h}^{\mathrm{ref},n} = \langle\mathbf{h}^{\mathrm{ref},n} \rangle$ its spatial mean, and $\bar{h}^{\mathrm{obs}}$ the mean of the last observed input frame. The relative $L_2$ discrepancy associated with the initial mean, $\sqrt{N}|\bar{h}^{\mathrm{ref},n} - \bar{h}^{\mathrm{obs}}|/\|\mathbf{h}^{\mathrm{ref},n}\|_2$, averages $2.27 \times10^{-10}$ over trajectories and forecast steps.

Paired field visualizations

Figure 7 and Figure 8 compare reference fields, model predictions, and absolute errors across six settings, including separate cavity and cylinder examples. Labels $u_x$, $u_y$ denote velocity components in CNS and CFDB. All scored fields use shared color scales within each row. Both 17M-active models use fixed E1000 checkpoints, the same ten-frame history, and autoregressive feedback.

Temporal error profiles. Figure 9 follows the same seven examples through every forecast step, complementing the final-frame maps with their temporal error profiles.

Space–time error structure

We examine the SWE and DR trajectories from Appendix § 9.9, retaining all 91 forecast steps and both DR fields. Fixing zero-based grid index $y = 64$ gives $x$–time cuts with forecast time on the vertical axis (Figure 10). SWE shows outward-moving and returning waves; DR shows evolving concentration patterns. Color scales are shared within each field.