1 Introduction↩︎

Spatiotemporal dynamics are central to many problems in science and engineering, including shock formation, transport and mixing, turbulence, and large-scale environmental prediction. Many of these systems are naturally posed as initial-boundary value problems for partial differential equations (PDEs). Obtaining numerically reliable solutions at useful spatial and temporal resolutions, however, remains computationally demanding. Stability restrictions such as the Courant–Friedrichs–Lewy condition, the need to resolve disparate scales, and repeated simulations for design, control, or uncertainty quantification can together lead to large wall-clock costs even on modern hardware. Classical finite-difference, finite-volume, finite-element, and spectral methods offer well-established accuracy and stability properties, but they often require dense meshes and small time steps to control truncation error over long horizons.

Deep-learning-based surrogate models seek to amortize this cost by replacing repeated numerical solves with a training stage followed by fast inference. Data-driven models have achieved strong predictive accuracy on canonical PDE benchmarks, including convolutional neural network (CNN) models for Cartesian grids [1], graph neural network (GNN) simulators for irregular meshes [2], [3], and neural operators for learning mappings between function spaces across discretizations [4], [5]. Despite these successes, purely data-driven surrogates are often black-box models that rely on large amounts of labeled trajectory data and can generalize poorly outside the training distribution. This limitation is especially consequential for autoregressive rollout, where small one-step errors can compound into unstable or physically inconsistent long-term predictions.

Physics-informed deep learning (PiDL) addresses these limitations by incorporating knowledge of governing laws into the learning problem. Physics-informed neural networks (PINNs) [6] and physics-constrained surrogate models [7] reduce the dependence on labeled data by penalizing violations of governing equations, conservation laws, initial conditions, and boundary conditions. These approaches improve data efficiency and physical consistency, but several practical challenges remain. Multi-objective loss balancing can be sensitive to normalization and weighting, residuals from stiff or strongly nonlinear operators can be difficult to optimize, and satisfying residual constraints during training does not by itself guarantee stable long-term autoregressive prediction.

Recent advances in differentiable solvers have enabled another class of hybrid PiDL methods that couple neural networks with numerical operators and train the combined system end to end. For example, Bar-Sinai et al. [8] learned data-driven discretizations for PDEs, Kochkov et al. [9] accelerated fluid simulation by replacing parts of the Navier–Stokes solver with neural networks, and Fan et al. [10] used differentiable hybrid neural modeling for fluid–structure interaction. Other studies have shown that solver-in-the-loop training and embedded numerical operators can improve long-horizon prediction by introducing useful inductive bias [11]–[13]. Beyond finite-difference-style operators, finite-volume and finite-element ideas have also been incorporated into learnable simulation frameworks to improve conservation, geometric flexibility, and operator reuse [14]–[16].

Although these hybrid methods have substantially improved the accuracy and robustness of neural surrogates, important gaps remain. In many approaches, the trainable component is still a generic neural architecture whose internal representation is only weakly aligned with the conservative structure of the underlying dynamics. Moreover, applying embedded numerical operators on coarse meshes can reduce cost but may lose fine-scale information, whereas applying them densely at high resolution can inherit the computational burden and stability limitations of the original solver. These issues motivate architectures that encode conservation and multiscale structure directly into the prediction mechanism while retaining the flexibility of learnable models.

In this work, we propose the Multi-Resolution Finite-Volume-inspired network (MuRFiV), a deep architecture for long-term prediction of spatiotemporal dynamics. The central idea is to separate the evolution of subdomain averages from the reconstruction of high-resolution local structure. At the global scale, MuRFiV uses a finite-volume-inspired global model to predict changes in subdomain averages through fluxes across subdomain interfaces. At the local scale, a shared local model reconstructs intra-subdomain local fluctuations and enforces a zero-mean constraint so that the reconstructed fine-scale field remains consistent with the global average. The framework can be used as MuRFiV-Eq, where an embedded numerical operator provides physics-based flux information corrected by the global network, or as MuRFiV-NoEq, where the global network learns the interface fluxes directly.

1.0.0.1 Contributions.

• Finite-volume-inspired architecture for stable roll We design MuRFiV around domain decomposition, flux-based global updates, and local residual reconstruction, embedding conservative structure directly into the prediction architecture. A lightweight local network is shared across subdomains to reconstruct fine-scale local fluctuations, making high-resolution feature prediction parameter efficient. out.

• High-quality embedded operators through interface-only flux prediction. In MuRFiV-Eq, the embedded numerical operator is evaluated on the original fine mesh only along subdomain interfaces, where the global network corrects the fluxes used to update subdomain averages. This design preserves high-quality physics-based flux information at the locations that control global evolution, while keeping the cost low by avoiding dense operator evaluation over the full domain.

• Robust prediction across canonical PDE systems. Across one- and two-dimensional Burgers’ equations, shallow water equations, and incompressible Navier–Stokes equations, MuRFiV achieves accurate long-term autoregressive prediction, improves stability relative to data-driven baselines, and generalizes to unseen physical parameters.

Together, these contributions show that combining finite-volume-inspired conservation with multi-resolution learning provides a practical route to accurate, stable, and efficient neural surrogates for spatiotemporal dynamics.

2 Methodology↩︎

2.1 Problem Formulation↩︎

We consider modeling spatiotemporal dynamics governed by general conservative PDEs of the following form:

\[\begin{align} \frac{\partial \boldsymbol{u}}{\partial t} = \nabla \cdot \boldsymbol{\mathcal{F}}(\boldsymbol{u}, \nabla \boldsymbol{u}, \nabla^2\boldsymbol{u}, \cdots, \boldsymbol{\lambda}) + \boldsymbol{\mathcal{Q}}(\boldsymbol{u}),& & & \boldsymbol{x}, t \in \Omega \times [0, T], \boldsymbol{\lambda} \in \mathbb{R}^{m}\\ \boldsymbol{\mathcal{I}}\left[\boldsymbol{x}, \boldsymbol{u}, \nabla_{\boldsymbol{x}}^2\boldsymbol{u}, \nabla_{\boldsymbol{x}}\boldsymbol{u}\cdot\boldsymbol{u}, \boldsymbol{\lambda}\right] = \boldsymbol{0},& & & \boldsymbol{x} \in \Omega, t =0,\boldsymbol{\lambda} \in \mathbb{R}^{m}\\ \boldsymbol{\mathcal{B}}\left[t, \boldsymbol{x}, \boldsymbol{u}, \nabla_{\boldsymbol{x}}^2\boldsymbol{u}, \nabla_{\boldsymbol{x}}\boldsymbol{u}\cdot\boldsymbol{u}, \boldsymbol{\lambda}\right] = \boldsymbol{0},& & & \boldsymbol{x}, t \in \partial\Omega \times [0, T], \boldsymbol{\lambda} \in \mathbb{R}^{m}\label{eq:pde} \end{align}\tag{1}\]

