Lagrangian Dynamic Mode Decomposition for Construction of Reduced-Order Models of Advection-Dominated Phenomena

TL;DR

This paper introduces a Lagrangian-based dynamic mode decomposition method to improve reduced-order modeling of advection-dominated phenomena, addressing limitations of traditional Eulerian POD and DMD techniques.

Contribution

The paper proposes a novel Lagrangian DMD approach that better captures advection effects, enhancing accuracy and efficiency over conventional methods.

Findings

Lagrangian DMD outperforms Eulerian methods in advection problems.

The new method demonstrates improved accuracy in numerical tests.

Lagrangian DMD is computationally efficient for complex flows.

Abstract

Proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD) are two complementary singular-value decomposition (SVD) techniques that are widely used to construct reduced-order models (ROMs) in a variety of fields of science and engineering. Despite their popularity, both DMD and POD struggle to formulate accurate ROMs for advection-dominated problems because of the nature of SVD-based methods. We investigate this shortcoming of conventional POD and DMD methods formulated within the Eulerian framework. Then we propose a Lagrangian-based DMD method to overcome this so-called translational issues. Our approach is consistent with the spirit of physics-aware DMD since it accounts for the evolution of characteristic lines. Several numerical tests are presented to demonstrate the accuracy and efficiency of the proposed Lagrangian DMD method.

Tables

Table 1. Table 1: Computational cost of the high-fidelity models and the corresponding Lagrangian DMD- and POD-based reduced-order models. Test 1 refers to advection problem; Test 2 to advection-diffusion problem; and Tests 3 and 4 to inviscid and viscous Burgers equations, respectively.

	Test 1	Test 2	Test 3	Test 4
Rank truncation $r$	3	10	3	14
DMD computational time (sec)	0.114718	0.046531	0.048450	0.055641
POD computational time (sec)	0.153869	0.320905	0.435342	0.795566
Eulerian HFM computational time (sec)	1.390459	29.782079	0.034519	55.713132
Lagrangian HFM computational time (sec)	0.023568	27.020414	0.039063	55.246262

Equations

$$\frac{\partial u}{\partial t}+f(u)\frac{\partial u}{\partial x}=\frac{\partial }{\partial x}\left(D(x,t,u)\frac{\partial u}{\partial x}\right),f(u)=\frac{\partial F(u)}{\partial u},$$

$$u(x,t=0)=u_0(x)$$

$$u_j^{n+1}=u_j^n-\frac{\Delta t}{\Delta x}(F_{j+1/2}^n-F_{j-1/2}^n)+\frac{\Delta t}{(\Delta x{)}^2}\left[D_{j+1/2}^{n+1}(u_{j+1}^{n+1}-u_j^{n+1})-D_{j-1/2}^{n+1}(u_j^{n+1}-u_{j-1}^{n+1})\right],$$

$$F_{j+1/2}^n=\frac{F(u_{j+1}^n)+F(u_j^n)}{2}-|a_{j+1/2}^n|\frac{u_{j+1}^n-u_j^n}{2},$$

$$a_{j+1/2}^n=\{\begin{array}{rlrlrl}& \frac{F_{j+1}^n-F_j^n}{u_{j+1}^n-u_j^n}& & \text{mbox}if& & u_{j+1}^n\ne u_j^n,\\ & f(u_j)& & \text{mbox}if& & u_{j+1}^n=u_j^n,\end{array}$$

$$D_{j+1/2}^{n+1}=\frac{1}{2}(D_j^{n+1}+D_{j+1}^{n+1}).$$

$$\mathbb{R}(u^{n+1})=u^{n+1}-u^n+\Delta t({\mathcal{D}}_1^u{\mathbb{F}}^n)-\Delta t({\mathcal{D}}_2{u}^{n+1})=0,$$

$$\mathbb{X}=\left[\begin{array}{cccc}|& |& & |\\ u^1& u^2& \cdots & u^m\\ |& |& & |\end{array}\right],\mathbf{X}\in {\mathbb{R}}^{N\times m}.$$

$$\mathbb{X}=\mathbb{U}\Sigma {\mathbb{V}}^{\ast },$$

$$r=\underset{k}{\min }\left\{\frac{{\sigma }_k}{{\sum }_{k=1}^K{\sigma }_k}<\epsilon \right\},$$

$$\mathbf{\Phi }=\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{\varphi }}_1& {\mathbf{\varphi }}_2& \cdots & {\mathbf{\varphi }}_r\\ |& |& & |\end{array}\right].$$

$$u_{\text{POD}}^{n+1}=\sum _{k=1}^r{\hat{u}}_k^{n+1}{\mathbf{\varphi }}_k=\mathbf{\Phi }{\hat{u}}^{n+1}.$$

$${\mathbf{\Phi }}^{\top }\mathbb{R}(\mathbf{\Phi }{\hat{u}}^{n+1})=0.$$

$$\frac{\text{d}u}{\text{d}t}=\mathcal{N}(u),$$

$${\mathcal{N}}_t(u(t_0))=u(t_0+t)=u(t_0)+{\int }_{t_0}^{t_0+t}\mathcal{N}(u(\tau ))d\tau ,$$

$$u^{n+1}={\mathcal{N}}_t(u^n).$$

$$\mathcal{K}g(u)=g(\mathcal{N}(u)).$$

$${\mathcal{K}}_{t}g(u^n)=g({\mathcal{N}}_t(u^n))=g(u^{n+1}).$$

$${\mathbb{X}}_1=\left[\begin{array}{cccc}|& |& & |\\ u^1& u^2& \cdots & u^{m-1}\\ |& |& & |\end{array}\right],{\mathbb{X}}_2=\left[\begin{array}{cccc}|& |& & |\\ u^2& u^3& \cdots & u^m\\ |& |& & |\end{array}\right],$$

