Physics-informed traditional and generative Mori-Zwanzig-reduced models for multiscale systems

Introduction

For complex multiscale phenomena the fine spatiotemporal resolution required to resolve all scales of motion is prohibitively expensive. Model reduction is a strategy toward computationally efficient simulations. Projection-based reduced-order models (ROMs) aim to derive modified equations for subsets of the necessary degrees of freedom by projecting the equations onto a basis of reduced dimensionality compared to the original solution space [6, 51, 17]. Proper orthogonal decomposition (POD) is a data-driven technique that uses the singular value decomposition (SVD) of solution snapshots to obtain a basis to be used in a POD-Galerkin method [36, 15]. These methods attempt to capture the essential large-scale dynamics with negligible computational cost relative to full resolution. ROMs can be augmented with a closure model to account for the influence of the truncated modes [5, 1, 2, 21, 10, 54].

The closure model is part of the ROM which accounts for the subgrid-scale (SGS) dynamics without explicitly resolving them. An accurate SGS model is essential because unresolved dynamics can significantly impact the degrees of freedom that are kept. Such effects can be represented in the ROM through memory terms, which account for the history of the interaction between resolved and unresolved degrees of freedom resulting from the application of the Mori-Zwanzig (MZ) formalism from nonequilibrium statistical mechanics [27, 46, 65]. Traditional methods assume that the energy production and dissipation of the small scales are in equilibrium [55, 60]. These deterministic SGS models, based on this equilibrium assumption, can be overly dissipative, leading to saturation of important dynamical features [52, 44, 9]. In some cases, traditional eddy-viscosity-based SGS models can be unstable in numerical simulations [27, 9, 26, 41]. These deficiencies arise in part from an inability to capture the backscatter of energy, especially in regions of chaotic and turbulent flow [38, 43].

Interest in developing data-driven and generative stochastic closure approaches has continued to grow, motivated by the need for efficient climate and weather modeling [28, 49, 4, 14], and developing generative models of underlying dynamics from observations is a rapidly expanding field [59, 24, 30, 56, 32]. There is ongoing work toward developing more accurate MZ-based ROMs using ideas from stochastic modeling [50, 47, 23]. Moreover, stochastic closure approaches can more accurately represent the interactions of the resolved and unresolved scales by accounting for the stochastic effects associated with eliminated degrees of freedom [48, 8, 25]. Statistical and machine-learning-based closure strategies for reduced-order modeling in simulating multiscale or chaotic systems remain an open area of research, particularly toward representing unresolved dynamics accurately while also maintaining stability and computational efficiency at fixed cost [2, 64, 57, 20, 40, 53, 34, 61].

Spectral methods constitute a valuable approach in scientific computing applications toward approximating partial differential equations (PDEs) [58, 7, 11, 29, 45]. The applicability and effectiveness of spectral methods depend crucially on the choice of basis functions used to expand the solution of a PDE [45]. In recent years, physics-informed deep learning has emerged as a contender in providing efficient representations of complex functions for the numerical solution of PDEs [3, 45, 33]. These efforts include identifying custom-made candidate basis functions on which to expand the solution of PDEs for evolution using a spectral approach [45]. Our previous work has demonstrated the benefit of physics-informed operator learning for improving traditional scientific computing methods [62], by using physics-informed custom basis functions in evolving the Markovian dynamics for model reduction. In this work, we further augment the ROM by applying the MZ formalism and developing a closure model for physics-informed custom basis functions.

Technical approach

We focus on the viscous Burgers’ equation, written as

on $(x,t) \in(0,2\pi) \times(0,1)$ with periodic boundary conditions. The initial condition is $s_0 = e^{\sin(x)}$ and the viscosity coefficient is $\nu= 0.1$.

Application of MZ formalism

POD–Galerkin-type ROMs provide the prototypical framework for projection-based model reduction. Galerkin projection enforces the governing equations in a weak sense on a reduced subspace. The resulting ROMs evolve a small set of modal coefficients representing the resolved dynamics.