where \(\boldsymbol{u}(\boldsymbol{x}, t)\) is the state variable while \(\boldsymbol{x}\) and \(t\) are the space and time coordinates, respectively. \(\boldsymbol{\mathcal{F}}\) is a non-linear operator parameterized by \(\boldsymbol{\lambda}\), while the \(\boldsymbol{\mathcal{Q}}\) represents any other terms that cannot be formulated as divergence (e.g. body force, volumetric source/sink term).

Traditionally, numerical methods are leveraged to resolve these spatiotemporal dynamics, where the numerical operator \(\boldsymbol{\mathcal{S}}\) can be formulated as: \[\begin{align} \boldsymbol{u}_{t_i} = \boldsymbol{\mathcal{S}}(\boldsymbol{u}_{t_{i - 1}},\boldsymbol{\lambda}),\quad i\in\{1,2,\cdots,N\}\\ \end{align}\] where \(\boldsymbol{u}_{t_i}\) denotes the state variable at the \(i^{th}\) time step (\(t_i\)). Please note that the time interval (\([t_{i-1}, t_i]\), \(\Delta T= t_i -t_{i-1}\)) is generally determined by the dynamics, which can be much longer than the actual numerical step size (\(\Delta t\), \(\Delta t \ll \Delta T\)) used by the solver due to numerical stability consideration, and the notation \(\boldsymbol{\mathcal{S}}\) may represent the composite operator of many numerical time iterations. In this work, we are interested in training a neural surrogate model \(\boldsymbol{\mathcal{N}}\) given the training dataset \(\mathbb{S}\) and a validation set \(\mathbb{V}\) such that: \[\begin{align} &\min_{\boldsymbol{\theta}} \mathcal{L}(\{\tilde{\boldsymbol{u}}_{t_i}\}, \mathbb{V})\\ \text{where, } & \begin{cases} \tilde{\boldsymbol{u}}_{t_i} = \boldsymbol{\mathcal{N}}(\tilde{\boldsymbol{u}}_{t_{i - 1}},\boldsymbol{\lambda};\boldsymbol{\theta})\\ \cdots\\ \tilde{\boldsymbol{u}}_{t_1} = \boldsymbol{\mathcal{N}}(\boldsymbol{u}_{t_{0}},\boldsymbol{\lambda};\boldsymbol{\theta}), \quad\boldsymbol{u}_{t_0} \in \mathbb{V} \end{cases} \end{align}\] where \(\boldsymbol{\theta}\) represents the trainable parameters of the neural surrogate model, and \(\mathcal{L}\) is the loss function. In this work, Mean Squared Error (MSE) was used as the loss function. Once the model is properly trained, the spatiotemporal trajectory can be predicted by rolling out the surrogate model auto-regressively.

2.2 Design Overview of MRFiV↩︎

MuRFiV consists of two major trainable components: the global network \(\boldsymbol{\mathcal{N}}_G\) and the local network \(\boldsymbol{\mathcal{N}}_L\), together with an optional component containing numerical operators. The detailed structure schematic of MuRFiV is shown in Fig. 1.

In particular, the global network is designed to capture large-scale dynamics while the local network focuses on learning small-scale features. Domain decomposition is leveraged to achieve multi-scale learning effectively and the entire computational domain is decomposed into connected and disjoint subdomains. In contrast to many previous hybrid neural solvers, where numerical operators work on down-sampled resolution for lower computational cost at the expense of degraded accuracy or missing details, the numerical operators in MuRFiV calculate the average change over each subdomain by evaluating the fluxes over subdomain interfaces based on the results at the original high resolution. In such a way, not only high accuracy is preserved but also it is achieved at a relatively low computational cost thanks to the spatial sparsity of the subdomain interfaces.

2.3 Domain Decomposition↩︎

For a \(D\)-dimensional computational domain defined in a Euclidean space \(\Omega = [l_d, h_d]^D\), where \(l_d\) and \(h_d\) represent the lower and higher boundary coordinates of the \(d^{th}\) dimension, respectively. We decompose the computational domain into \(N_s = \prod_{d=1}^D n_d\) connected and disjoint subdomains \({\Omega}_{\boldsymbol{j}}\), where \(\boldsymbol{j} = [j_1,\cdots,j_d,\cdots, j_D]\), \(\boldsymbol{j}\in \boldsymbol{J}\), \(j_d \in \{1,2,\cdots,n_d\}\) represents the index of the subdomain \(\Omega_{\boldsymbol{j}}\). In this work, we use uniform decomposition for ease of implementation, where every subdomain is an orthotope (e.g., line segment in 1D, rectangle in 2D, cuboid in 3D, ...) and is identical to the others.

2.4 Global Model↩︎

For conservative systems, the Finite Volume Method calculates the integral average of the state variable by applying Gauss’s Law to the governing equation (Eq. 1 ): \[\frac{\partial\bar{\boldsymbol{u}}_{\boldsymbol{j}}}{\partial t} V_{\boldsymbol{j}} = \iiint_{\Omega_{\boldsymbol{j}}} \frac{\partial \boldsymbol{u}}{\partial t}\text{ d}V = \oiint_{\partial \Omega_{\boldsymbol{j}}} \boldsymbol{\mathcal{F}}(\boldsymbol{u}, \nabla \boldsymbol{u}, \nabla^2\boldsymbol{u}, \cdots,\,\boldsymbol{\lambda})\cdot\boldsymbol{n} \text{ d}S + \iiint_{\Omega_{\boldsymbol{j}}}\boldsymbol{\mathcal{Q}}_i(\boldsymbol{u})\text{ d} V \label{eq:consPDE}\tag{2}\] where \(\bar{\boldsymbol{u}}_{\boldsymbol{j}}\) is the integration average over the subdomain \(\Omega_{\boldsymbol{j}}\) (\(\bar{\boldsymbol{u}}_{\boldsymbol{j}} = \frac{1}{V_{\boldsymbol{j}}}\iiint_{\Omega_{\boldsymbol{j}}} \boldsymbol{u} \text{ d}V\)). \(\boldsymbol{\mathcal{F}}\) is the numerical operator calculating flux at each subdomain interface and \(\boldsymbol{n}\) is the unit normal vector of a subdomain interface with an area of \(S\). For a given training set \(\mathbb{S}\), all state variables are stored at discretized time steps with a fixed time interval \(\Delta T\) (i.e. learning step size), where: \[\Delta \bar{{\boldsymbol{u}}}_{t_i, \boldsymbol{j}} = \bar{\boldsymbol{u}}_{t_i, \boldsymbol{j}} - \bar{\boldsymbol{u}}_{t_{i-1}, \boldsymbol{j}} = \int_{t_{i-1}}^{t_{i}} \frac{\partial\bar{\boldsymbol{u}}_{\boldsymbol{j}}}{\partial t} \text{ d}t +\epsilon\] where \(\epsilon\) is the discretization and truncation error of the numerical solver \(\boldsymbol{\mathcal{S}}\). Since this work aims to train a surrogate model, we assume the solution obtained from numerical solvers is the ground truth (\(\epsilon\rightarrow0\)). MuRFiV uses a global model \(\boldsymbol{\mathcal{M}}_G\) to predict the average change: \(\widetilde{\Delta \bar{{\boldsymbol{u}}}}_{t_i} = \boldsymbol{\mathcal{M}}_G(\boldsymbol{u}_{t_{i-1}}, \boldsymbol{\lambda};\,\boldsymbol{\theta}_G)\), where \(\boldsymbol{\theta}_G\) are the trainable parameters. In MuRFiV, a global neural network \(\boldsymbol{\mathcal{N}}_G\) together with optional embedded numerical operators \(\boldsymbol{\mathcal{S}}_e\) consists of the global model \(\boldsymbol{\mathcal{M}}_G\).

