Dynamic Mode Decomposition for Construction of Reduced-Order Models of Hyperbolic Problems with Shocks
Hannah Lu Daniel M. Tartakovsky
math.NA cs.NA


 
               Figures/
              





mymainaddress]Hannah Luemail1
[email1]email: 
mymainaddress]Daniel M. Tartakovskymycorrespondingauthor
[mycorrespondingauthor]Corresponding author
tartakovsky@stanford.edu
[mymainaddress]Department of Energy Resources Engineering, Stanford University, Stanford, CA 94305, USA


Construction of reduced-order models (ROMs) for hyperbolic conservation laws is notoriously challenging mainly due to the translational property and  nonlinearity of the governing equations. While the Lagrangian framework for ROM construction resolves the translational issue, it is valid only before a shock forms. Once that occurs, characteristic lines cross each other and projection from a high-fidelity model space onto a ROM space distorts a moving grid, resulting in numerical instabilities. We address this grid distortion issue by developing a  physics-aware dynamic mode decomposition (DMD) method based on hodograph transformation. The latter provides a map between the original nonlinear system and its linear counterpart, which coincides with the Koopman operator. This strategy is consistent with the spirit of physics-aware DMDs in that it retains information about shock dynamics. Several numerical examples are presented to validate the proposed physics-aware DMD approach to constructing accurate ROMs.



Conservation law Koopman operator Proper orthogonal decomposition Hodograph transformation









§ INTRODUCTION


Since introduction of Euler equations, hyperbolic conservation laws play a significant role in gas dynamics, astrophysics, plasma, traffic flow, multiphase flow in porous media <cit.> and other fields of science and engineering. Wave-like solutions of hyperbolic equations can exhibit various rarefaction and shock behaviors, whose occurrence  strongly depends on a functional form of the flux function. Discontinuity and uniqueness of such solutions pose challenges in theoretical treatment of hyperbolic conservation laws <cit.>. Theoretical advances, such as entropy conditions and the concept of a weak solution <cit.>, ameliorate this difficulty by providing physical interpretation to these solutions. Likewise, numerical  high-resolution methods have been designed to resolve nonlinearities and accurately capture shocks <cit.>.  Although continued developments in scientific computing have improved the performance of high-resolution simulations, their computational cost is often too high to model complex systems at spatiotemporal resolutions and scales of interest. The cost can become prohibitive when used in the context of uncertainty quantification or data assimilation, both of which require a large number of repeated forward model runs.

Reduced-order models (ROMs) provide an efficient alternative to their high-fidelity, physics-based counterparts that can be deployed in large-scale multiphysics simulations. Robust tools for construction of ROMs for problems described by ordinary differential equations or parabolic partial differential equations (PDEs) include proper orthogonal decomposition (POD) <cit.> and dynamic mode decomposition (DMD) <cit.>. 

The challenge of extending these techniques to hyperbolic or advection-dominated parabolic PDEs with smooth solutions was met in <cit.> through development of the physics-aware DMD and POD approaches within a Lagrangian framework.  However, in the presence of strong shocks and/or sharp gradients, Lagrangian POD methods can generate numerical instability caused by grid distortion <cit.>. Once characteristics of a nonlinear hyperbolic PDE intersect each other, the projection from a high-dimensional manifold of the high-fidelity model (HFM) onto a low-dimensional subspace of the low-fidelity model (LFM), e.g., ROM, is not guaranteed and typically fails to preserve topological properties of the original HFM. We elaborate on this point in section <ref>, in terms relevant to DMD.

We use hodograph transformation <cit.> to resolve this outstanding issue in construction of ROMs for PDEs with discontinuous solutions and shocks. Hodograph diagrams have originated in meteorology to plot wind from soundings of the Earth's atmosphere. Since then,  hodograph transformation morphed into a technique designed to transform nonlinear PDEs into linear ones by  interchanging the dependent and independent variables in a PDE. Hodograph-type transformations have been used to find quasilinear analogues of  semi-linear equations, and to derive new analytical solutions to special classes of PDEs <cit.>. 

Advantages of mapping nonlinear PDEs onto their linear counterparts are self-evident: analytical tools available for linear PDEs provide better understanding of a solution's behavior, and numerical solvers for linear systems are both easier to implement and computationally cheaper. 

The Koopman operator theory <cit.> shares the goal of hodograph transformation: a Koopman operator is an infinite-dimensional linear operator that represents the underlying finite-dimensional nonlinear dynamic system by judiciously choosing observable functions. It is also similar in its goal to integral transformations that map certain classes of nonlinear PDEs onto their linear counterparts. For example, the Cole-Hopf transformation and the Kirchhoff transformation map, respectively, Burgers' equation and a class of nonlinear diffusion (heat conduction) equations onto a linear diffusion equation. 

These integral transformations have been used in the context of the Koopman operator theory and DMD/POD to constructed ROMs for Burgers' equation <cit.> and a nonlinear diffusion equation <cit.>. A major goal of our study is to establish clear connections between hodograph transformation and the Koopman theory. This relationship between the two is then used both to identify observables for a Koopman operator via hodograph transformation and to construct ROMs for hyperbolic conservation laws with shocks via DMD.