Applying the Galerkin projection to a time-dependent PDE results in the weak formulation

for $1 \le k \le p$. We assume we have $p$ orthonormal basis functions $\{\widetilde{\phi}_k\}_{k=1}^{p}$ with respect to a suitable inner product $\langle\cdot,\cdot\rangle$ on $(0,2\pi)$. We can approximate the solution of the PDE using a modal expansion

Substituting this expansion into the weak formulation and treating the boundary terms appropriately, as described in Appendix A, yields a system of ODEs for $k=1,\ldots,p$ governing the full-order expansion coefficients

For some cutoff $m<p$, we decompose the full-order coefficients into resolved and unresolved modes as $a=\{\widehat{a},\widetilde{a}\}$, where $\widehat{a}=\{a_k\}_{k=1}^{m}$ and $\widetilde{a}=\{a_k\}_{k=m+1}^{p}$ denote the resolved and unresolved modes, respectively. The goal is to develop a reduced model for the dynamics in the resolved set of variables $\widehat{a}$ for $m<p$ and compensate for the unaccounted information with memory and noise in addition to the reduced Markovian dynamics.

The most optimal projection operator $P$ in an $L^2$ sense is the expectation conditioned on the resolved variables with respect to a suitable weight function. We take this to be the product

where $w(\widehat{a})$ is an arbitrary distribution for the $\widehat{a}$’s and $\delta(\widetilde{a})$ is the Dirac measure for the $\widetilde{a}$. We specifically choose $w$ as the Gaussian measure, defined by

For viscous Burgers’, since all non-constant modes eventually converge to zero, we take $\mu_k=0$ and $\sigma_k$ to be rapidly decaying, e.g., $\sigma_k=e^{-k}$ for all $k$. We note that this conditional expectation is equivalent to setting the unresolved modes to zero since

The projection operator $P$ can also be viewed as an orthogonal projection, with respect to the inner product $\langle\cdot,\cdot\rangle_C$, on the subspace of all functions of $\widehat{x}$ defined by

We can approximate this integral using

where $\{(z_i,\omega_i)\}_{1\leq i\leq N_q}$ is the Gauss-Hermite quadrature rule for $N_q$ quadrature points.

After projecting a high-dimensional dynamical system onto resolved coordinates with initial conditions $\varepsilon$, the exact reduced dynamics can be expressed as a generalized Langevin equation using the MZ formalism

where $R_k$ is the Markovian term, $K_k$ accounts for memory effects, and $\zeta_k$ represents the influence of unresolved initial conditions via the orthogonal dynamics [27, 46, 65].

Applying the projection again eliminates the noise, since it is equivalent to only considering initial conditions of the form $\varepsilon=(\hat{\varepsilon},0)$. The projection operator commutes with nonlinear functions, and so the MZ identity for a resolved variable $a_k(t)$ with initial conditions $a(0)=(\hat{\varepsilon},0)$ is

for $k=1,\ldots,m$. The Markovian dynamics are defined by

The SGS term is therefore the convolution integral representing the memory effects

where we have defined

The technical approach is summarized in the schematic in Figure 1 and given in more detail in Appendix § 7. In this work, we demonstrate stochastic generative modeling approaches toward representing this history-dependent term.

Generalized Langevin dynamics

Our goal is to develop models that can accurately evolve a range of initial conditions. Following [18, 19, 21], we can expand the memory integrand in terms of kernel functions: if $\{h_l(\hat{x})\}_{l\in I}$ is a basis for functions of the resolved variables and $I$ is the corresponding index set of the chosen basis, then we can write

for $k=1,\ldots,m$, where $\{\kappa_{lk}(t)\}$ are termed the memory kernels. The calculation of these functions is carried out by leveraging the fluctuation-dissipation theorem (FDT) [35, 12, 39]. In practice, we can only use finitely many basis functions so the above expansion is necessarily replaced by an approximation. The technical approach is summarized in the schematic in Figure 2.