2.4.1 Modeling Flux between Subdomain↩︎

In MuRFiV, the optional embedded numerical operator (\(\boldsymbol{\mathcal{S}}_e\)) approximates the change of the integral average \(\widetilde{\Delta\bar{{\boldsymbol{u}}}}^s_{t_i, \boldsymbol{j}}\) using forward Euler scheme based on the previous learning step: \[\begin{alignat}{2} \widetilde{\Delta \bar{{\boldsymbol{u}}}}^s_{t_i, \boldsymbol{j}} &\approx \frac{\Delta T}{V_{\boldsymbol{j}}} \cdot \left. \frac{\partial\bar{\boldsymbol{u}}_{\boldsymbol{j}}}{\partial t}\right |_{t_{i-1}}\tag{3}\\ &= \frac{\Delta T}{V_{\boldsymbol{j}}}\cdot \left[\left.\oiint_{\partial \Omega_{\boldsymbol{j}}} \boldsymbol{\mathcal{F}}(\boldsymbol{u}, \nabla \boldsymbol{u}, \nabla^2\boldsymbol{u}, \cdots,\, \boldsymbol{\lambda})\cdot \tilde{\boldsymbol{n}} \text{ d}S\right|_{t_{i-1}} + \left.\iiint_{\Omega_{\boldsymbol{j}}}\boldsymbol{\mathcal{Q}}_i(\boldsymbol{u})\text{ d} V\right|_{t_{i-1}}\right]\tag{4}\\ &= \frac{\Delta T}{V_{\boldsymbol{j}}} \cdot \boldsymbol{\mathcal{S}}_e \left(\boldsymbol{u}\left(\boldsymbol{x}, t_{i-1}\right),\, \boldsymbol{\lambda}\right)\,,\quad \boldsymbol{x}\in \partial \Omega_{\boldsymbol{j}}\cup\partial \Omega^A_{\boldsymbol{j}} \end{alignat}\tag{5}\] where \(\Omega^A_{\boldsymbol{j}}\) denotes the adjacent subdomains of \(\Omega_{\boldsymbol{j}}\). The only error introduced in Eq. 5 is due to the large time step (\(\Delta T\)) used in the Euler discretization (Eq. 3 ) compared to the ground truth solver \(\boldsymbol{\mathcal{S}}\). While the spatial discretization used in Eq. 4 is mathematically equivalent to \(\boldsymbol{\mathcal{S}}\) attributed to applying Gauss’s Law on the fine mesh. Meanwhile, such operators \(\boldsymbol{\mathcal{S}}_e\) are only applied to the boundaries of each subdomain. By leveraging Gauss’s Law, MuRFiV’s numerical operator \(\boldsymbol{\mathcal{S}}_e\) can achieve high accuracy at a significantly lower computational cost.

2.4.2 Global Neural Network↩︎

When MuRFiV is equipped with the numerical operators \(\boldsymbol{\mathcal{S}}_e\) (MuRFiV-Eq), a global neural network \(\boldsymbol{\mathcal{N}}_G\) is designed to correct the error introduced by Eq. 3 due to large time step size \(\Delta T\): \[\boldsymbol{F}^{ \mathcal{N}}_c(\boldsymbol{x}, t_i) = \boldsymbol{\mathcal{N}}_G(\boldsymbol{u}_{t_{i-1}},\, \boldsymbol{\lambda};\,\boldsymbol{\theta}_G)\,,\quad\boldsymbol{x} \in \partial \Omega^b,\] where \(\boldsymbol{F}^{ \mathcal{N}}_c(\boldsymbol{x}, t_i)\) represents the correction of the integration of flux over time \([t_{i-1}, t_i]\) predicted by the global neural network, while \(\partial \Omega^b\) denotes the interfaces of all subdomains. Then the change of the integral average \(\Delta \bar{{\boldsymbol{u}}}_{t_i, \boldsymbol{j}}\) for each subdomain \(\boldsymbol{j}\) is approximated by \(\widetilde{\Delta \bar{{\boldsymbol{u}}}}_{t_i, \boldsymbol{j}}\) where: \[\widetilde{\Delta \bar{{\boldsymbol{u}}}}_{t_i, \boldsymbol{j}} = \boldsymbol{\mathcal{M}}_G(\boldsymbol{u}_{t_{i-1}}, \boldsymbol{\lambda};\,\boldsymbol{\theta}_G) = \frac{\Delta T}{V_{\boldsymbol{j}}} \boldsymbol{\mathcal{S}}_e \left(\boldsymbol{u}\left(\boldsymbol{x}, t_{i-1}\right),\,\boldsymbol{\lambda}\right) + \sum_{\boldsymbol{x}\in\partial \Omega_{\boldsymbol{j}}} \boldsymbol{F}^{ \mathcal{N}}_c(\boldsymbol{x}, t_i)\cdot \boldsymbol{n}\] When no equations are leveraged (MuRFiV-NoEq), the Global Network \(\boldsymbol{\mathcal{N}}_G\) will directly predict the fluxes (\(\boldsymbol{F}_{\boldsymbol{\mathcal{N}}}\)) instead of a correction of the fluxes. Thus the integral average change for each subdomain \(\boldsymbol{j}\) is calculated by: \[\widetilde{\Delta\bar{{\boldsymbol{u}}}}_{t_i, \boldsymbol{j}} = \boldsymbol{\mathcal{M}}_G(\boldsymbol{u}_{t_{i-1}}, \boldsymbol{\lambda};\,\boldsymbol{\theta}_G) = \sum_{\boldsymbol{x}\in\partial \Omega_{\boldsymbol{j}}} \boldsymbol{F}^{ \mathcal{N}}(\boldsymbol{x}, t_i)\cdot \boldsymbol{n}\]