$$u_{\text{DMD}}^{n+1}=\mathbf{\Phi }{\mathbf{\Lambda }}^{n+1}b,n>m,$$

$$\left\{\begin{aligned} &\frac{\text{d}X(t)}{\text{d}t}=f(u(X(t),t)),\\ &\frac{\text{d}u(x,t)}{\text{d}t}\bigg{\rvert}_{x=X(t)}=\left[\frac{\partial}{\partial x}\left(D(x,t,u)\frac{\partial u(x,t)}{\partial x}\right)\right]\bigg{\rvert}_{x=X(t)}.\end{aligned}\right.$$

$$\{\begin{array}{rl}& {\tilde{u}}_j^n={\mathcal{P}}_0(u_j^n),\\ & {\tilde{u}}_j^{n+1}={\tilde{u}}_j^n+\Delta t\frac{D_{j+1/2}^{n+1}({\tilde{u}}_{j+1}^{n+1}-{\tilde{u}}_j^{n+1})-D_{j-1/2}^{n+1}({\tilde{u}}_j^{n+1}-{\tilde{u}}_{j-1}^{n+1})}{(\Delta x{)}^2},\\ & u_j^{n+1}={\mathcal{P}}_n({\tilde{u}}_j^{n+1}),\\ & x_j^{n+1}=x_j^n+\frac{\Delta t}{2}(f(u_j^n)+f(u_j^{n+1})),\end{array}$$

$$\{\begin{array}{rl}& {\mathbb{R}}_x(x^{n+1})\equiv x^{n+1}-x^n-\frac{\Delta t}{2}(f(u^n)+f(u^{n+1}))=0,\\ & {\mathbb{R}}_u({\tilde{u}}^{n+1})\equiv {\tilde{u}}^{n+1}-{\tilde{u}}^n-\Delta t{\mathcal{D}}_2{\tilde{u}}^{n+1}=0.\end{array}$$

$$\mathbb{X}=\left[\begin{array}{cccc}|& |& & |\\ x^1& x^2& \cdots & x^m\\ |& |& & |\\ u^1& u^2& \cdots & u^m\\ |& |& & |\end{array}\right].$$

$$\left[\begin{array}{c}|\\ x_{\text{POD}}^{n+1}\\ |\\ u_{\text{POD}}^{n+1}\\ |\end{array}\right]=\mathbf{\Phi }\left[\begin{array}{c}|\\ {\hat{x}}^{n+1}\\ |\\ {\hat{u}}^{n+1}\\ |\end{array}\right].$$

$${\mathbf{\Phi }}^T\mathbb{R}\left(\mathbf{\Phi }\left[\begin{array}{c}|\\ {\hat{x}}^{n+1}\\ |\\ {\hat{u}}^{n+1}\\ |\end{array}\right]\right)=0.$$

$$\mathbb{y}^{n}=\mathbb{g}(\mathbb{u}^{n})=\begin{bmatrix}g_{1}(\mathbb{u}^{n})\\ g_{2}(\mathbb{u}^{n})\\ \vdots\\ g_{p}(\mathbb{u}^{n})\end{bmatrix},\qquad g_{j}:\mathcal{M}\to\mathbb{C}\ \mbox{are observable function ($j=1,2,\cdots,p$)}.$$

$${\mathbb{Y}}_1=\left[\begin{array}{cccc}|& |& & |\\ y^1& y^2& \cdots & y^{m-1}\\ |& |& & |\end{array}\right],{\mathbb{Y}}_2=\left[\begin{array}{cccc}|& |& & |\\ y^2& y^3& \cdots & y^m\\ |& |& & |\end{array}\right],$$

$$y_{\text{DMD}}^{n+1}=\mathbf{\Phi }{\mathbf{\Lambda }}^{n+1}b,n>m$$

$$u_{\text{DMD}}^n=g^{-1}(y_{\text{DMD}}^n).$$

$$y^n=g(u^n)=\left[\begin{array}{c}{g}_1^n\\ g_2^n\end{array}\right],\text{mbox}with{g}_1=x^n,g_2=u^n.$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

Advection-diffusion equations are routinely used as a high-fidelity representation of mass conservation at a variety of spatiotemporal scales in a plethora of applications [1]. These equations become highly nonlinear when advection velocity and/or diffusion coefficient depend(s) on a system state, e.g., in the case of multiphase flows in porous media [1]. High-dimensional complex dynamics described by such nonlinear advection-diffusion equations often posses low-dimensional structures, suggesting the possibility of their replacement with reduced-order models (ROMs) [2, 3, 4].

Singular-value decomposition (SVD) can be utilized to extract a low-dimensional structure from the data generated with a high-fidelity model (HFM), i.e., to construct a conventional ROM. An example of such a ROM is built by deploying Galerkin projection to map a HFM onto a much smaller subspace; projection-based ROMs are referred to as proper orthogonal decomposition (POD) [5, 6, 7]. The efficiency and accuracy of POD in nonlinear setting are increased by combining it with either the empirical interpolation method (EIM) [8] or discrete empirical interpolation method (DEIM) [9]. Another example is the Dynamic Mode Decomposition (DMD) method [10, 11], which is used to discover a spatiotemporal coherent structure in the HMF-generated data. DMD’s connection to the Koopman operator theory of nonlinear dynamic systems [12] is of theoretical interest [13, 11, 14, 15], while its equation-free spirit facilitates its use in conjunction with machine learning techniques [16, 17].

While the robustness of DMD for parabolic problems has been established (including numerical [18] and theoretical  [19, 20] analysis of its accuracy and convergence), both DMD and POD are known to fail in translational problems, such as wave-like phenomena, moving interfaces and moving shocks [10].It can be explained by the intuition that the dominating advection behavior is traveling through the whole high-dimensional domain, making it impossible to determine a global spatiotemporal basis confined in a low-dimensional subspace. We facilitate this intuitive explaination with a concrete example in section 2. In terms of the Koopman operator theory, important physical observables (e.g., advection speed, shock speed, shock formation time) are unaccounted for in the standard DMD algorithm. Remedies for POD include the deployment of local basis [21], domain decomposition [22], or basis splitting [23]. A similar extension of DMD consists of multi-resolution DMD [24], which separates frequencies of different scales by filtering windows. Unfortunately, these remedies often compromise the ROM’s efficiency by increasing its computational complexity. Alternative generalizations of DMD and POD explore symmetry and self-similarity properties to eliminate the translational issue using analytical tools [25, 26, 6, 27, 28, 29]. However, such tools are usually problem-dependent and mostly applicable to single-wave dominated problems.

Motivated by the recent work on Lagrangian POD [30], we propose a physics-aware DMD method to construct a ROM within the Lagrangian framework. We chose the temporally evolving characteristic lines, a crucial physical quantity in advection-dominated systems, as a key observable of the underlying Koopman operator. Then, the DMD algorithm is used to identify, from sufficient data, a low-dimensional structure in the Lagrangian framework and thus to construct a physics-based ROM by approximating the underlying Koopman operator. The Lagrangian DMD can be applied to general advection-diffusion nonlinear flows. Furthermore, DMD outperforms POD in terms of computational costs due to the feature of iteration free. With the error analysis in [20], one can also estimate the accuracy of the ROM.

In Section 2, we review conventional POD and DMD in the Eulerian framework and demonstrate the issue of translational problems using advection-dominated linear flow as an example. Section 3 provides a brief illustration of the Lagrangian POD in [30], introduces our new Lagrangian DMD, and demonstrates its connection to the Koopman operator. Section 4 contains several computational experiments used to validate the accuracy and efficiency of the proposed approach to ROM construction. It also compares the Lagrangian POD and DMD in terms of their accuracy and computational costs. Finally, in Section 5, we draw conclusions and discuss the related ongoing work.

2 Conventional Eulerian Reduced-Order Models

Consider a scalar state variable $u(x,t):[a,b]\times[0,T]\rightarrow\mathbb{R}^{+}$, whose dynamics is described by a one-dimensional nonlinear advection-diffusion equation

[TABLE]

subject to the initial condition

[TABLE]

and appropriate (arbitrary) boundary conditions at $x=a$ and $x=b$. Within the Eulerian framework, the space is fixed and the interval $[a,b]$ is discretized with a uniform grid $\mathbb{x}=[x_{1}=a,x_{2}\cdots,x_{N-1},x_{N}=b]^{\top}$ of mesh size $\Delta x\equiv x_{j+1}-x_{j}=(b-a)/N$ and $N$ nodes. Likewise, the time interval $[0,T]$ is discretized uniformly with time step $\Delta t\equiv t^{n+1}-t^{n}=T/M$ and $M+1$ nodes so that $t^{0}=0<t^{1}<\cdots<t^{M}=T$. At the $n$th time node, the state variable $u(x,t)$ is represented by a vector $\mathbb{u}^{n}=[u_{1}^{n},\cdots,u_{N}^{n}]^{\top}$ for $n=0,\cdots,M$. For simplicity, (2.1) is solved with a conservative first-order upwind scheme with forward Euler for the advection part and center difference with backward Euler for the diffusion part,

[TABLE]

where

[TABLE]

In vector form, the above scheme reads

[TABLE]

where $\mathcal{D}_{1}^{u}\in\mathbb{R}^{N\times N}$ and $\mathcal{D}_{2}\in\mathbb{R}^{N\times N}$ are discrete approximations of the first derivative (using upwind) and second derivative (using center difference), respectively. Here $\mathbb{F}^{n}=[F_{1/2}^{n},\cdots,F_{N-1/2}^{n}]^{\top}$ and $\mathbb{R}$ is the vectorized residual of the scheme. Certain CFL condition needs to be satisfied to ensure the stability of the scheme depending on the functional forms of $f$ and $D$. Simulation results obtained with the above method constitute our high-fidelity model (HFM).

A reduced-order, low-fidelity model (ROM) is constructed from a data set comprising a sequence of solution snapshots collected from the HFM. Let $\mathbb{X}$ denote the data matrix, consisting of $m$ snapshots of $\mathbf{u}$,

[TABLE]

Two alternative strategies for building a ROM from these data, both grounded in Singular Value Decomposition (SVD), are described below.

2.1 POD

The conventional POD method generates a ROM by using a low-dimensional basis extracted from the data $\mathbf{X}$ in (2.5) to project the dynamics $\mathbf{u}(t)$ onto a lower-dimensional hyperplane. If the data matrix $\mathbf{X}\in\mathbb{R}^{N\times m}$ has rank $K\leq\min\{N,m\}$, then the POD modes are constructed by using a reduced SVD,

[TABLE]

where $\mathbb{U}\in\mathbb{C}^{N\times K}$ is the matrix of $K$ orthonormal columns of length $N$; $\mathbb{\Sigma}=\mathbb{R}^{K\times K}$ is the diagonal matrix with real diagonal elements $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{K}>0$; $\mathbb{V}\in\mathbb{C}^{m\times K}$ is the matrix of $K$ orthonormal columns of length $m$; and the superscript ∗ denotes its conjugate transpose. A reduced-order model is constructed by choosing a rank $r$ ($r\ll K$), which satisfies the energy criteria

[TABLE]

where $\varepsilon$ is a user-supplied small number ($\varepsilon=10^{-4}$ in all our numerical examples). Next, the matrix $\mathbb{U}\in\mathbb{C}^{N\times K}$ is replaced with a matrix $\boldsymbol{\Phi}\in\mathbb{C}^{N\times r}$ comprising $r$ orthogonal columns of length $N$,

[TABLE]

The orthonormal vectors $\{\boldsymbol{\phi}_{1},\cdots,\boldsymbol{\phi}_{r}\}$ form a POD basis. Finally, a ROM (low-fidelity solution) is constructed by the Galerkin projection of $\mathbf{u}$ onto the low-dimensional space spanned by the POD basis,

[TABLE]

Substituting (2.9) into (2.4) and projecting onto the low-dimensional space, yields an equations for the vector of coefficients $\hat{\mathbb{u}}^{n+1}$:

[TABLE]

To deal with the nonlinearity of (2.10) numerically, one might use Newton iteration or other efficient methods [8, 9].

2.2 DMD

We start by recasting the spatially discretized (2.1) in the form of a general nonlinear dynamic system

[TABLE]

where $\mathbb{u}\in\mathcal{M}\subset\mathbb{R}^{N}$, $\mathcal{M}$ is a smooth $N$-dimensional manifold, and $\mathcal{N}$ is a finite dimensional nonlinear operator. Given a flow map $\mathcal{N}_{t}:\mathcal{M}\to\mathcal{M}$,

[TABLE]

the time-discretized counterpart of (2.11) is

[TABLE]

The DMD method approximates the modes of the so-called Koopman operartor:

Definition 2.1** (Koopman operator [10]).**

The Koopman operator $\mathcal{K}$ for the nonlinear dynamic system (2.11) is an infinite-dimensional linear operator that acts on all observable functions $g:\mathcal{M}\to\mathbb{R}$ such that

[TABLE]

The discrete-time Koopman operator $\mathcal{K}_{t}$ for the discrete dynamic system (2.13) is defined by

[TABLE]

In practice, the most accessible observable is usually the state itself. Thus conventional DMD method generates a ROM by seeking a truncated finite approximation of the Koopman operator $\mathcal{K}$ coorperating with the chosen observable function $g$ as an identity map. This is done by splitting the data matrix $\mathbb{X}$ into two,

[TABLE]

and using these two datasets to approximate the eigenvalues and eigenvectors of $\mathcal{K}$ by means of the following algorithm [10].

Algorithm 2.1**.**

DMD algorithm

Compute SVD of matrix $\mathbb{X}_{1}\approx\mathbb{U}\boldsymbol{\Sigma}\mathbb{V}^{*}$ with $\mathbb{U}\in\mathbb{C}^{N\times r}$, $\boldsymbol{\Sigma}\in\mathbb{R}^{r\times r}$ and $\mathbb{V}\in\mathbb{C}^{r\times m}$, where $r$ is the truncated rank chosen by a certain criteria, e.g. (2.7). 2. 2.

Compute $\tilde{\mathbb{K}}=\mathbb{U}^{*}\mathbb{X}_{2}\mathbb{V}\boldsymbol{\Sigma}^{-1}$ as an $r\times r$ low-rank approximation of $\mathcal{K}$. 3. 3.

Compute eigendecomposition of $\tilde{\mathbb{K}}$: $\tilde{\mathbb{K}}\mathbb{W}=\mathbb{W}\boldsymbol{\Lambda}$, $\boldsymbol{\Lambda}=(\lambda_{k})$. 4. 4.

Reconstruct eigendecomposition of $\mathcal{K}$, whose eigenvalues and eigenvectors are $\boldsymbol{\Lambda}$ and $\boldsymbol{\Phi}=\mathbb{U}\mathbb{W}$, respectively.

Each column of $\boldsymbol{\Phi}$ is a DMD mode corresponding to a particular eigenvalue in $\boldsymbol{\Lambda}$. With the approximated eigenvalues and eigenvectors of $\mathcal{K}$ in hand, the projected future solution can be constructed analytically for all times in the future. In particular, at each future time $t=t^{n+1}$,

[TABLE]

where $\mathbb{b}=\Phi^{-1}\mathbb{u}^{0}$ is an $r\times 1$ vector representing the initial amplitude of each mode. Note that no more iterations are needed in the prediction. The solution at any future time can be approximated directly with (2.17) using only information from the first $m$ snapshots of the HFM.

2.3 Challenge Posed by Translational Problems

Both POD and DMD have been used to construct LFMs for a wide range of problems with high accuracy. However, ROMs constructed with such SVD-based methods are known to have poor performance for translational problems, such as an advection-dominated version of (2.1). To illustrate this phenomenon, we consider a linear advection-diffusion equation, i.e., (2.1) with constant $f$ and $D$, defined on $(x,t)=[0,2]\times[0,1]$. This equation is subject to the initial condition $u(x,0)=0.5\exp[-(x-0.3)^{2}/0.05^{2}]$ and the boundary conditions $u(0,t)=u(2,t)=0$. The space domain $[0,2]$ is discretized into $N=2000$ intervals and time domain $[0,1]$ is discretized into $M=1000$ steps. Both DMD and POD algorithms use the same dataset consisting of $m=250$ snapshots.

To achieve a diffusion-dominated regime, we set $f=10^{-4}$ and $D=10^{-2}$ in some consistent units. Figure 1 provides a visual comparison of the reference solution with its counterparts obtained with DMD and POD, both with $r=20$ SVD rank truncation. Although not shown here, and consistent with the earlier findings reported in [20], the DMD- and POD-based LFMs are of spectral accuracy in a relatively small subspace of time ($t<0.3$), with the relative error increasing with time. POD has slightly better accuracy than DMD due to the iterations in the subspace, but DMD is considerably faster because of its iteration-free nature.

An advection-dominated regime is achieved by setting $f=1.0$ and $D=10^{-3}$. Figure 2 reveals that both DMD and POD fail to capture the system dynamics, yielding unphysical (oscillatory and negative) predictions. This failure cannot be remedied by increasing the SVD rank truncation $r$: increasing $r$ from 20 to 30 does not improve the prediction’s accuracy, either quantitatively or qualitatively. These results highlight the main challenge translational problems pose for the SVD-based methods. Given the first 250 snapshots of the high-fidelity solution, SVD extracts dominant DMD/ POD modes $\phi_{i}$ from the region the wave has encountered; in our example, the subdomain $[0,1]$. As time increases, the wave solution encounters other parts of the computational domain; in our example, at later times the dominant signal lies mostly in the subdomain $[1,2]$. Specifictly, one can observe that the dominant DMD/POD modes have fluctuations only in the subdomain $[0,1]$ and stay flat [math] in the subdomain $[1,2]$ in Figure 3. It is therefore not surprising that a ROM constructed from dominant modes in $[0,1]$ does not serve as an accurate surrogate for the rest of the computational domain.

3 Lagrangian Reduced-Order Models

Motivated by construction of a POD-based ROM for the advection-diffusion equation (2.1) within the Lagrangian framework [30], we propose a Lagrangian DMD. In the semi-Lagrangian frame, (2.1) is written as

[TABLE]

Applying a first-order discretization to this system gives

[TABLE]

where $\mathcal{P}_{n}$ stands for the interpolation in the grid $\mathbb{x}^{n}=(x_{1}^{n},\cdots,x_{N}^{n})^{\top}$ and $\mathbb{x}^{0}$ is the starting uniform grid.

Or, in vector form,

[TABLE]

Here $\mathbb{x}^{n}=[x_{1}^{n},\cdots,x_{N}^{n}]^{\top}$ denotes the locations of the Lagrangian computational grid at the $n$th time step, $\mathbb{\tilde{u}}^{n}$ is the interpolation from the Lagrangian grid to the Eulerian grid, and $\mathcal{D}_{2}$ represents the discrete approximation of the second derivative on the uniform Eulerian grid at the $n$th time step.

3.1 POD

We arrange $m$ snapshots of the HFM in the Lagrangian framework into a data matrix of size $2N\times m$,

[TABLE]

Applying the conventional POD of section 2.1 to the data matrix in (3.4), one obtains a POD basis $\boldsymbol{\Phi}$ analogous to 2.8 for the space-solution vector $[\mathbf{x};\mathbf{u}]^{\top}$. Then the Lagrangian solution is approximated by

[TABLE]

Inserting (3.5) into (3.3) and projecting onto the subspace spanned by $\boldsymbol{\Phi}$, one obtains the solution vector $[\hat{\mathbb{x}}^{n+1};\hat{\mathbb{u}}^{n+1}]^{\top}$ by solving the following equation:

[TABLE]

Several complications can arise when applying Lagrangian POD in practice. If only Eulerian HFM data are available, i.e., in the absence of the grid deformation $\mathbb{x}^{n}$ computed with an Eulerian HFM, one can construct an optimal Lagrangian basis by following the strategy proposed in [30, Sec. 3.3]. Another potential complication is a Lagrangian grid entanglement. There is no guarantee that an approximation of the Lagrangian moving grid in the low-dimensional subpspace preserves the topological properties of the original HFM simulation. In many cases, e.g., when characteristic lines intersect each other, the Lagrangian grids in the projected space are severely distorted, inducing numerical instabilities. One strategy for ameliorating this problem is to solve the diffusion step back to stationary Eulerian grid by interpolation between the Eulerian and Lagrangian grids [30, Sec. 3.4]. This procedure can reduce the method’s efficiency and accuracy.

3.2 DMD

The fundamental concept behind the Koopman operator theory is to transform the finite-dimensional nonlinear problem (2.13) in the state space into the infinite-dimensional linear problem (2.15) in the observable space. Compared to Eulerian framework, Lagrangian framework contains more informative physical quantanties as candidates of the obervables, making Lagrangian DMD to fit intuitively and naturally into the Koopman operator theory. We briefly review the ethics in choosing observable functions below for the sake of completeness, followed by a description of our approach for approximating the underlying Koopman operator with Lagrangian DMD.

Since $\mathcal{K}_{t}$ in Definition 2.1 is an infinite-dimensional linear operator, it is equipped with infinitely many eigenvalues $\{\lambda_{k}\}_{k=1}^{\infty}$ and eigenfunctions $\{\phi_{k}\}_{k=1}^{\infty}$. In practice, one deals with a finite number of the eigenvalues and eigenfunctions. The following assumption underpins the finite approximation and is essential to the choice of observables.

Assumption 3.1**.**

Let $\mathbf{y}$ denote a vector of observables,

[TABLE]

If the chosen observable $\mathbb{g}$ is restricted to an invariant subspace spanned by eigenfunctions of the Koopman operator $\mathcal{K}_{t}$, then it induces a linear operator $\mathbb{K}$ that is finite-dimensional and advances these eigen-observable functions on this subspace [31].

Based on the above assumption, DMD algorithm [10] is applied to approximate the eigenvalues and eigenfunctions of $\mathbb{K}$ using the collected snapshot data in the observable space:

Algorithm 3.1**.**

Physics-aware DMD algorithm:

Create data matrices of observables $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$ as

[TABLE]

where each column is given by $\mathbb{y}^{k}=\mathbb{g}(\mathbb{u}^{k})$.

• 1.

Compute SVD of the matrix $\mathbb{Y}_{1}\approx\mathbb{U}\boldsymbol{\Sigma}\mathbb{V}^{*}$ with $\mathbb{U}\in\mathbb{C}^{p\times r},\boldsymbol{\Sigma}\in\mathbb{R}^{r\times r},\mathbb{V}\in\mathbb{C}^{r\times m}$, where $r$ is the truncated rank chosen by certain criteria.

• 2.

Compute $\tilde{\mathbb{K}}=\mathbb{U}^{*}\mathbb{X}^{\prime}\mathbb{V}\boldsymbol{\Sigma}^{-1}$ as an $r\times r$ low-rank approximation for $\mathbb{K}$.

• 3.

Compute eigendecomposition of $\tilde{\mathbb{K}}$: $\tilde{\mathbb{K}}\mathbb{W}=\mathbb{W}\boldsymbol{\Lambda}$, $\boldsymbol{\Lambda}=(\lambda_{k})$.

• 4.

Reconstruct eigendecomposition of $\mathbb{K}$. Eigenvalues are $\boldsymbol{\Lambda}$ and eigenvectors are $\boldsymbol{\Phi}=\mathbb{U}\mathbb{W}$.

• 5.

Future $\mathbb{y}_{\text{DMD}}^{n+1}$ can be predicted by

[TABLE]

with $\mathbb{b}=\boldsymbol{\Phi}^{-1}\mathbb{y}_{1}$.

• 6.

Transform from observables back to state-space:

[TABLE]

In data-driven modeling, judicious selection of the observables is crucial to the accuracy and efficiency of a Koopman operator’s approximation. Identification of general rules for choosing the observables continues to be a subject of ongoing research. For example, the use of measurements of the state variable $u(x,t)$ as an observable led to the poor DMD performance in the advection-dominated regime (Figure 2). A Lagrangian formulation of the problem provides a means of identification of optimal observables. Indeed, the physics of advection-dominated systems suggests that the location of a moving particle is a key quantity, which is as important as the value of the state variable at that location. It is therefore natural to introduce an observable function that keeps track of both. Thus we choose our observable to be

[TABLE]

Then, we follow Algorithm 3.1.

4 Numerical Experiments

To ascertain the accuracy and robustness of the Lagrangian DMD, we use it to construct ROMs for a series of linear and nonlinear advection-dominated problems. In all tests, the reference solutions are computed in the Eulerian framework using (2.3). The space domain $[0,2]$ is discretized into $N=2000$ intervals and the time domain $[0,1]$ is discretized into $M=1000$ steps. Both Lagrangian DMD and Lagrangian POD algorithms use $m=250$ snapshots (up to $t=0.25$) as a training dataset. The rank truncation criteria (2.7) with $\varepsilon=10^{-8}$ is used. The error bound derived in [20] is reviewed in the following theorem and plotted in each example as an estimate of the observable.

Theorem 4.1**.**

Define the global truncation error

[TABLE]

Then, for $n\geq m$,

[TABLE]

where $\varepsilon_{m}$ is a constant only depending on the number of snapshots $m$.

4.1 Linear Advection Equation

We start by considering (2.1) with $f\equiv 1$ and $D\equiv 0$. The resulting linear advection equation is defined on $(x,t)\in(0,2)\times(0,1]$, and is subject to the initial condition

[TABLE]

Figure 4 provides a visual comparison between the reference solution, obtained with the numerical scheme (2.3), and solutions of the ROMs constructed with either Lagrangian DMD or Lagrangian POD. Unlike their conventional (Eulerian) counterparts (see Figure 2), both Lagrangian DMD and Lagrangian POD capture the solution dynamics in the extrapolating mode, i.e., for $t>0.25$.

A more quantitative comparison of the relative performance of the two SVD-based strategies is presented in Figure 5 in terms of the global truncation error defined in (4.1). Both Lagrangian DMD and Lagrangian POD capture the advection with high accuracy. Due to the linearity and conservation property of this problem, the ROMs constructed by the two methods are of machine error. Thus, the error bound developed in [20] is not tight but sufficient to serve as an indicator of successful approximation.

4.2 Linear Advection-Diffusion Equation

Next, we consider (2.1) with $f\equiv 1$ and $D\equiv 0.01$. The resulting linear advection-diffusion equation is defined on $(x,t)\in(0,2)\times(0,1]$, and is subject to the initial and boundary conditions (4.3). The choice of the parameter values $f$ and $D$ ensures that the system is in the advection-dominated regime, for which the conventional POD and DMD fail.

Figure 6 provides a visual comparison between the reference solution $u(x,t)$ and those obtained with the ROMs. The latter capture the system’s dynamics, although their estimates of the solution tails become less accurate with time. This suggests that Lagrangian DMD and POD are capable of identifying the low-rank structure of the advection-diffusion dynamics in the advection-dominated regime.

Figure 7 indicates that the Lagrangian DMD and POD have a near identical accuracy, which deteriorates with extrapolation time $t>0.25$. The error bound for the DMD estimate of the observable $g(u)$ is appreciably tighter than in the case of advection (Figure 5). With the error bounds, one can design an algorithm combining short-term computation of HFM with long-term computation of LFM.

4.3 Inviscid Burgers Equation

The inviscid Burgers equation is recovered from (2.1) by setting $f\equiv u$ and $D\equiv 0$. We define this equation on $(x,t)\in(0,2\pi)\times(0,1]$, subject to the initial conditions

[TABLE]

Figure 8 provides a graphical illustration of the Lagrangian ROMs’ ability to capture the system state dynamics in this nonlinear hyperbolic problem. Since the Lagrangian description treats first-order hyperbolic conservation laws, such as the inviscid Burgers equation, exactly via the method of characteristics, the addition of the particle trajectories $x(t)$ to the set of observables ensures that the Lagrangian POD and DMD are both of machine error accuracy (Figure 9). Again the error bound serves as an indicator of accurate ROMs.

The level-set method provides an alternative way to interpret the first-order hyperbolic conservation laws. In the Appendix, we report our experiments with the level-set DMD, which is essentially a Lagrangian DMD for two-dimensional linear advection equation.

4.4 Viscous Burgers Equation

The viscous Burgers equation is obtained from (2.1) be setting $f\equiv u$ and $D\equiv 0.1$. Again, this equation is defined on $(x,t)\in(0,2\pi)\times(0,1]$ and is subject to the initial and boundary conditions (4.4).

For this nonlinear problem, Lagrangian DMD is visually more accurate that Lagrangian POD (Figure 10), especially at later times. This is confirmed by plotting the error $\mathcal{E}$ in Figure 11. As mentioned in [32], the Lagrangian grid might become distorted (especially in the presence of large gradients $\partial_{x}u$) during the compressing process of ROM in the POD algorithm. The error estimate of the observable does a good job in evaluating the bound.

4.5 Computational costs

Table 1 collates the rank of the ROMs and the computational times of the HFM and the Lagrangian DMD- and POD-based ROMs. In some cases, the SVD dominates the computational time of the ROM. Once the basis is constructed, the computation in the low rank subspace is much faster. This explains why the POD-based ROM of Test 3 takes more time to compute than the HFM. In other cases, the ROMs are much more efficient than the HFM computations. DMD is the most efficient methods because of its iteration-free nature.

5 Conclusions

In this paper, we investigate the issue of translational problem for conventional proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD) in the Eulerian framework. A new physic-aware DMD, based on the Lagrangian framework, is proposed to overcome the shortcomings of reduced order models (ROMs) of advection-dominated nonlinear phenomena. Characteristic lines, an important physical quantity in such systems, are taken into account in order to learn the Koopman operator of the underlying dynamics. The Lagrangian framework provides an optimal choice of observable functions for learning the Koopman operator. It allows one to construct a ROM in a relatively small subspace by using the DMD algorithm with satisfactory accuracy. Compared to the Lagrangian POD method, physics-aware DMD is more efficient computationally thanks to its iteration-free nature.

One possible direction for future work is to investigate the advection-diffusion system in Lagrangian coordinates [33]. Interpolation between Eulerian grid and Lagrangian grid will not be needed anymore but careful discretization of the diffusion operator will need to be handled. Existing numerical studies in Lagrangian coordinates and related methods [34, 35] could be explored as guidelines of choosing physical observables in reduced order modeling.

All the numerical tests presented in this paper are shock-free. Once shock is formed, the Lagrangian formulation (3.1) becomes invalid. Although one can still make the scheme (3.2) work by numerical remedies, instability or unphysical solutions could appear when sharp gradients or shocks occur. The instability could become more severe in the compressed low-dimensional space [32]. The modifications in [32] bypass this issue by compensating computational costs in projecting back to the Eulerian grid. From the perspective of physic-aware data-driven modeling, we realize that significant information like shock formation time, shock location and shock speed is not interpreted well enough from data. In another word, other quantities should be chosen as essential observables in order to learn the underlying Koopman operator. We leave this line of research for a follow-up study.

Acknowledgements

This work was supported in part by Defense Advanced Research Project Agency under award number 101513612 and by Air Force Office of Scientific Research under award number FA9550-17-1-0417.

Appendix A Level-set DMD for Hyperbolic Conservation Laws

The level-set approach [36] provides another way to interpret conservation laws. Supposed that a state variable $u(x,t)$ satisfies the one-dimensional conservation law

[TABLE]

with $f(u)\geq 0$. Its corresponding level-set formulation,

[TABLE]

is a linear two-dimensional transport equation for the dependent variable $c(x,y,t):\mathbb{R}^{2}\times\mathbb{R}^{+}\to\mathbb{R}$. Together with a Lipschitz-continuous initial function $c_{0}$, which embeds the initial data $u_{0}$ (see the example below), the zeroth-level set of $c(x,y,t)$, i.e., the $x-y$ contour of the solution to $c(x,y,t)=0$, gives the solution to the conservation law (A.1), $u(x,t)$.

By way of example, we consider the inviscid Burgers equation (4.4). Its level-set formulation is

[TABLE]

We apply Lagrangian DMD to construct a ROM for this two-dimensional linear advection equation from $m=250$ snapshots. Figure 12 demonstrates that the ROM with the $r=3$ rank truncation approximates the HFM solution $u(x,t)$ with satisfactory accuracy. Although solving a two-dimensional linear problem takes more computational time than solving the nonlinear one-dimensional problem in this case, the level-set DMD provides another venue for investigation of physics-aware DMD that might have efficient applications to other problems.

Reference

• [1]

D. M. Tartakovsky, M. Dentz, Diffusion in porous media: Phenomena and mechanisms, Transp. Porous Media (2019).

doi:10.1007/s11242-019-01262-6.

• [2]

P. Benner, S. Gugercin, K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Rev. 57 (4) (2015) 483–531.

• [3]

A. Quarteroni, G. Rozza, et al., Reduced order methods for modeling and computational reduction, Vol. 9, Springer, 2014.

• [4]

S. Acharjee, N. Zabaras, A concurrent model reduction approach on spatial and random domains for the solution of stochastic PDEs, Int. J. Num. Meth. Engrg. 66 (12) (2006) 1934–1954.

• [5]

P. Holmes, J. L. Lumley, G. Berkooz, C. W. Rowley, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge Univ. Press, 2012.

• [6]

J. L. Lumley, Stochastic tools in turbulence, Courier Corporation, 2007.

• [7]

S. Volkwein, Model reduction using proper orthogonal decomposition, Lecture Notes, Institute of Mathematics and Scientific Computing, University of Graz. see http://www. uni-graz. at/imawww/volkwein/POD. pdf 1025 (2011).

• [8]

M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera, An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations, Comptes Rendus Mathematique 339 (9) (2004) 667–672.

• [9]

S. Chaturantabut, D. C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM J. Sci. Comput. 32 (5) (2010) 2737–2764.

• [10]

J. N. Kutz, S. L. Brunton, B. W. Brunton, J. L. Proctor, Dynamic mode decomposition: data-driven modeling of complex systems, Vol. 149, SIAM, 2016.

• [11]

I. Mezić, Analysis of fluid flows via spectral properties of the Koopman operator, Annu. Rev. Fluid Mech. 45 (2013) 357–378.

• [12]

B. O. Koopman, Hamiltonian systems and transformation in Hilbert space, Proc. Natl. Acad. Sci. U.S.A. 17 (5) (1931) 315.

• [13]

I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlin. Dyn. 41 (1-3) (2005) 309–325.

• [14]

I. Mezić, A. Banaszuk, Comparison of systems with complex behavior, Physica D 197 (1-2) (2004) 101–133.

• [15]

C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, D. S. Henningson, Spectral analysis of nonlinear flows, J. Fluid Mech. 641 (2009) 115–127.

• [16]

S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proc. Natl. Acad. Sci. U.S.A. 113 (15) (2016) 3932–3937.

• [17]

M. Schmidt, H. Lipson, Distilling free-form natural laws from experimental data, Science 324 (5923) (2009) 81–85.

• [18]

D. Duke, J. Soria, D. Honnery, An error analysis of the dynamic mode decomposition, Exper. Fluids 52 (2) (2012) 529–542.

• [19]

M. Korda, I. Mezić, On convergence of extended dynamic mode decomposition to the Koopman operator, J. Nonlin. Sci. 28 (2) (2018) 687–710.

• [20]

H. Lu, D. M. Tartakovsky, Prediction accuracy analysis for dynamic mode decomposition, ArXiv (2019).

• [21]

D. Amsallem, M. J. Zahr, C. Farhat, Nonlinear model order reduction based on local reduced-order bases, Int. J. Num. Meth. Engrg. 92 (10) (2012) 891–916.

• [22]

D. J. Lucia, Reduced order modeling for high speed flows with moving shocks, Tech. rep., Air Force Inst. of Tech., Wright-Patterson Air Force Base, OH (2001).

• [23]

K. Carlberg, Adaptive h-refinement for reduced-order models, Int. J. Num. Meth. Engrg. 102 (5) (2015) 1192–1210.

• [24]

J. N. Kutz, X. Fu, S. L. Brunton, Multiresolution dynamic mode decomposition, SIAM J. Appl. Dyn. Syst. 15 (2) (2016) 713–735.

• [25]

J.-F. Gerbeau, D. Lombardi, Approximated Lax pairs for the reduced order integration of nonlinear evolution equations, J. Comput. Phys. 265 (2014) 246–269.

• [26]

M. E. Kavousanakis, R. Erban, A. G. Boudouvis, C. W. Gear, I. G. Kevrekidis, Projective and coarse projective integration for problems with continuous symmetries, J. Comput. Phys. 225 (1) (2007) 382–407.

• [27]

M.-L. Rapún, J. M. Vega, Reduced order models based on local POD plus Galerkin projection, J. Comput. Phys. 229 (8) (2010) 3046–3063.

• [28]

C. W. Rowley, I. G. Kevrekidis, J. E. Marsden, K. Lust, Reduction and reconstruction for self-similar dynamical systems, Nonlinearity 16 (4) (2003) 1257.

• [29]

C. W. Rowley, J. E. Marsden, Reconstruction equations and the Karhunen–Loève expansion for systems with symmetry, Physica D 142 (1-2) (2000) 1–19.

• [30]

R. Mojgani, M. Balajewicz, Lagrangian basis method for dimensionality reduction of convection dominated nonlinear flows, arXiv preprint arXiv:1701.04343 (2017).

• [31]

S. L. Brunton, B. W. Brunton, J. L. Proctor, J. N. Kutz, Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control, PloS One 11 (2) (2016) e0150171.

• [32]

J. Morton, A. Jameson, M. J. Kochenderfer, F. Witherden, Deep dynamical modeling and control of unsteady fluid flows, in: Advances in Neural Information Processing Systems, 2018, pp. 9278–9288.

• [33]

J.-L. Thiffeault, Advection-diffusion in Lagrangian coordinates, Phys. Lett. A 309 (5-6) (2003) 415–422.

• [34]

M. Shashkov, B. Wendroff, A composite scheme for gas dynamics in Lagrangian coordinates, J. Comput. Phys. 150 (2) (1999) 502–517.

• [35]

S. E. Hieber, P. Koumoutsakos, A Lagrangian particle level set method, Journal of Computational Physics 210 (1) (2005) 342–367.

• [36]

Y.-H. Tsai, Y. Giga, S. Osher, A level set approach for computing discontinuous solutions of Hamilton-Jacobi equations, Math. Comput. 72 (241) (2003) 159–181.