Traditional MZ-based closure model

Memory kernels

Following [19], the memory kernels obey

where $f_{lk}$ and $g_{lj}$ are auxiliary quantities defined, for $l,j\in I$ and $1\leq k\leq m$, as

and where $\mathcal{L}=\sum_{j=1}^{p}R_j(\varepsilon)\frac{\partial}{\partial\varepsilon_j}$ is the Liouville operator. Denoting by $\{H_l(z)\}_{0\leq l\leq d}$ the one-dimensional Hermite polynomials corresponding to weight functions $e^{-z^2/2}$, we take $\{h_j(\hat{\varepsilon})\}_{j\in I}$ to be the tensor products of $H_l\left(\frac{\varepsilon_k}{\sigma_k}\right)$. In order to establish a proof of concept, we limit ourselves to $d=1$, so we have $|I|=\binom{m+d}{d}=m+1$, where $I$ denote the index set of the multivariate Hermite basis functions.

The inner products in (18) are approximated using (10), with the quadrature points furnishing the initial conditions for generating the full-order solutions. More precisely, these are of the form

We can compute the physical initial conditions that correspond to these modes via

The collection of initial conditions for $N_q = 2$ quadrature points is shown in Figure 3.

Figure 4 plots auxiliary quantities $f(t) = (f_{lk}(t))_{l\in I,1\le k\le m}$ and $g(t) = (g_{lj}(t))_{l,j\in I}$. Figure 5 plots the memory kernels $K(t) = (\kappa_{lk}(t))_{l\in I,1\le k\le m}$, calculated by discretizing the convolution integral using the trapezoidal rule. See Appendix § 8 for full derivations.

Reduced model

For $1 \le l \le m$, we model the reduced coefficients $u_k$ and incorporate the memory kernel

for $k=1,\ldots,m$. The system is integrated forward in time using an explicit trapezoidal rule via Heun’s method.

where

Then, the reduced field $\hat{s}$ is reconstructed as

Figure 6 plots the evolution of the true full-order coefficients $a$, Markovian coefficients $b$, and reduced coefficients $u$. The reduced dynamics better match the true coefficients than the Markovian dynamics alone. Figure 7 shows the reduced and Markovian fields $\hat{s}$ with absolute errors compared to the true filtered field. We see smaller absolute error magnitudes using the memory kernels compared to the Markovian dynamics alone. We also see increased error towards the end of the temporal domain. These results show that for longer time, we may require higher-order projections (i.e., beyond linear) for representing the memory kernels. We also require a much finer timestep when discretizing the convolution integral, typically smaller by a factor of 100 compared to the full order and Markovian models.

Figure 8 plots the relative $\ell_2$ error of the reduced and Markovian approximations compared to the true filtered field. We see a decreased error across the temporal domain using the memory kernels, except for longer time. Table 1 shows the average relative $\ell_2$ error, only 2.9% for the reduced model compared to almost 7% for the Markovian model.

	Average relative $\ell_2$ error
Markovian	6.6%
Reduced	2.9%

Generative closure model

Diffusion modeling

A generative approach samples the SGS term from the data distribution $p_{\mathrm{data}}$ along the trajectory in time via learning a denoising diffusion process [31, 56, 30]. We assume access to some initial distribution $p_{\mathrm{init}}$ that we can easily sample from, such as the Gaussian $p_{\mathrm{init}} = \mathcal{N}(0, I_d)$. The goal of our model is then to transform samples from $t_0 \in p_{\mathrm{init}}$ into samples from $t_1 \sim p_{\mathrm{data}}$. We assume access to finite independent samples from $p_{\mathrm{data}}$, which serve as a proxy for the true distribution. The simulation of a stochastic differential equation (SDE) can be used for converting a simple distribution $p_{\mathrm{init}}$ into a more complex target distribution $p_{\mathrm{data}}$. The diffusion model consists of a flow model $u^\theta: \mathbb{R}^d \times[0,1] \to\mathbb{R}^d$, $(t,\gamma) \mapsto u^\theta_\gamma(t)$ with parameters $\theta$ and a fixed diffusion coefficient $\sigma_\gamma: [0,1] \to[0,\infty)$, $\gamma\mapsto\sigma_\gamma$, such that