2.5 Local Model↩︎

The local model \(\boldsymbol{\mathcal{M}}_L\) is designed to predict the variance of the state variable within each subdomain \(\boldsymbol{u}'_{t_i, \boldsymbol{j}}\) and the local model is shared among all the subdomains for efficient use of training data as well as trainable parameters. \(\boldsymbol{u}'_{t_i, \boldsymbol{j}}\) is defined as: \[\boldsymbol{u}'(\boldsymbol{x}, t_i) = \boldsymbol{u}(\boldsymbol{x}, t_i) - \bar{\boldsymbol{u}}_{t_i, \boldsymbol{j}} = \boldsymbol{u}(\boldsymbol{x}, t_{i-1}) + \Delta \bar{{\boldsymbol{u}}}_{t_i, \boldsymbol{j}}\,,\quad \boldsymbol{x}\in \Omega_{\boldsymbol{j}} \label{eq:g95l}\tag{6}\] Based on the definition above, the variable \(\boldsymbol{u}'_{t_i, \boldsymbol{j}}\) has a “zero-mean” property for each subdomain: \[\iiint_{\Omega_{\boldsymbol{j}}} \boldsymbol{u}'_{t_i, \boldsymbol{j}} \text{ d}V = \iiint_{\Omega_{\boldsymbol{j}}} \boldsymbol{u}(\boldsymbol{x}, t_i) \text{ d}V - \bar{\boldsymbol{u}}_{t_i, \boldsymbol{j}}\cdot V_{\boldsymbol{j}} = 0 \label{eq:zero95mean}\tag{7}\] The local model takes advantage of the “zero-mean” property as an extra normalization during the model forward pass when predicting the local variance \(\tilde{\boldsymbol{u}}'_{t_i, \boldsymbol{j}}\), which significantly improves the robustness of MuRFiV in the long-term rollout. \[\begin{alignat}{2} \tilde{\boldsymbol{u}}'(\boldsymbol{x}_{\text{out}}, t_i) &= \boldsymbol{\mathcal{M}}_L\left(\boldsymbol{u}(\boldsymbol{x}_{\text{in}}, t_{i-1}), \boldsymbol{\lambda};\boldsymbol{\theta}_L\right) \\ &= \tilde{\boldsymbol{u}}'^{ \mathcal{N}}(\boldsymbol{x}_{\text{out}}, t_i) - \frac{1}{V_{\boldsymbol{j}}} \iiint_{\Omega_{\boldsymbol{j}}} \tilde{\boldsymbol{u}}'^{ \mathcal{N}}(\boldsymbol{x}_{\text{out}}, t_i)\text{ d}V\label{eq:localnn95zm95c} \\ \text{where}\quad \tilde{\boldsymbol{u}}'^{ \mathcal{N}}(\boldsymbol{x}_{\text{out}}, t_i) &= \boldsymbol{\mathcal{N}}_L\left(\boldsymbol{u}(\boldsymbol{x}_{\text{in}}, t_{i-1}), \boldsymbol{\lambda};\boldsymbol{\theta}_L\right) + \boldsymbol{\mathcal{I}}n(\widetilde{\Delta\bar{{\boldsymbol{u}}}}_{t_i, \boldsymbol{j}}) \,,\; \boldsymbol{x}_{\text{in}}\in \Omega_{\boldsymbol{j}}\cup\Omega^A_{\boldsymbol{j}}\,,\; \boldsymbol{x}_{\text{out}}\in \Omega_{\boldsymbol{j}} \end{alignat}\tag{8}\] The \(\boldsymbol{\mathcal{N}}_L\) denotes the local neural network, with trainable parameters \(\boldsymbol{\theta}_L\), and \(\boldsymbol{\mathcal{I}}n\) is a deterministic interpolation function (e.g. nearest, linear, cubic), which up-samples the predicted average change over subdomain \(\Omega_{\boldsymbol{j}}\) to the mesh resolution. When the local neural network is predicting \(\tilde{\boldsymbol{u}}'_{}\) at a discretized mesh grid, the integral mean in Eq. 8 will degrade to algebraic average. For each subdomain, both the historical information (\(\tilde{\boldsymbol{u}}'_{t_{i-1},\boldsymbol{j}}\)) and its adjacent subdomain information (\(\tilde{\boldsymbol{u}}'(\boldsymbol{x}, t_{i-1})\), \(\boldsymbol{x}\in\Omega^A_{\boldsymbol{j}}\)) will pose influence on the local variance, thus the local model takes the input of both the query subdomain but also its adjacent subdomains \(\Omega^A_{\boldsymbol{j}}\). In one-dimensional cases, \(\Omega_{\boldsymbol{j}}^A\) is defined as \(\{\Omega_{j_1-1}, \Omega_{j_1+1}\}\), while for two-dimensional cases, \(\Omega_{\boldsymbol{j}}^A\) is given by: \(\begin{Bmatrix} \Omega_{j_1-1, j_2-1},& \Omega_{j_1-1, j_2} & \Omega_{j_1-1, j_2+1}\\ \Omega_{j_1, j_2-1},& \Omega_{j_1, j_2} & \Omega_{j_1, j_2+1}\\ \Omega_{j_1+1, j_2-1},& \Omega_{j_1+1, j_2} & \Omega_{j_1+1, j_2+1} \end{Bmatrix}\).

2.6 Long-term Prediction of Spatiotemporal Dynamics at High Resolution↩︎

Once we have the Global Model and Local Model ready, MuRFiV predicts the spatiotemporal dynamics at the original high resolution in an autoregressive manner: \[\begin{align} \tilde{\boldsymbol{u}}(\boldsymbol{x}, t_1) &= \boldsymbol{\mathcal{M}}_G(\boldsymbol{u}_{t_0, \boldsymbol{j}},\; \boldsymbol{\lambda};\;\boldsymbol{\theta}_G) \,+\, \boldsymbol{\mathcal{M}}_L(\boldsymbol{u}_{t_0,\boldsymbol{j}},\;\boldsymbol{\lambda};\;\boldsymbol{\theta}_L)\\ \tilde{\boldsymbol{u}}(\boldsymbol{x}, t_i) &= \boldsymbol{\mathcal{M}}_G(\tilde{\boldsymbol{u}}_{t_{i-1}, \boldsymbol{j}}, \boldsymbol{\lambda}; \boldsymbol{\theta}_G) + \boldsymbol{\mathcal{M}}_L(\tilde{\boldsymbol{u}}_{t_{i-1},\boldsymbol{j}},\boldsymbol{\lambda};\boldsymbol{\theta}_L)\,,\; i\geq2 \end{align}\]

3 Results↩︎

In this section, we demonstrate the performance advantages of two MuRFiV variants in predicting long-term spatiotemporal dynamics, comparing them against U-Net baselines with slightly larger model size. To ensure a fair comparison, all models undergo the same training procedure, consisting of one-step teacher-forcing followed by autoregressive rollout training. Models are evaluated using the checkpoint with the lowest validation loss. Rollout results on training initial conditions and training physical parameters are provided in Appendix 7.

3.1 One-dimensional Burgers’ Equation↩︎

We begin by evaluating MuRFiV on the one-dimensional viscous Burgers’ equation with periodic boundary conditions: \[\frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2}\,,\quad x\in[0,2\pi],\,t\in[0,\infty)\] where \(\nu\) denotes viscosity. The Reynolds number is defined as \(Re = 1/\nu\). The training dataset comprises 48 trajectories across four Reynolds numbers \(\{20, 80, 140, 200\}\), stored on a uniform mesh with \(1024\) cells. Each trajectory starts from an initial condition randomly sampled from a 16-dimensional function space and contains 129 temporal snapshots with a time interval \(\Delta T = 0.01\) seconds. Further details on dataset generation are provided in Appendix 5.1. We train the three models on the same dataset with 100 epochs of one-step teacher-forcing training followed by 50 epochs of eight-setp rollout training. In MuRFiV, we divide the domain into \(32\) subdomains of the same size, each with \(32\) cells.

We then investigate the prediction capability of MuRFiV by evaluating its performance on \(36\) testing trajectories with unseen Reynolds numbers (\(Re\in\{30, 50, 70\}\)) and unseen initial conditions drawn from the same distribution. The test dataset therefore probes both long-term rollout accuracy and parameter generalization beyond the training set.

Figure 2 presents the forecasting results on this extrapolation task. Visualizations of the predicted velocity field for representative low and high Reynolds numbers are shown in Figure 2 (a) and (b). Both MuRFiV variants demonstrate robust generalization, accurately capturing shock formation and propagation despite the shift in physical parameters. In contrast, the U-Net baseline exhibits significant performance degradation, failing to produce physically consistent predictions in the extrapolation regime.

The quantitative superiority of MuRFiV is reinforced by the error analysis in Figure 2 (c-f). Panel (c) reports the average error across all testing cases on a linear scale, while panels (d-f) show the error evolution for specific Reynolds numbers on a log scale. Although U-Net exhibits a lower error \(\epsilon\) than MuRFiV during the initial rollout steps, it quickly accumulates error as autoregressive prediction proceeds. In contrast, both MuRFiV variants maintain stable and accurate predictions over long-term integration, indicating superior robustness against error accumulation. Furthermore, MuRFiV-Eq demonstrates lower prediction error than MuRFiV-NoEq across the unseen Reynolds numbers, underscoring the efficacy of embedding physical operators within the model structure. This enhanced generalizability suggests that the finite-volume backbone and multi-resolution hierarchy enable MuRFiV to learn the underlying physical dynamics rather than merely overfitting to the training data distribution.

3.2 Two-dimensional Burgers’ Equation↩︎

We extend our evaluation to the two-dimensional scalar Burgers’ equation with periodic boundary conditions: \[\frac{\partial u}{\partial t} = -\left(u\frac{\partial u}{\partial x} +u\frac{\partial u}{\partial y}\right)+ \nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right),\; (x, y)\in[0,2\pi]^2,\;t\in[0,\infty) \label{eq:sburgers}\tag{9}\] where \(\nu\) represents the viscosity. Similar to the 1D case, the dataset includes trajectories with varying Reynolds numbers (\(Re=1/\nu \in \{20, 40, 60, 80\}\)), stored on uniform mesh of size \(1024\times 1024\). For each Reynolds number, the training trajectories starts from 36 different ICs randomly sampled from a 72-dimensional distributions and contains 129 temporal snapshots with a time interval \(\Delta T = 0.01\) seconds. Further details on dataset generation are provided in Appendix 5.2. The surrogate models are trained on the same dataset with 100 epochs of one-step teacher-forcing training followed by 50 epochs of 8-step rollout training. In MuRFiV, we divide the domain into \(32\times32\) sub-domains, each consists with \(32\times32\) cells.