A general procedure of the physics-aware DMD algorithm and its connection with Koopman operators are reviewed in section <ref>. This section also contains a numerical demonstration of the failure of Lagrangian-based ROMs to capture the dynamics described by conservation laws with shocks. We illustrate the use of hodograph transformation by analyzing the inviscid Burgers' equation (<ref>) and more general hyperbolic PDEs with a convex flux function (appendix <ref>). In section <ref>, we combine hodograph transformation with the Koopman operator theory to design a physics-aware DMD algorithm for construction of ROMs for conservation laws with shocks. Several numerical tests are presented in section <ref> to validate the proposed physics-aware DMD approach. Main conclusions drawn from our study are summarized in section <ref>.




§ CONSTRUCTION OF ROMS AND THEIR FAILURE FOR PROBLEMS WITH SHOCKS



Consider a state variable u(x,t): [a,b] × [0,T] →ℝ, where the constants a, b, T ∈ℝ. The dynamics of u(x,t) is described by a one-dimensional scalar conservation law 

    ∂ u/∂ t+∂ F(x,t,u)/∂ x =0
      or  ∂ u/∂ t+f(x,t,u)∂ u/∂ x =0,    f(x,t,u) = ∂ F(x,t,u)/∂ u.

This hyperbolic PDE is subject to the initial condition u(x,t=0) = u_0(x) and, when appropriate (i.e., when |a|, |b| < ∞), boundary conditions at a and/or b. 











The intervals [0,T] and [a,b] are discretized with (N+1) and J nodes separated, respectively, by Δ t and Δ x. To be specific, we solve (<ref>) with a conservative first-order upwind scheme <cit.> 

    u_j^n+1 = u_j^n-Δ t/Δ x(F_j+1/2^n-F_j-1/2^n),