If we set $\sigma_\gamma= 0$, we recover a flow model. The flow training target $u_\gamma^{\mathrm{target}}$ is the marginal vector field. To construct a training target $u_\gamma^{\mathrm{target}}$ for a flow model, we choose a conditional probability path $p_\gamma(t|\tau)$ that fulfills $p_0(\cdot|\tau) = p_{\mathrm{init}}$, $p_1(\cdot|\tau) = \delta_\tau$ (Dirac delta distribution). The probability path specifies a gradual interpolation between $p_{init}$ and $p_{data}$. We consider a Gaussian conditional probability path $p_\gamma(t|\tau)=\mathcal{N}(t;\alpha_\gamma\tau,\beta_\gamma^2 I_d)$. This is the probability path used by denoising diffusion models [30, 31]. Let $\alpha_\gamma:[0,1]\to\mathbb{R}$ and $\beta_\gamma:[0,1]\to\mathbb{R}$ be noise schedulers defined as monotonic, continuously differentiable functions satisfying $\alpha_1=\beta_0=1$ and $\alpha_0=\beta_1=0$. Flow matching consists of training a model $u^\theta_\gamma$ via minimizing the conditional flow matching (CFM) loss

where $u^{\mathrm{target}}_\gamma(t|\tau)$ is the conditional vector field. This can be any vector field such that the corresponding ODE gives the conditional probability path, i.e.,

for $0\le\gamma\le1$. For a Gaussian probability path, this expression can be found analytically as

After training, samples are generated by simulating the SDE with stochastic dynamics at each physical timestep for some $\sigma_\gamma$. This form has the same marginal probability path as the corresponding SDE. The score correction ensures that the learned distribution has the correct marginal density evolution toward $p_{data}$ [56, 30, 42]. We can write the form for the diffusion model as

where we train a score network $s^\theta_\gamma$ via denoising score matching (DSM)

where the conditional score distribution is

Further training details and hyperparameters for flow and score networks are in Appendix § 9.

At inference time, we can combine $s^\theta_\gamma(t)$ with $u^\theta_\gamma(t)$ and sample via simulating the SDE from $\gamma=0$ to $1$. Often, our goal is to generate samples conditioned on some data $S$. Guided (conditional) generation involves sampling from $\tau\sim p_{data}(\cdot|S)$, where $S$ is a conditioning variable and $p_{data}(\cdot|S)$ is the conditional data distribution. Guided flow matching consists of training a model $u^\theta_\gamma$ via minimizing the classifier-free guidance (CFG) CFM loss

To incorporate score matching, we additionally train a score network $s^\theta_\gamma$ by minimizing the CFG DSM loss

These objectives account for the possibility of the “label” not being in the training set (i.e., $S=\varnothing$). At inference time, for a fixed choice of guidance scale, we can combine $s^\theta_\gamma(t|S)$ with a guided vector field $u^\theta_\gamma(t|S)$ and define

Then we can sample via

simulating the SDE from $\gamma=0$ to 1. The goal is for $t_1$ to adhere to the guiding variable $S$.

Generative model for memory effects

Next, we demonstrate the applicability of generative SGS modeling for the memory term directly, improving upon traditional analytical approaches. We aim to develop a reduced model for the coefficients $r_k$ to approximate the true resolved coefficients $a_k$ for $k=1,\ldots,m$. The reduced dynamics are modeled as