We next evaluate the prediction capability of MuRFiV on unseen Reynolds numbers (\(Re\in\{30, 50, 70\}\)) and randomly sampled, unseen initial conditions. The test dataset therefore examines whether the learned dynamics remain accurate when both the physical parameter and the trajectory realization lie outside the training set.

Figure 3 presents the results for this extrapolation task. Visualizations of the predicted scalar field \(u\) for representative low (\(Re=30\)) and high (\(Re=70\)) Reynolds numbers are shown in Figure 3 (a) and (b). These contours show that MuRFiV robustly handles the shift in physical parameters. Even with unseen Reynolds numbers, MuRFiV correctly predicts the shock speed and sharpness, while preserving the complex interference patterns formed by interacting fronts. The U-Net baseline, lacking explicit physical constraints, struggles to generalize and exhibits noticeable phase errors together with overly diffused shock structures.

The error analysis in Figure 3 (c-f) confirms the superior generalization of MuRFiV. Panel (c) plots the average error across all test cases on a linear scale, while panels (d-f) show the error evolution for specific Reynolds numbers on a log scale. These error curves show that both MuRFiV variants deliver much more consistent performance across Reynolds numbers, whereas the baseline degrades substantially outside the training regime. Furthermore, MuRFiV-Eq achieves the lowest error throughout most of the rollout, further highlighting the efficacy of embedding physical operators into the model structure. This enhanced generalizability suggests that the embedded finite-volume structure and multi-resolution hierarchy enable MuRFiV to learn the underlying physical operators rather than simply overfitting to the training data distribution.

3.3 Shallow Water Equation↩︎