where n = 0,…,N indicates the nth time step, with n=0 corresponding to t=0, and n=N to t = T;  j = 1,…,J denotes the jth spatial node, such that j=1 and J coincide with x = a and b, respectively; and

    F_j+1/2^n =  F(·, u_j+1^n)+F(·, u_j^n) /2 - |a_j+1/2^n| u_j+1^n-u_j^n/2,
       a_j+1/2^n = {   F_j+1^n-F_j^n/u_j+1^n-u_j^n         u_j+1^n≠ u_j^n,
       f(·, u_j)         u_j+1^n=u_j^n.
    .

A numerical solution provided by (<ref>) with sufficiently small Δ t and Δ x are referred to as a reference HFM throughout the paper.

Standard (Eulerian) approaches to construction of a ROM for (<ref>) often fail due to the traveling-wave nature of its solution <cit.>. In a shock-free scenario, the Lagrangian framework can resolve the translational issue in the POD or DMD approaches to ROMs by keeping track of the characteristic lines. 

In the Lagrangian framework, (<ref>) becomes

    {   dx/dt = f(x,t,u),      x(0) = η
       du/dt =0,      u(η,0) = u_0(η),
    .

where η∈ℝ is a label of the characteristic x(t). As in the Eulerian case, we use the uniform discretization of the time interval [0,T], such that 0=t^0<t^1<…<t^N=T with time step Δ t = t^n+1-t^n.  At time t=0, the space, [a,b], is discretized with a uniform mesh 𝐱^0 = [x^0_1,…, x^0_J]^⊤ of mesh size Δ x^0  = x^0_j+1-x^0_j. Unlike its Eulerian counterpart, the spatiotemporal discretization of u(x,t) in the Lagrangian framework, 𝐮^n = [u_1^n,…, u_J^n]^⊤ for n = 0,…,N, may be nonuniform in space due to the temporal evolution of the grid nodes x_j(t). The backward Euler discretization, used in <cit.>, transforms (<ref>) into


    {   x_j^n+1 = x_j^n+Δ t f(x_j^n+1, (n+1) Δ t, u_j^n+1),
       u_j^n+1 = u_j^n
    .


or, in vector form,

    {   𝐑_x(𝐱^n+1) ≡𝐱^n+1-𝐱^n-Δ t𝐟^n+1(·, 𝐮 ^n+1) = 0, 
       𝐑_u(𝐮^n+1) ≡𝐮^n+1-𝐮^n = 0,
    .

where 𝐱^n = [x_1^n, …, x_J^n]^⊤ denotes the nodes of the Lagrangian computational grid at the nth time step. This numerical scheme involves N iterations in the two high-dimensional J × 1 vectors, 𝐱^n+1 and 𝐮^n+1. It provides a Lagrangian HFM. 

To construct a ROM, a data set consisting of a sequence of M solution snapshots (M ≤ N and, ideally, M ≪ N) is collected from the HFM. Since 𝐮^n is conservative and invariant in time, we only need the data matrix 𝐗 with M snapshots of the Lagrangian grid 𝐱^n:


    𝐗 = [   |   |   |; 𝐱^1 𝐱^2   ⋯ 𝐱^M;   |   |   | ].


In the next two subsections, we briefly revisit the algorithms of Lagrangian POD <cit.> and Lagrangian DMD <cit.> used to construct a ROM.




 §.§ Lagrangian POD


Identification of the POD modes is based on a reduced singular value decomposition (SVD),

    𝐗 = 𝐔Σ𝐕^*,

where 𝐔∈ℂ ^J × K,Σ = ℂ ^K× K,𝐕∈ℂ^M × K, K is the rank of the matrix 𝐗 approximated by the reduced SVD. Further rank truncation can be achieved by using the energy criterion,


    r = min_k{σ_k/∑_k=1^K σ_k<ε},


where σ_k are the diagonal elements of Σ, and ε is a small number (tolerance), chosen to be ε = 10^-4 in all our numerical examples. After the truncation, one gets the POD modes

    Φ = 𝐔(:,1:r) =  [   |   |   |; ϕ_1 ϕ_2   ⋯ ϕ_r;   |   |   | ].
 
Notice that r≪ K≤min{J,M}, and the basis {ϕ_1,⋯,ϕ_r} is orthonormal. Galerkin projection in the low-dimensional space spanned by the POD basis provides a ROM (low-fidelity solution),

    𝐱_POD^n+1 = ∑_k = 1^rx̂_k^n+1ϕ_k = Φ𝐱̂^n+1.

The r× 1 vector 𝐱̂^n+1 of coefficients is computed as a solution of 

    Φ^⊤𝐑(
    Φ[      |; 𝐱̂^n+1;      | ])=0.

that is obtained by substituting (<ref>) into the first equation in (<ref>) and projecting onto the subspace spanned by Φ.




 §.§ Lagrangian DMD



The DMD algorithm <cit.> is applied to the Lagrangian grid matrix 𝐗 in (<ref>).

Lagrangian DMD algorithm


  0. Create data matrices of (M-1) observables, 𝐗_1 and 𝐗_2,

    𝐗_1 = [     |     |     |;   𝐱^1   𝐱^2     ⋯ 𝐱^M-1;     |     |     | ],     𝐗_2 = [   |   |   |; 𝐱^2 𝐱^3   ⋯ 𝐱^M;   |   |   | ].


  1. Apply SVD of matrix 𝐗_1 ≈𝐔Σ𝐕^* with 𝐔∈ℂ^J × r, Σ∈ℂ^r× r, 𝐕∈ℂ^r× M, where r is the truncated rank chosen by a certain criterion, e.g., (<ref>).

  2. Compute 𝐊̃=𝐔^*𝐗'𝐕Σ^-1 as an r× r low-rank approximation of 𝐊.

  3. Compute eigen-decomposition of 𝐊̃: 𝐊̃𝐖 = 𝐖Λ, Λ = (λ_k).

  4. Reconstruct eigen-decomposition of 𝐊. Eigenvalues are Λ and eigenvectors are Φ  = 𝐔𝐖.

  5. Future 𝐱_DMD^n+1 is predicted by

    𝐱_DMD^n+1 =ΦΛ^n+1𝐛,   n > M

with 𝐛 =Φ^-1𝐱_1.

  6. Interpret the solution in the moving grid:

    u_DMD(x_j^n, n Δ t) = u_0(x_j^0).








 §.§ ROM Failure for Problems with Shocks: Inviscid Burgers Equation



One of the most studied examples of (<ref>) is the inviscid Burgers equation


    ∂ u/∂ t+u∂ u/∂ x = 0,      u(x,0)=u_0(x);

which we define on the space-time domain (x,t) ∈ [0,2π]× [0,1]. Depending on the boundary and initial conditions, this problem admits both smooth and discontinuous solutions u(x,t). For example, a smooth solution is obtained for the periodic boundary conditions, u(0,t) = u(2π,t), and the initial data u_0(x) = 1+sin(x). In this setting, standard (Eulerian) ROMs fail, while the ROMs based on either Lagrangian POD or Lagrangian DMD perform well in terms of both accuracy and computational efficiency <cit.>.
















A solution to (<ref>) develops shocks in finite time for, e.g., a Gaussian-type initial data, 

    u_0(x) = 1/2 + 1/2exp[- (x-0.3)^2/0.01].


In the pure Lagrangian approach (<ref>), the discretization has to account for shock formation. Once the characteristic lines cross each other, the Lagrangian mesh becomes sensitive to the choice of discretization of u[x(t),t]. 

For instance, a discretization of (<ref>) with f(·, u) =u,


    {   u_j^n+1 = u_j^n,
    
       x_j^n+1 = x_j^n +Δ t/2(u_j^n+u_j^n+1),
    .

would lead to the so-called “overshoot" that admits multi-value solutions (Figure <ref>(a)), which contradicts the entropy condition. This is a typical problem with the Lagrangian framework. It should come as no surprise that an attempt to build a ROM with the Lagrangian DMD based on the faulty discretization (<ref>) likewise results in failure (Figure <ref>(b)). The Lagrangian DMD faithfully reproduces the unphysical solution obtained with the faulty discretization scheme (<ref>). In other words, the resulting unphysical ROM is not caused by the DMD algorithm itself; the data from the full Lagrangian model (<ref>) provide inaccurate and incomplete (without shock) information from the very beginning.




 

To isolate the performance of the Lagrangian DMD, we consider a numerical scheme that is known for its ability to handle shocks: the backward semi-Lagrangian method (BSLM)


    {   u_j^n+1 = u_j^n,
       x^* = x_j^n+Δ t/2u_j^n,
       x_j^n+1 = x_j^n +Δ t/2(u(x^*,t^n)+u(x^*,t^n+1)).
    .

We employed the (explicit) mid-point rule to avoid implicit iterations. As one can see from the above figure, (<ref>) gives the right physical solution with shock features. Figure <ref>(a) reveals that this numerical scheme is indeed capable to accurately approximate the solution of the inviscid Burgers equation with shocks. However, the Lagrangian DMD algorithm using snapshots from the full solution (<ref>) suffers from instability once a shock is about to form (Figure <ref>(b)). The grid becomes severely distorted once the characteristic lines intersect each other at the interface where sharp gradients of u(x,t) occur.  At the intersect, one arrival location of x corresponds to two different departure values of u. However, the DMD modes projection from the HFM to the ROM does not keep the topological information about this multivalued mapping in the ROM process, resulting in the Lagrangian grid distortion. 




The Lagrangian POD approach suffers from similar problems <cit.>. Moreover, the POD projection on the accurate Lagrangian scheme (<ref>) would still require interpolation in the high-dimensional space. One might need techniques such as DEIM <cit.> to keep the resulting ROM's efficiency. Nevertheless, extensions of POD are beyond the scope of our study; we focus on DMD-based ROMs due to their iteration-free nature.






§ HODOGRAPH TRANSFORMATION





We start with a mathematical definition of hodograph transformation reproduced from <cit.>.

A pure hodograph transform is a transformation of the form

    τ = t,   ξ = u(x,t).





For the inviscid Burgers equation (<ref>), we first consider a scenario where only one shock is developed from the initial data u_0(x) in finite time. This necessitates the following assumption.

The function u_0(x) satisfies four conditions:


  * u_0(x) is smooth;

  * u_0(x) decreases monotonically, i.e., u_0^'(x) <0 for all x; and 
 lim_x→ +∞ u_0(x) = u_R, lim_x→ -∞u_0(x) = u_L with constants u_R < u_L;

  * u_0(x) has a unique inflection point (x^*, u^*) with u^* = u_0(x^*), meaning u_0^''(x^*) = 0;

  * u_0^''' (x^*) >0.



This assumption ensures existence of an inverse function x(t,u) : [0,T]× [u_R,u_L] → [a,b] of the monotonic function u(t,x) : [0,T] × [a,b] → [u_R,u_L]. 



It follows from Definition <ref> that the inverse function x(t,u) is a pure hodograph transform. 



 §.§ Solution Before Shock Formation



With u acting as the independent variable and x as the dependent variable, hodograph transformation x = x(t,u) maps the equation for characteristics (<ref>) of the inviscid Burgers equation (<ref>), before the shock formation time t^* (defined later), onto  


    d x/dt(t,u) = u,    x(0,u) = x_0(u);     for   (t,u) ∈ [0,t^*) × [u_R,u_L].


Assumption <ref> translates into conditions on the function x_0(u):


  * x'_0 (u)<0;

  * x_0(u) has a unique inflection point at (u^*,x^*);

  * x”'_0(u^*)<0.

Differentiation of (<ref>) with respect to u gives


    ∂^2 x/∂ t ∂ u(t,u) = 1,

from which

    ∂ x/∂ u(t,u) = x'_0(u) + t,     for   (t,u) ∈ [0,t^*) × [u_R,u_L].


Since x'_0(u)<0, 



∂_u x (t,u) < 0 as long as 

t^* = min_u[-x'_0(u)] = -x'_0(u^*).

Thus defined t^* determines the time of shock formation. The shock location is x^* = x(u^*,x^*).




 §.§ Solution After Shock Formation


At times t larger than t^* = -x'_0(u^*), i.e., once the shock forms, (<ref>) is no longer valid. In the (x,u) plane, one would use the entropy (Rankine-Hugoniot) condition to construct a weak formulation of Burgers' equation. Its analog in the (u,x) plane gives an equation for the shock speed s: 

    s = 1/2u_1^2 - u_2^2/u_1-u_2 =u_1+u_2/2,

where u_1(t) and u_2(t) are defined as the limits of u(t) from the top and bottom of the shock, respectively. They are computed as solutions of a system of ordinary differential equations <cit.>

    {   du_1/dt =1/2u_1-u_2/g(u_1)-t,
       du_2/dt =-1/2u_1-u_2/g(u_2)-t,
    .

where g(u) ≡ -x'_0(u). These ODEs are subject to initial conditions u_1(t^*) =u^* and u_2(t^*) = u^*.
Since s = d x^* / dt, an equation for the shock trajectory x^*(t) is


    d x^* /dt =u_1+u_2/2.





 §.§ Summary of Hodograph Solution


Under Assumption <ref>, the hodograph-transformed Burgers equation (<ref>) in the Lagrangian framework (<ref>) takes the form of the following ODEs for x(t,u): 

    {   t<t^*:       
       t>t^*:      {               u∈(u_R,u_2)∪(u_1,u_L)
                   u∈(u_2,u_1),
    .
    .

where t^* = -x'_0(u^*), and u_1 and u_2 are solutions of (<ref>).


One can show that u_1(t) is a monotonically increasing function and u_2(t) is a monotonically decreasing function, such that

    u_1≥ u^*,    u_2≤ u^*,    x'_0(u_1)+t≤ 0,    x'_0(u_2)+t≤ 0.

In many cases of interests, and in all our numerical experiments, either u_2 = u_R and u_1 = u_L or |u_2 - u_1| ≪Δ t (so that u_2 ≈ u_R and u_1 ≈ u_L). This allows one to focus on shock propagation, i.e., on (<ref>), without having to solve (<ref>).







Functions u_0(x) that do not satisfy Assumption <ref>, such as (<ref>), require a decomposition of the initial data into monotonic parts. Each monotonic piece of u_0(x) would have a unique inverse function x_0(u). The entropy condition implies that the increasing x_0(u), i.e. x'_0(u)>0, results in a rarefaction solution, which satisfies (<ref>). The union of  the rarefaction pieces and shock pieces would give the full solution.



The inviscid Burgers equation is an example of hyperbolic conservation laws with monotonically increasing flux functions F(·, u). Generalization to hyperbolic conservation laws with a convex flux is presented in appendix <ref>.





§ PHYSICS-AWARE DMD FOR CONSERVATION LAWS WITH SHOCKS



Previous theoretic investigations, e.g., <cit.>, demonstrated that the key to a DMD's success in capturing nonlinear dynamics is to identify the underlying Koopman operator. Several numerical studies <cit.> confirmed this finding. The Koopman operator theory ensures that a DMD algorithm utilizes all relevant physical information to learn the dynamics. We refer to this approach as physics-aware DMD to distinguish it from the conventional DMD that learns only from (simulations-generated) data.

We review the Koopman operator theory and analyze its connection with hodograph transformation in the context of hyperbolic conservation laws.  Then, we present our general framework for physics-aware DMD for problems with discontinuous solutions and shocks.

A suitable spatial discretization of (<ref>) leads to a nonlinear dynamical system 

    d𝐮/dt = 𝒩(𝐮),

where 𝐮(t) = (u_1,…,u_J)^⊤∈ℳ⊂ℝ^J is the solution vector with u_j(t) = u(x_j,t) and discretization nodes x_j (j=1,…,J); and 𝒩 is a finite-dimensional nonlinear operator.  A flow map 𝒩_t :ℳ→ℳ,


    𝒩_t [𝐮(t_0)] ≡𝐮(t_0+t) = 𝐮(t_0)+∫_t_0^t_0+t𝒩[𝐮(τ)] dτ,

induces the corresponding discrete-time dynamical system

    𝐮^n+1 = 𝒩_t(𝐮^n).



Consider a state 𝐮 on a smooth J-dimensional manifold ℳ, whose dynamics are described by (<ref>). The Koopman operator 𝒦 is an infinite-dimensional linear operator that acts on all observable functions g: ℳ→ℂ such that

    𝒦 g(𝐮) = g[𝒩(𝐮)].

The discrete-time Koopman operator 𝒦_t for the discrete dynamical system (<ref>) is defined as

    𝒦_t g(𝐮^n) = g[𝒩_t(𝐮^n)] =g(𝐮^n+1).


A major benefit of the Koopman operator is the transformation of the finite-dimensional nonlinear problem (<ref>) in the state space into an infinite-dimensional linear problem (<ref>) in the observable space. Since 𝒦_t is an infinite-dimensional linear operator, it has an infinite number of eigenvalues {λ_k}_k=1^∞ and eigenfunctions {ϕ_k}_k=1^∞.  In practice, one has to make a finite approximation of the eigenvalues and eigenfunctions. The following assumption is essential to the finite approximation and the choice of observables:


Consider a vector of P observables 𝐲,

    𝐲^n = 𝐠(𝐮^n) = [ g_1(𝐮^n);        ⋮; g_P(𝐮^n) ],    
    g_p: ℳ→ℂ ,
 
Let 𝐠 be restricted to an invariant subspace spanned by eigenfunctions of the Koopman operator 𝒦_t. 

Under this assumption, 𝐠 induces a linear operator 𝐊 that is finite-dimensional and advances these eigen-observable functions on this subspace <cit.>. The physics-aware DMD Algorithm <ref> <cit.> can be applied to approximate the eigenvalues and eigenfunctions of 𝐊 from snapshots data collected in the observable space. 

There is no principled way to choose the observables without expert knowledge of a system under consideration. Selection of observables remains a grand challenge and an active research area, e.g., machine learning and deep learning techniques were recently employed to identify the underlying Koopman operator <cit.>. 

In the context of conservation laws with shocks, the equivalency between hodograph transformation and the Koopman operator, established in this study, facilitates a “smart" choice of the observables. It is implemented via the following algorithm.

Physics-Aware DMD algorithm


  0. Create data matrices of (M-1) observables, 𝐘_1 and 𝐘_2,

    𝐘_1 = [     |     |     |;   𝐲^1   𝐲^2     ⋯ 𝐲^M-1;     |     |     | ],   𝐘_2 = [   |   |   |; 𝐲^2 𝐲^3   ⋯ 𝐲^M;   |   |   | ].

Each column of these matrices is given by 

    𝐲^n = 𝐠(𝐮^n) = [ g_1^n; g_2^n ],

where g_1^n =  x(n Δ t, 𝐮^n) is the inverse function of u(t,x), and g_2^k is a problem-dependent recording of shock information. 


  1. Apply SVD: 𝐘_1 ≈𝐔Σ𝐕^*, with 𝐔∈ℂ^P× r, Σ∈ℂ^r× r, 𝐕∈ℂ^r× M, and r denoting the truncated rank chosen by certain criteria.

  2. Compute 𝐊̃=𝐔^*𝐗'𝐕Σ^-1 as an r× r low-rank approximation for 𝐊.

  3. Compute eigen-decomposition of 𝐊̃: 𝐊̃𝐖 = 𝐖Λ, Λ = (λ_k).

  4. Reconstruct eigen-decomposition of 𝐊. Eigenvalues are Λ and eigenvectors are Φ  = 𝐔𝐖.

  5. Future 𝐲_DMD^n+1 is predicted by

    𝐲_DMD^n+1 =ΦΛ^n+1𝐛,    n>M

with 𝐛 =Φ^-1𝐲_1.

  6. Transform from observables back to state-space:

    𝐮_DMD^n =𝐠^-1(𝐲_DMD^n)= g_1^-1(𝐱_DMD^n).





Numerically, g_1 can be obtained by interpolation from a uniform mesh in the (x,u) plane to a uniform mesh in the (u,x) plane, and so can g_1^-1. The monotonicity assumption <ref> ensures that the observable functions are one-to-one maps.



The challenge of incorporating the shock information into the Lagrangian DMD algorithm <cit.> is the dependence of shock speed on the dependent variable u. Hodograph transformation facilitates the incorporation of this  implicitly nonlinear information by turning u an independent variable and by rendering the shock speed given by the Rankine-Hugoniot condition  linear and, in fact, constant.



For problems with shocks, one needs to collect snapshots until and after a shock forms. Otherwise, the Koopman operator cannot learn the shock dynamics.



For mixed wave problems, one needs to collect snapshots until and after all forms of propagation occurs. This requires pre-observation, pre-processing and understanding of the data. General initial data u_0(x) has to be separated into monotonic sub-regions. Physical quantities, such as shock speed and intersection point of shock and rarefaction propagation, must be understood from given data features. They give an explicit form of the shock observable function g_2; although problem-dependent, all the shock information is linear with respect to u.





§ NUMERICAL TESTS



We apply the physics-aware DMD to construct ROMs of scalar conservation laws in different scenarios, including a shock, rarefaction and a mixture of both. The conservative first-order upwind scheme (<ref>) is employed as reference solution, except when an analytical solution is available. The rank truncation criterion (<ref>) with ε = 10^-4 is used in all cases.





 §.§ Riemann problem for Burgers Equation with Shock


Consider the inviscid Burgers equation (<ref>) defined for (x,t) ∈ [-0.5,1.5]× [0,1] and with initial data

    u_0(x) = {    2        -0.5 ≤ x < 0
        0        0 ≤ x ≤ 1.5.
    .


This problem admits an analytical solution

    u(x,t) = {   2             -0.5 ≤ x < st
       0             st < x ≤ 1.5,
    .

where the shock speed s = 1 is determined from the Rankine-Hugoniot condition. The data needed for our DMD algorithm comes from the temporal snapshots of u(x,t) in (<ref>). 

The discontinuous initial data u_0(x) in (<ref>) do not satisfy Assumption <ref>. Thus, we approximate the step function u_0(x) with a smooth function, e.g., the hyperbolic tangent

    u_0≈ 1- tanh(x/δ),     δ≪ 1,

which satisfies Assumption <ref>. In the (u,x) plane,

    x_0≈δ/2log(2-u_0/u_0),     δ≪ 1.


This approximation is valid in the neighborhood of the shock interface. Away from it, (<ref>) is used. Snapshots of x(t,u) on a uniform mesh of u, which consists of J = 2000 equidistant points, are collected at M = 250 times until T=0.25. The ROM is used to predict the solution u(x,t) for larger times, t > T.



Figure <ref> demonstrates the the physics-aware Lagrangian DMD algorithm with hodograph transformation captures the behavior of the shock propagation. Only r=2 modes are needed to construct the ROM, which remains accurate for relatively long time in the extrapolation mode. Hodograph transformation converts the nonlinear conservation law (<ref>) with discontinuous initial data (<ref>) into a linear shift with constant speed, which is readily learned from data.




 §.§ Riemann Problem for Burgers Equation with Rarefaction Wave


Consider the inviscid Burgers equation (<ref>) defined for (x,t) ∈ [-1,1]× [0,1] and with initial data

    u_0(x) = {    -1        -1 ≤ x < 0
        1        0 ≤ x ≤ 1.
    .

This problem admits an analytical solution in the form of a rarefaction wave,

    u(x,t) = {   -1             -1 ≤ x < - t
       x/t            -t<x<t
       1             t < x ≤ 1.
    .


A hyperbolic-tangent approximation analogous to (<ref>) is used to deal with the discontinuity in the initial data u_0(x). And the same structure of data matrix is used in the DMD algorithm with J=2000 and M=250 until t=0.25.



Figure <ref> shows the same satisfactory ROM results for this problem with a rarefaction wave. Only r=2 modes are needed to obtain  accurate predictions. As the hodograph transform, x = x(t, u) satisfies a linear ODE, with u acting as an independent variable. Given an accurate approximation of the initial discontinuity, i.e., selecting δ to be sufficiently small, the ROM trained on the data from generated with the analytical solution (<ref>) is even more accurate than the HFM solution. The upwind scheme (<ref>) has first-order accuracy, O(Δ t), while the Lagrangian DMD algorithm can have spectral accuracy. Visually (in Figure <ref> and Figure <ref>), the DMD solution has a much sharper interface than that estimated with the first-order upwind scheme (<ref>). 




 §.§ Smooth Solution of Riemann Problem for Burgers Equation with Nonmonotonic Initial Data


Consider the inviscid Burgers equation (<ref>) defined for (x,t) ∈ [0,2π]× [0,1] and with initial data


    u_0(x) =1+sin(x).


Since these initial data violate Assumption <ref>, we decompose the interval [0, 2π] into two parts: in the left part, x ∈ [0, π], u_0(x)  monotonically increases; in the right part, x ∈ [π, 2π], it monotonically decreases. Each part has a unique inverse function of x_0(t,u_0). Since the shock formation time is infinite, the equation of characteristics for this Riemann problem is equivalent to (<ref>) on any finite-time interval [0,T]. Although this is a shock-free scenario, the two parts have different wave propagation behaviors. The numerical scheme (<ref>) with J = 2000 spatial discretization points and N =1000 time steps provides the reference solution. The data used to inform our DMD method consist of M = 250 snapshots of this solution.














Figure <ref> demonstrates the ability of the ROM based on our physics-aware DMD algorithm to capture these nonlinear dynamics. The ROM was trained on the early (t ≤ 0.25) data, which exhibit smooth gradients. Yet, it is capable of accurately predicting sharp gradients at later times, e.g., t=1. That is because, in the (u,x) domain of the hodograph transform, higher gradients of u(·,x) translate into flatter horizontal plots of x(·, u).








 §.§ Riemann Problem for Burgers Equation with Rarefaction and Shock


Consider the inviscid Burgers equation (<ref>) defined for (x,t) ∈ [0,2]× [0,1] and with the Gaussian-type initial data in (<ref>). This is the setting we used to illustrate the failure of the Lagrangian DMD in section <ref> (Figure <ref>). The numerical scheme (<ref>) with J =2000 spatial discretization points and N = 10^5 time steps provides the reference solution. The finer time discretization is needed to satisfy the CFL constraints. The data used to inform our DMD method consist of M = 3000 snapshots of this solution. These data are sufficiently rich to identify the rarefaction and shock behavior of the solution. The decomposition of initial data u_0(x) is needed to enforce monotonicity. The increasing branch of u_0(x) is responsible for the rarefaction and its decreasing branch gives rise to the shock.



Figure <ref> shows that the physics-aware DMD based on hodograph transformation provides  an accurate ROM for this Riemann problem, which could not be treated with the original physics-aware DMD. The physical shock  information, which is needed for the observable function g_2, includes the shock speed and the intersection point of the rarefaction wave and the shock trajectory. In this setting, the shock speed varies with time but is  still linear with respect to u. The physics-aware DMD algorithm can learn this linear relationship from the data with no difficulties. All of the advantages of linearity are achieved with the hodograph transform.



 §.§ Riemann problem for Buckley-Leverett Equation


In the last numerical experiment, we consider the hyperbolic conservation law (<ref>) with a nonmonotonic flux function,


    F = u^2/u^2+a(1-u)^2,    a=0.5,

that is defined for (x,t) ∈ [0,2]× [0,0.5] and is equipped with initial data

    u_0(x) = {    1        0 ≤ x <1
        0        1 ≤ x ≤ 2.
    .

Similar to (<ref>), the initial discontinuity is approximated with the hyperbolic tangent function. The hodograph treatment of this more general problem is provided in appendix <ref>.  The numerical scheme (<ref>) with J =2000 spatial discretization points and N = 1000 time steps provides the reference solution up to t = 0.5. The data used to inform our DMD method consist of M = 250 snapshots of this solution until t = 0.125. This set of snapshots is sufficiently rich to reveal a self-similar structure of the solution. 

Although the initial data u_0 are monotonic, their decomposition is needed according to the convex hull construction of the flux function (appendix <ref>). The reformulation involves two branches of different linear equations with two sets of the disjoint initial data. The shock observation function g_2 comprises the shock speed as well as the intersection point of the rarefaction wave and the shock trajectory. This intersection point defines the magnitude of the shock and informs the convex hull construction of the flux function. 



Figure <ref> demonstrates that the physics-aware DMD accurately captures the future states in relatively long time. Hodograph transformation allows one to determine the underlying linear Koopman operator in the nonlinear conservation laws. The iteration-free feature of DMD enhances its effectiveness and efficiency.




§ SUMMARY AND CONCLUSIONS



The Lagrangian physics-aware DMD <cit.> provides a robust tool to construct ROMs of hyperbolic conservation laws, a class of problems for which standard (Eulerian) DMD methods fail. However, this algorithm is limited to problems that admit smooth strong solutions. We extended it to problems with shocks and rarefaction waves, thus addressing a long-standing challenge in ROM construction. 
This challenge stems from sever grid distortion typical of Lagrangian POD and DMD algorithms. Lacking information about shocks and discontinuities, DMD mode projection from the HFM to a ROM does not preserve the topological structure of the interface where characteristic lines cross each other. We resolved this issue by combining hodograph transformation with physics-aware DMD algorithm <cit.>. 

Hodograph transforms are consistent with the Koopman operator theory in that both aim to identify linear structures in the underlying nonlinear dynamics. Our Lagrangian physics-aware DMD algorithm enhanced by hodograph transformation is capable of predicting the dynamics of weak solutions, which satisfies the entropy condition.  We demonstrated the accuracy and robustness of our algorithm on several numerical tests.

To the best of our knowledge, our study is the first to establish a connection between hodograph transformation and the Koopman operators. By providing a principled way for identifying the observables needed by the Koopman operator theory, this connection opens a door to construct ROMs for a wide range of nonlinear PDEs that are linearizable by hodograph transformation <cit.>. There is an algorithmic method to do the linearization via extended hodograph transforms. As a result, one can take advantage of the linearity and design robust iteration-free physics-aware DMD. Moreover, data-driven modeling and uncertainty quantification can be further explored using this framework. Our numerical experiments demonstrated that many physical quantities, such as the shock speed in Burgers equation and the mobility constant in Buckley-Leverett equation, can be learned from (simulation) data as long as one analyzes them in a “smart" way.



§ ACKNOWLEDGEMENTS


This research was supported in part by Air Force Office of Scientific Research under award numbers FA9550-17-1-0417 and FA9550-17-1-0417, and by a gift from Total.




§ SCALAR CONSERVATION LAWS WITH CONVEX FLUXES



Burgers' equation has a monotonically increasing flux function. Here, we extend our analysis to smooth, strictly convex flux functions F(u). We consider a hyperbolic conservation law (<ref>) defined for (x,t) ∈ℝ× [0,T]. It is subject to the initial condition u(x,0) = u_0(x), where the initial data u_0(x) satisfy the following assumption.


The real-valued function u_0(x) is such that


  * lim _x →±∞ u_0(x)=∓ 1, and

  * u_0(x) is non-increasing and, therefore, the inverse function x(u_0) is well-defined on -1 ≤ u_0 ≤ 1.




The domain of definition, x ∈ℝ, can be generalized to a finite-length interval (u_R,u_L). The derivation is similar.





 §.§ Solution Before Shock Formation


Similar to section <ref>, hodograph transformation yields an equation for x(t,u):

    d x/dt(t, u) = f(u),    x(0, u) = x_0(u);      u ∈(-1,1).

The convexity of F(u) ensures that its derivative f(u) is an increasing function. Let G denote the inverse function of f:
 
    G[f(u)] = f[G(u)] = u.

Then, defining y(t,u) = x[t,G(u)], (<ref>) becomes


    d y/d t (t, u) 
    =u,    y(0,u) = y_0(u) = x_0[G(u)];      u ∈(-1,1).


Differentiating both sides of this equation with respect to u,

    ∂^2 y/∂ t ∂ u = 1,

which gives

    d u/d t (t, u) = y_0^'(u)+t.


Therefore the shock formation time is determined by

    t^*=-min_u y_0^'(u)=-y_0^'[f(u^*)]





 §.§ Solution After Shock Formation


The shock speed s is given by the Rankine-Hugoniot condition,

    s=F(u_1)-F(u_2)/u_1-u_2,

where u_1(t) and u_2(t) are defined as the limits of u(t) from the top and bottom of the shock, respectively. Since s = d x^* / dt, this gives an equation for the shock trajectory x^*(t),


    d x^*/d t = F(u_1)-F(u_2)/u_1-u_2.


A system of coupled ODEs for u_1(t) and u_2(t) is derived in <cit.>,


    d u_1/d t = F_1(u_1, u_2) ≡1/g(u_1) - f^'(u_1) t [f(u_1) -  F(u_1) - F(u_2) / u_1 - u_2] 
       d u_2/d t = F_2(u_1, u_2) ≡1/g(u_2) - f^'(u_2) t[f(u_2) - F(u_1) - F(u_2)/u_1 - u_2].

where g(u) = -x'_0(u). These ODEs are subject to initial conditions u_1(t^*) = u^* and u_2(t^*) = u^*.



 §.§ Summary of Hodograph Solution

In summary, the reformulation for general scalar conservation law with convex flux is 

    {   t<t^*:       
       t>t^*:      {               u∈(u_R,u_2)∪(u_1,u_L),
                   u∈(u_2,u_1).
    .
    .

where t^* = -x'_0(u^*).


One can show that u_1(t) is monotonically increasing in time and u_2(t) is monotonically decreasing, so that

    u_1 ≥ u^*,    u_2≤ u^*,    x'_0(u_1)+t ≤ 0,    x'_0(u_2)+t ≤ 0.

In many cases of interests, and in our numerical experiment, either u_2 = u_R and u_1 = u_L or |u_2 - u_1| ≪Δ t (so that u_2 ≈ u_R and u_1 ≈ u_L). This allows one to focus on shock propagation, i.e., on (<ref>), without having to solve (<ref>).







For more general initial condition u_0, one needs to decompose u_0(x) into regions of monotonicity. Each monotonic piece of u_0 would have a unique inverse function x_0(u_0). Then, based on the generalized entropy condition, one constructs the convex hull for the flux function F(u), providing a way to decompose the initial data. Shock propagating initial data and rarefaction propagating initial data are determined afterwards. Then,  the full solution is the combination of the rarefaction pieces and the shock pieces.