We choose to represent the memory effects using a conditional distribution

This approach differs from the previous formulation in Eq. (15), where we explicitly construct intrusive memory kernels and then evaluate the memory integral. We use diffusion modeling to learn the conditional distribution $p_{\mathrm{data}}(\cdot|\hat{r})$, i.e., the resolved modes are conditioned on the complete set of resolved variables. This conditioning adheres with the Markovian dynamics, which are dependent on the set of resolved variables. Then, the reduced approximation $\hat{s}$ is compared to the filtered full-order field $\bar{s}$ (i.e., using only $m$ modes in the reconstruction).

The technical approach is summarized in the schematic in Figure 9.

For the reduced model, we take $m = 4$, compared to the full-order system with $p = 67$. The generative SGS terms conditioned on the test data compared to the true SGS terms are plotted in Figure 10. With the incorporation of stochastic dynamics through score matching and any $\sigma_\gamma> 0$, the resulting predictions are expected to exhibit additional variability relative to flow matching alone. For these experiments, we empirically select $\sigma_\gamma= 0.1$, which controls the magnitude of the injected diffusion and therefore the stochasticity of the predictions.

For deploying the reduced model, we use a Runge-Kutta (RK) time integration scheme at each time $t_n$, with a factor of 20 fewer timesteps compared to the full-order simulation. We sample $\tau^{\mathrm{SGS}}(\hat{r}^n,\xi) \sim p_{\mathrm{data}}(\cdot\mid\hat{r}^n)$ by simulating the SDE and then integrate forward as

where

For fixed function evaluation cost, we use 20 diffusion timesteps for the SDE integration at each physical timestep. Figure 11 plots the true full-order coefficients $a$, Markovian coefficients $b$, and reduced coefficients $r$. The reduced dynamics better match the true coefficients than the Markovian dynamics alone. Figure 12 shows the reduced and Markovian fields $\hat{s}$ with absolute errors compared to the true filtered field. We see smaller absolute error magnitudes using the generative SGS model compared to the Markovian dynamics alone.

Figure 13 plots the relative $\ell_2$ error of the reduced and Markovian approximations compared to the true filtered field. For this realization, we see a decreased error across the temporal domain using the generative SGS model. Table 2 shows the average relative $\ell_2$ error, less than 2.5% for the reduced model compared to almost 7% for the Markovian model across 10 samples.

	Average relative $\ell_2$ error
Markovian	6.6%
Reduced	2.4%

Conclusions

In this work, we built reduced models by applying the MZ formalism to the system of ODEs for the physics-informed expansion coefficients. We implemented the FDT to solve for memory kernels applicable to evolving a collection of initial conditions and demonstrate lower errors on average for evolution of an out-of-distribution initial condition. These results show that higher-order approximations for the memory kernels become relevant and may be required to continue accounting for accurate memory for longer times. So, future work might include developing different representations for the memory kernels themselves, through generative approaches and higher-order projections. Finally, we used generative modeling to learn a conditional distribution to represent the memory effects on the resolved variables. We see robust predictions for evolution of a particular initial condition, with lower errors than evolving the Markovian dynamics alone. We see benefits in representing the memory itself stochastically. Future work should include accounting for memory by conditioning on previous states for some memory horizon length, as well as training across a collection of initial conditions.

Acknowledgements

This project was completed with support from the U.S. Department of Energy, Advanced Scientific Computing Research program, under the Scalable, Efficient and Accelerated Causal Reasoning Operators, Graphs and Spikes for Earth and Embedded Systems (SEA-CROGS) project (Project No. 80278). The computational work was performed using PNNL Institutional Computing at Pacific Northwest National Laboratory. Pacific Northwest National Laboratory (PNNL) is a multi-program national laboratory operated for the U.S. Department of Energy (DOE) by Battelle Memorial Institute under Contract No. DE-AC05-76RL01830.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC0023112. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Weak formulation derivation