We further assess the efficacy of MuRFiV on the two-dimensional shallow water equations, which govern the flow of a thin layer of fluid under the influence of gravity and viscosity. The system is described by: \[\begin{align} \frac{\partial \eta }{\partial t} + \frac{\partial \eta u}{\partial x}+\frac{\partial \eta v}{\partial y} &= 0\\ \frac{\partial \eta u}{\partial t} + \frac{\partial}{\partial x}\left(\eta u^2 + \frac{1}{2}g\eta^2\right) & + \frac{\partial (\eta uv)}{\partial y}= \nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)\\ \frac{\partial \eta v}{\partial t} + \frac{\partial}{\partial y}\left(\eta v^2 + \frac{1}{2}g\eta^2\right) &+ \frac{\partial (\eta uv)}{\partial x}= \nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) \end{align}\,,\, (x, y)\in[0,2\pi]^2,\, t\in[0,\infty)\] where \(\eta\) denotes the fluid column height, and \((u, v)\) are the velocity components. \(g\) is the gravitational acceleration (set to \(10\)), and \(\nu\) represents the viscosity. The training set comprises \(80\) different trajectories initialized from random states, with \(5\) different viscosities (\(1/\nu \in \{1, 2, 4, 8, 16\}\)). Further details on dataset generation are provided in Appendix 5.3. The training data is saved on uniform mesh with \(256\times256\) cells, and we divide the mesh equally into \(16\times16\) subdomins in MuRFiV.

We next evaluate the generalization capability of MuRFiV on \(4\) unseen viscosities (\(1/\nu \in \{1.5, 3, 6, 12\}\)) adn \(64\) randomly sampled initial conditions. The test dataset therefore examines whether the learned dynamics remain reliable when the viscous regime changes beyond the training set.

Figure 4 summarizes the forecasting performance in this extrapolation setting. Panels (a) and (b) visualize representative test cases with the lowest and highest viscosities, respectively. The results indicate that MuRFiV generalizes well to new physical parameters, maintaining physical consistency in wave propagation speeds, amplitudes, and interaction patterns. In contrast, U-Net exhibits significant error growth and unphysical wave damping when applied to out-of-distribution cases.

The error analysis in Figure 4 (c-f) confirms the robustness of MuRFiV. Panel (c) presents the aggregate error on a linear scale, while panels (d-f) depict the error evolution for specific viscosities on a log scale. These plots show that MuRFiV outperforms the baseline by a large margin in this challenging extrapolation setting. Notably, MuRFiV-Eq maintains the lowest error profile throughout most of the rollout, further validating that the embedded conservation structure improves long-term stability under unseen viscosities. Similar to the Burgers’ cases, the results suggest that the physics-informed structural bias enhances the model’s ability to extrapolate to unseen regimes rather than simply memorizing the training distribution.

3.4 Incompressible Navier-Stokes Equation - Kolmogorov Flow↩︎

Finally, we evaluate MuRFiV on the 2D incompressible Navier-Stokes equations configured for Kolmogorov flow, a classical benchmark for turbulence modeling. The governing equations are: \[\begin{align} & \frac{\partial \boldsymbol{u}}{\partial t} = -\left(\boldsymbol{u}\cdot\nabla \right)\boldsymbol{u} + \nu\nabla^2\boldsymbol{u}-\nabla p + \boldsymbol{F}_s\\ & \nabla\cdot \boldsymbol{u} = 0 \end{align}\quad \boldsymbol{x}\in[0,2\pi]^2,\, t\in[0,\infty)\] where \(\boldsymbol{u}\) represents the velocity vector, \(p\) is the pressure field, and \(\nu\) denotes the kinematic viscosity. The flow is driven by a sinusoidal external force \(\boldsymbol{F}_s = \sin(k_f y)\hat{\boldsymbol{x}}\) with forcing wavenumber \(k_f=4\). This setup generates a chaotic flow field characterized by a cascade of energy from large to small scales. The Reynolds number is defined as \(Re = 1/\nu\). The training dataset consists of \(80\) trajectories of fully developed flow across \(5\) Reynolds numbers (\(Re\in\{128, 256, 512, 1024, 2048\}\)), each with \(129\) time steps and saved on a \(256\times256\) mesh. More details on the dataset generation can be found in Appendix 5.4. In MuRFiV, we equally split the \(256\times256\) mesh into \(16\times16\) subdomains, each with a \(16\times16\) mesh.

We next evaluate the generalization capability of MuRFiV on \(4\) unseen Reynolds numbers (\(Re\in\{192, 384, 768, 1536\}\)), starting from \(64\) random initial conditions. The test dataset therefore examines whether the learned dynamics can preserve both instantaneous flow structures and spectral statistics under a viscosity shift.

Figure 5 presents the prediciton results on this extrapolation task. Panels (a) and (b) show instantaneous vorticity fields \(\omega\) (\(\omega = \partial_x v - \partial_y u\)) for representative test cases. These visualizations show that MuRFiV robustly adapts to the change in viscosity, producing physically consistent flow patterns and preserving coherent vortical structures better than the U-Net baseline.

Given the Kolmogorov flow is a chaotic system, we examine the statistics turbulence kinetic energy spectrum instead of absolute error in this section. The spectrum analysis in Figure 5 (c) and (d) further highlights the superiority of MuRFiV. These panels compare the TKE spectra of the predictions against the DNS reference. Even in this extrapolation regime, MuRFiV maintains a TKE spectrum that closely follows the DNS reference, whereas the baseline fails to capture the correct energy distribution, especially at higher wavenumbers, leading to inaccurate representation of the turbulent statistics. Among the two variants, MuRFiV-Eq provides the closest match to the target flow statistics, further demonstrating the benefit of embedding physical structure into the model. This result shows that the finite-volume-inspired multiresolution design improves not only instantaneous prediction quality but also the preservation of turbulent flow statistics under unseen Reynolds numbers.

3.5 Design Insight and Ablation Study↩︎

In this section, we justify the design decision of each MuRFiV component by gradually adding components to the baseline model until reach the full MuRFiV design. All the ablation experiments are tested on the shallow water equation with out-of-distribution initial conditions and viscosity and all the variants has comparable model size (in terms of trainable parameters) as the baseline model (Please refer to Table 4)

Figure 6 shows the box plot of rollout error at different timesteps, compared across all the structures covered by the ablation study, while the detailed error at the last prediction step is reported in Table 1. In this test, we use the same test set used in the Section 3.3.

3.5.0.1 Simple Multi-Resolution Structure Improves long-term prediction

We start improving the baseline method by introducing a multi-resolution structure, noted as SR. SR also consists of two neural networks, one is \(\boldsymbol{\mathcal{N}}_G\) which predicts the solution field at next step on a coarse mesh (\(16\times16\), in this case). The predicted field on the coarse mesh is then up-sampled via bi-cubic to the original resolution and then corrected by a super-resolution network \(\boldsymbol{\mathcal{N}}_\text{SR}\). In Figure 6, the multi-resolution structure (SR) already brings noticeable improvement over the U-Net baseline, especially for long-term rollout predictions. Especially when no numerical operator is embedded (SR-NoEq), we lowered the mean prediction error more than four orders of magnitude evaluated at the \(512^\text{th}\) rollout steps (See Table 1). However, we do observe many “blow-up” cases with equation operators embedded, indicating the instability of this structure when couple with numerical operators.

3.5.0.2 Shared Local Network is Better Than a Global Super-Resolution Network

Next, we leverage the domain decomposition idea used in MuRFiV by replacing the super-resolution network \(\boldsymbol{\mathcal{N}}_\text{SR}\) with local variance predictor \(\boldsymbol{\mathcal{N}}_L\). This variant is noted as LN in Figure 6. Instead of naively send the whole field predicted by \(\boldsymbol{\mathcal{N}}_G\) to \(\boldsymbol{\mathcal{N}}_\text{SR}\) for super-resolution, we decompose the prediction into subdomains which are handled by \(\boldsymbol{\mathcal{N}}_L\) individually. As reported in Table 1, the domain-decomposition structure alone brings over 7 times lower mean prediction error when no numerical operators is embedded (NoEq) and significantly booted the stability (i.e. fewer NaN values) when numerical operators are embedded (Eq).

3.5.0.3 Guaranteed Zero-Mean for Local Network Prediction Improves Accuracy

We further enhance the LN structure by enforcing a zero-mean output for each subdomain (noted as LNSM in Figure 6 and Table 1). Such design choice is inspired by FVM, where the solution in each mesh cell represents the averaged value within the cell. In practice, we guarantee the zero-mean output via subtracting the mean value from the \(\boldsymbol{\mathcal{N}}_L\), so that \(\boldsymbol{\mathcal{N}}_L\) only predicts the variance within each subdomain while \(\boldsymbol{\mathcal{N}}_G\) predicts the mean value of the subdomain. From a neural network perspective, such zero-mean enforcement also serve as a de-bias normalization layer, which also significantly contributes to the long-term rollout stability. The zero-mean/normalization layer design significantly boosted the accuracy over 5 times and eliminated the ‘blow-up’ case when numerical operators are embedded.

3.5.0.4 Predicting Flux Instead of Solutions

Lastly, we make the \(\boldsymbol{\mathcal{N}}_G\) to predict the flux to form the complete structure of MuRFiV. Instead of directly predicting the mean value of each subdomain, \(\boldsymbol{\mathcal{N}}_G\) predicts the flux at each subdomain interfaces, which ensures adjacent subdomains satisfy conservation: the flux leave one domain must go into its adjacent domains¹, while the mean value within the subdomain is then obtained via calculating the net fluxes of the given subdomain. From neural network perspective, such design also serves as a normalization, resembling differential output. The flux prediction structure contribute to another \(23\%\) accuracy improvement over the LNSM when numerical operators are embedded, though the improvement for NoEq variants is relatively minor.

From all the ablation tests, we can verify that all the major design components of MuRFiV improves the long term prediction accuracy, under the same model size.

3.6 Compare with More Baselines↩︎

In this section, we compare MuRFiV with other baseline methods, including ViT[17] and FNO[4]. We also studied MuRFiV with ViT as the major learning structure.

3.6.0.1 Prediction accuracy

Figure 7 (a) shows the prediction error averaged on 64 testing trajectories governed by shallow water equation (the same test set as in Section 3.3). The solid lines indicate the convolution neural network based models, while dashed lines represent the ViT based models. Blue lines are classic models like U-Net and ViT, while the grey line represents FNO. The yellow and red lines represent the MuRFiV-NoEq and MuRFiV-Eq, respectively. The doted blue line represents the error of a larger U-Net, noted as U-Net-L.

On the test set, all the baseline models, including ViT, FNO, and U-Net of two model sizes, show a significantly error accumulation during long term rollout at different pace. While all the MuRFiV variants managed to maintain the prediction error at a low level, even a slight decrease of the mean absolute error can be observed around step 64-256 due to the intrinsic dissipation property of the dynamics itself. Meanwhile, FNO does achieve an impressive lower prediction error for the first few rollout steps on the testing set, and error accumulation rate is relatively slow than other baselines, it still outperformed by all MuRFiV variants by about one order of magnitude lower error at the 512 rollout step.

3.6.0.2 Advantage of MuRFiV is independent of neural network types

Interestingly, the increase of the model size does not help improve the rollout accuracy on the testing set. The U-Net-L even shows a steeper error accumulation curve than the smaller U-Net. Although ViT shows a slower error accumulation rate around step 128-300, it still fall far behind the ViT based MuRFiV models. Despite the ViT and U-Net baseline show noticeable error accumulation pattern due to their significant difference in both the model size and model architecture, the corresponding ViT and CNN based MuRFiV variants show a very similar error curve. This indicate the merit of MuRFiV is highly related to the overall architecture design and is independent of specific neural network formulations.

3.6.0.3 Cost comparison

Figure 7 (b) compares the computational cost of all the models studied in this section. Among the baseline methods, FNO has a model size comparable to ViT and is only moderately larger than MuRFiV-Eq and MuRFiV-NoEq, while its one-step inference cost is substantially higher than the CNN-based MuRFiV variants and the U-Net baselines. Compared with the standard U-Net, the larger U-Net-L also requires noticeably higher computational cost, yet this increase does not translate into better rollout accuracy. The ViT-based models are the most expensive in terms of FLOPs, yet the corresponding ViT-based MuRFiV variants still deliver much better rollout accuracy than the plain ViT baseline. The detailed model sizes (numbers of trainable parameters) are listed in Table 4.

4 Conclusion↩︎

In this work, we introduced MuRFiV, a multi-resolution finite-volume-inspired framework for long-term prediction of spatiotemporal dynamics at the high resolution. The method separates the evolution of subdomain averages from the reconstruction of intra-subdomain structure: a global model updates subdomain averages through interface fluxes, while a lightweight local model shared across subdomains reconstructs zero-mean local fluctuations in a parameter-efficient manner. This decomposition embeds the finite-volume view of conservative dynamics into the prediction architecture and provides a direct mechanism for coupling coarse-scale evolution with fine-scale spatial features.

The framework also supports two complementary variants. In MuRFiV-NoEq, the global network learns interface fluxes directly, showing that the finite-volume-inspired architecture itself provides a strong inductive bias for stable rollout. In MuRFiV-Eq, an embedded numerical operator evaluates flux information on the original fine mesh only along subdomain interfaces, and the global network corrects these fluxes before updating subdomain averages. This interface-only use of the embedded operator preserves high-quality physics-based information at the locations that govern subdomain exchange while avoiding dense high-resolution numerical evaluation over the full domain.