Rewrite the viscous Burgers’ equation as

Set $f(s)=s^2/2$ and substitute $f(s)$ into (57) to write

Rewriting as a system of first-order equations gives

For the weak formulation, we can write

Integrating by parts on Eq. (61) gives

Writing the flux terms as surface integrals gives

where $\hat{\mathbf{n}}$ is the outward-facing normal, $f^*$ represents the numerical flux term [45], and the subscripts $R$ and $L$ refer to values of the functions at the right and left boundaries. The numerical flux $f^*(c^-,c^+)$ used is defined by

Integrating by parts on Eq. (62) gives

Discretize the approximations using linear expansions in terms of the basis functions

Substituting the summations into Eq. (65) to write

where $x_{out}=2\pi$ and $x_{in}=0$. Expanding the surface integrals gives

As a result, the expansion coefficients for Eq. (61) are given by

for $1\leq k\leq p$. Substituting the summations into Eq. (70)

Expanding the surface integrals yields

So the expansion coefficients for Eq. (62) are given by

for $1\leq k\leq p$. Substituting Eq. (77) into Eq. (74) yields a system of ODEs for the expansion coefficients

with the initial condition

for $1\leq k\leq p$. We use a quadrature rule to calculate the inner products.

Let $\mathcal{A}_k[\hat{s}_p]$ represent the nonlinear terms, $\mathcal{D}$ represent the viscous terms, and $\mathcal{Q}$ represent the derivative terms, allowing us to write

where

The Markovian term is given by

with the boundary conditions incorporated as described above. Using $\mathcal{A}$, $\mathcal{D}$, and $\mathcal{Q}$, we can write

for $1\leq k\leq m$.

Fluctuations and kernels for MZ-reduced model

From Appendix A, we have the system of ODEs for the expansion coefficients

with the initial condition

where

For $m < p$, we want to build a ROM based on only $\hat{a} = \{a_k\}_{k=1}^{m}$. With initial condition $a(0) = \{\{\hat{a}_k\}_{k=1}^{m},\mathbf{0}_{k=m+1}^{p}\}$. The MZ identity is

where $\{h_j(\hat{\varepsilon})\}_{j\in I}$ is a basis for functions of the resolved variables and $\{\kappa_{jk}(t)\}_{j\in I,1\leq k\leq m}$ are the memory kernels. To calculate the memory kernels [19], we solve

where

where

Denoting by $\{H_l(z)\}_{l\geq0}$ the one-dimensional Hermite polynomials corresponding to weight functions $e^{-z^2/2}$, we take $\{h_j(\hat{\varepsilon})\}_{j\in I}$ to be their tensor products. We consider the case of linear polynomials $d=1$, so have $|I| = \binom{m+d}{d}=m+1$. Write out the terms in $g_{lj}$

Write out the terms in $f_{lk}$

Then we have

with

where

Then we can find $f_{lk}$ by computing the inner product with $\{h_l\}$. To calculate the memory kernels $K(t)=(\kappa_{lk}(t))_{l\in I,1\le k\le m}$, discretize the convolution integral using the trapezoidal rule. At each time $t_n$, we have

with $K(t_0)=f(t_0)$.

Generative model training details and hyperparameters

For generative models, such as diffusion models, future work will include more rigorous assessment of hyperparameter settings, including network sizes, learning rates, and training iterations, as well as inference parameters during deployment, such as guidance scales and noise constants, as in [13, 16, 63, 22, 37]. Future studies will investigate implications of parameter choices and sensitivities, as well as assessing optimal tuning strategies for training and inference hyperparameters.

	Flow	Score
Network size	$3 \times128$	$4 \times256$
Learning rate	0.001	0.0001
Activation function	SiLU	SiLU
Training epochs	5000	5000
Batch size	100	50
Unconditional training probability	0.1	0.2
Guidance scale	1.5	0.8