Across one- and two-dimensional Burgers’ equations, shallow water equations, and incompressible Navier-Stokes equations in the Kolmogorov-flow setting, MuRFiV consistently improves long-term autoregressive prediction relative to blackbox data-driven baselines. The results show robust generalization to unseen Reynolds numbers, viscosities, and initial conditions, with MuRFiV maintaining shock propagation, wave interactions, coherent vortical structures, and turbulent kinetic-energy spectra more reliably than the baselines. The ablation study further confirms that the main design choices: multi-resolution coupling, shared local reconstruction, zero-mean local outputs, and flux-based global updates, all contribute to reducing rollout error and improving stability. Additional comparisons with FNO, ViT, and larger U-Net baselines indicate that the advantage of MuRFiV comes primarily from the proposed architecture rather than from a specific neural-network backbone or model size.

The present formulation is especially well aligned with conservation-dominated dynamics in which subdomain-average changes can be represented through flux exchange across interfaces. Future work will extend this idea to more general domain decompositions, complex geometries, non-periodic boundary conditions, and systems requiring additional structure such as strict divergence-free constraints. It will also be valuable to study adaptive subdomain layouts and more expressive embedded operators for problems where elliptic effects, boundary interactions, or multiscale source terms play a dominant role. Overall, the results demonstrate that combining finite-volume-inspired global updates with lightweight local reconstruction provides a practical path toward accurate, stable, and efficient neural surrogates for complex spatiotemporal systems.

Code Availability↩︎

All the source codes to reproduce the results in this study will be openly available on GitHub at https://github.com/jx-wang-s-group/MuRFiV upon publication.

Appendix↩︎

5 Training and Testing Dataset↩︎

All the training and testing datasets are generated by our finite volume solvers based on Jax[18]. In this section, we provide more details about how we set up the numerical solver for generating datasets. For ease of implementation, all the cases apply periodic boundary conditions at all the boundaries.

5.1 One-dimensional Burgers’ Equation↩︎

5.1.0.1 Initial Conditions

To make the dataset challenging and more representative of the dynamics in Burgers’ equation, the training set contains three groups (Group \(D\), \(D^-\), and \(D^+\)) of trajectories and each group contains \(12\) different trajectories with \(4\) different values of Reynolds number \(Re \in \{20, 80, 140, 200\}\). Each trajectory within group \(D\) starts from different initial conditions \(u_{t_0}\), which are randomly sampled from a 16-dimensional random space: \[u_{t_0} = \sum_{k\in\mathcal{K}} A_k\sin\left(k x + \varphi_k\right),\quad \mathcal{K} = \{-4, -3, -2, -1, 1, 2, 3, 4\}\label{eq:ic951db}\tag{10}\] where \(A_k\) and \(\varphi_k\) are random variables sampled from two uniform distribution \(U(0,1)\) and \(U(0,2\pi)\), respectively. Since the initial conditions in group \(D\) are linear combinations of sinusoidal functions, the average values are all strictly zeros, which will eventually develop into a standing wave regardless of the initial condition. To cover a more comprehensive dynamics pattern of Burgers’ equation, we also add groups \(D^-\) and \(D^+\), where the initial conditions are a simple shift of the ICs in group \(D\). Specifically, we generate the ICs in the group \(D^-\) by subtracting \(2\) from ICs in group \(D\) and adding \(2\) to form group \(D^+\)’s ICs. In such a way, the training dataset covers both standing waves and moving waves in both directions in one-dimensional space. Meanwhile, it is still a challenging dataset for neural network training, because merely \(12\) different patterns from a \(16\)-dimensional random space are included in the training set.

The solution is solved with WENO5 scheme and central difference, with first order forward Euler scheme for time marching. The solver used a time step size \(\Delta t = 1\times10^{-4}\)s, and we collect the solution every \(100\) numerical steps. The training and testing set are sampled from the same initial condition distribution but with different random seed and viscosity.

5.2 Two-dimensional Burgers’ Equation↩︎

Two dimensional Burgers’ equation case shares almost the same procedure of generating the initial condition as the one-dimensional Burgers’. The only slight difference exists in how we formulate \(\boldsymbol{u}_{t_0}\): \[u_{t_0} = \sum_{k_j\in\mathcal{K}}\sum_{k_i\in\mathcal{K}} A_{k_i, k_j}\sin\left(k_i x + k_j y + \varphi_{k_i,k_j}\right),\quad \mathcal{K} = \{-3, -2, -1, 1, 2, 3\}\label{eq:ic952db}\tag{11}\] where \(A_{k_i, k_j}\) and \(\varphi_{k_i, k_j}\) are random variables sampled from two uniform distribution \(U(0,1)\) and \(U(0,2\pi)\), respectively. In this case, the initial condition is sampled from a 72-dimensional space. All the other parameters like discretization schemes, time step size are identical as the one-dimensional Burgers’ equation case.

5.3 Shallow Water Equation↩︎

In this case, the initial condition for the fluid column height (\(\eta\)) is also randomly sampled using the same formula as in the two-dimensional Burgers’ equation (Eq.@eq:eq:ic952db ). The only difference is \(\mathcal{K} = \{-2, -1, 0, 1, 2\}\). To ensure we do have a negative fluid column height, the randomly sampled \(\eta\) is then linearly shift so that the minium value of the sampled \(\eta\) is 1. While the velocity is initialized as zero field.

We then use first order schemes for both the spatial and temporal derivatives, with a numerical time step size \(1\times10^{-6}\)s. The data is collected every \(5\times10^3\) numerical steps. We discard the first 16 collection steps, to allow the solution field to develop moving waves.

To resolve the near-discontinuous wave in low viscosity cases while maintain numerical stability in high viscosity cases, the SWE is solved on difference mesh size for different viscosity. The specific mesh size used in preparing the dataset can be found in Table2. The training and testing set are then down-sampled from these fine mesh results to a \(256\time256\) mesh.

5.4 Navier-Stocks Equation↩︎

Similarly, we sample the initial condition for velocity components from a distribution resembles Eq.@eq:eq:ic952db . Here \(\mathcal{K} = \{-1, -1, 0, 1, 1\}\) and \(A_{k_i, k_j}\sim U(0, 0.01)\), both x- and y-directional components are sampled from the same distribution, while the pressure field are initialized as zero. We also add a randomly translation in the external force to enrich the dataset.

We then solve the equation with first order scheme with pressure projection. The numerical solver uses a time step size \(\Delta t = 5\times10^{-6}\) and we only collect data after the turbulence is fully developed, at a frequency of every \(10^4\) numerical steps.

Similar to the shallow water equation, we also use different mesh size under different Reynolds numbers (refer to Table.2 for the exact mesh size used in solvers). Then the results are all down sampled to the \(256\times256\) mesh.

6 Model Size↩︎

7 More Rollout Test Results↩︎

This section summarizes the rollout performance on training initial conditions and training physical parameters. The figures are moved here to keep the main results section focused on generalization to unseen conditions.

References↩︎