Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

rspa.royalsocietypublishing.org

Research Cite this article: Gouasmi A, Parish EJ, Duraisamy K. 2017 A priori estimation of memory effects in reduced-order models of nonlinear systems using the Mori–Zwanzig formalism. Proc. R. Soc. A 473: 20170385. http://dx.doi.org/10.1098/rspa.2017.0385 Received: 2 June 2017 Accepted: 1 September 2017

Subject Areas: mathematical modelling, applied mathematics, computational mathematics Keywords: Mori–Zwanzig formalism, reduced-order modelling, closure modelling, orthogonal dynamics Author for correspondence: Ayoub Gouasmi e-mail: [email protected]

A priori estimation of memory effects in reduced-order models of nonlinear systems using the Mori–Zwanzig formalism Ayoub Gouasmi, Eric J. Parish and Karthik Duraisamy Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI, USA AG, 0000-0003-0647-6723 Reduced models of nonlinear dynamical systems require closure, or the modelling of the unresolved modes. The Mori–Zwanzig procedure can be used to derive formally closed evolution equations for the resolved physics. In these equations, the unclosed terms are recast as a memory integral involving the time history of the resolved variables. While this procedure does not reduce the complexity of the original system, these equations can serve as a mathematically consistent basis to develop closures based on memory approximations. In this scenario, knowledge of the memory kernel is paramount in assessing the validity of a memory approximation. Unravelling the memory kernel requires solving the orthogonal dynamics, which is a high-dimensional partial differential equation that is intractable, in general. A method to estimate the memory kernel a priori, using full-order solution snapshots, is proposed. The key idea is to solve a pseudo orthogonal dynamics equation, which has a convenient Liouville form, instead. This ersatz arises from the assumption that the semi-group of the orthogonal dynamics is a composition operator for one observable. The method is exact for linear systems. Numerical results on the Burgers and Kuramoto–Sivashinsky equations demonstrate that the proposed technique can provide valuable information about the memory kernel.

1. Introduction Complex dynamical systems are typically prohibitively expensive to solve using direct simulation approaches, 2017 The Author(s) Published by the Royal Society. All rights reserved.

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

φ(0) = x,

(1.1)

where φ = φ(x, t) ∈ RN is the vector of modal coefficients. A small subset φˆ = (φ1 , . . . , φm ) ∈ Rm , m  N, of this system is then evolved in time. The underlying assumption is that, at every time instant, the essential or dominant dynamics are captured within these m modes. The construction of accurate reduced models for nonlinear problems is, however, rendered challenging by the interaction between the resolved and unresolved physics. In many cases, the evolution equation of φˆ involves the discarded modes φ˜ = (φm+1 , . . . , φN ) ∈ RN−m . The process of modelling this contribution is referred to as closure. This is an outstanding challenge that reduced-fidelity approaches such as large eddy simulation (LES) [12,13] must also deal with. This work is inspired from the optimal prediction framework developed by the Chorin group [14–17]. This framework is a reformulation of the Mori–Zwanzig (M-Z) formalism [18,19] of non-equilibrium statistical mechanics. The M-Z approach provides a mathematically consistent framework to express reduced-order models of dynamical systems. A key element is the use of projection operators P in phase space. Functions f (x) of the initial state are projected onto a ˆ = xˆ only. The formalism consists of separating the subspace of functions of the resolved part φ(0) state variables into a resolved and unresolved set. The nonlinear dynamical system is recast as a Liouville equation; a linear PDE in the space of the initial conditions. A generalized Langevin equation (GLE) that describes the evolution of the resolved variables as a function of the resolved variables only can then be derived. In this setting, the effect of the unresolved modes on the resolved modes appears as a non-local memory integral and a noise term. The M-Z approach—by itself—does not lead to a reduction in computational complexity as the evaluation of the memory and noise terms requires the solution of the orthogonal dynamics equation, but it does provide a starting point for the development of closure models [20,21]. It is important to emphasize that, while this work is inspired from the Chorin framework, and the M-Z formalism by transitivity, the authors are considering applications in a specific context. The M-Z procedure originated in the non-equilibrium statistical mechanics community, where the goal is to solve for probability density functions and time correlation functions of non-equilibrium systems [22]. Originally, the procedure was limited to Hamiltonian dynamical systems. Chorin extended this result to general time-dependent and non-Hamiltonian systems [14, ch. 9], leading to applications to problems such as the viscous Burger’s equation (VBE) [23] and the Euler equations [24]. Chorin’s framework was developed for optimal prediction, the goal of which is to approximate the solution of nonlinear time-dependent problems such as equation (1.1) where a full-order solution is not affordable and the unresolved part x˜ of the initial conditions is uncertain. Chorin and co-workers developed a method to compute mean solutions with respect to initial probability distributions. In the same vein, Venturi, Cho and Karniadakis considered the M-Z formalism for uncertainty quantification [25–27]. This work considers nonlinear, non-Hamiltonian, time-dependent systems in which the initial conditions ˜ = x˜ = 0, with no uncertainty. are fully resolved, i.e. φ(0) In optimal prediction, the projection operator is a conditional expectation (Pf )(x) = E[f (x) | xˆ ]. It can also be approximated by a finite rank projection. In our context, the projector is a simple truncation (Pf )(ˆx, x˜ ) = f (ˆx, 0). It has been used in [20,28]. Although the M-Z approach is independent of the projector, its viability is dependent upon the well-posedness of the orthogonal dynamics equation [29], which depends on the projector. Givon et al. [29] studied the existence of solutions to the orthogonal dynamics for Hamiltonian systems. They proved the existence of

...................................................

d φ(t) = R(φ(t)); dt

2

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

insofar as real-world science and engineering applications are concerned. This fact has motivated the development of manifold reduced-order modelling (ROM) techniques based on proper orthogonal decomposition [1–4], balanced proper orthogonal decomposition [5,6], balanced truncation [7,8], the reduced basis method [9,10] and Krylov subspaces [11]. In projection-based ROM of partial differential equations (PDEs), the governing equations are projected onto an alternate space of solutions, in which the state of the system can be represented by an equivalent set of distinct modes. This typically results in a system of ordinary differential equations (ODE):

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a) Model reduction and memory As an introductory example [22], consider a linear dynamical system:

and

⎫ d ˜⎪ ⎪ φˆ = A11 φˆ + A12 φ, ⎪ ⎪ dt ⎪ ⎬ d 21 22 φ˜ = A φˆ + A φ˜ ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎭ φ(0) = x,

(2.1)

with A11 ∈ Rm×m , A12 ∈ Rm×(N−m) , A21 ∈ R(N−m)×m and A22 ∈ R(N−m)×(N−m) . To solve directly for ˆ one has to model the ‘unclosed’ term A12 φ˜ as a function of φˆ only. This can be accomplished by φ, ˆ first solving the equation for φ˜ assuming knowledge of φ: ˜ = eA22 t x˜ + φ(t)

t 0

22

eA

(t−s)

ˆ ds. A21 φ(s)

(2.2)

...................................................

2. The Mori–Zwanzig formalism

3

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

classical solution for finite-rank projections and the existence of weak solutions for the conditional expectation. Even though the truncation projector can be seen as a limiting case of the conditional expectation, non-Hamiltonian systems are considered in this work. Solution techniques for the orthogonal dynamics have been developed for finite-rank projection operators. One of the most notable has been devised by Chorin [14]. The technique involves the decomposition of the solution using a finite-rank orthonormal basis and the evolution in time of a set of Volterra equations for the basis coefficients. One of the main issues with this approach lies in the high expense associated with the inner product and the number of basis functions required for the decomposition to be representative. Twenty-one basis functions were used for the 4-d.f. Hamiltonian system studied in [14]. For the Burgers equation, Bernstein [23] showed that unless the number of resolved modes is very small, a similar approach for the truncation projection operator would not be computationally tractable. The number of basis functions needed would typically escalate for higher-dimensional problems. In a molecular dynamics setting, Darve [30] developed a discrete version of the M-Z formalism and applied it to compute the kernel in GLEs of Hamiltonian systems. This methodology was developed for the conditional expectation projector, and it was shown that the solution to the discrete orthogonal dynamics is the same as that of the continuous version in the asymptotic limit. Darve also points out the high computational expense if the dimension of the resolved space is large. At this juncture, it is important to emphasize that the M-Z formalism does not lead to an ROM directly, as it does not address the issue of how to evaluate the outstanding memory term in a tractable fashion, or whether it is even possible. A priori knowledge of the memory kernel can help answer these two questions and is the primary focus of this work. A method is developed to numerically estimate the memory kernel a priori for nonlinear systems. The key idea is to consider an ersatz of the orthogonal dynamics equations that can be solved using the method of characteristics. This ersatz arises from the assumption that the semi-group generated by the orthogonal dynamics is a composition operator. Furthermore, it will be shown that evaluating the memory kernel amounts to computing sensitivities of the solution to the orthogonal dynamics with respect to the initial conditions in a direction that is imposed by the ROM solution. The remainder of this paper will be organized as follows: §2 will introduce the M-Z formalism. In §3, an approach to estimate the memory kernel will be presented. In §4, our procedure will be applied to the Brusselator, the VBE and Kuramoto–Sivashinsky (K-S). Conclusions and perspectives will be provided in §5.

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

By injecting equation (2.2) in the unclosed term, one obtains ˆ ds + A12 eA t x˜ . A12 e(t−s)A A21 φ(s) 22

22

(2.3)

Three distinct terms feature in this equation: a Markovian term, a memory involving the past history of the resolved physics and a term that encapsulates the influence of the initial conditions. In an optimal prediction context, the latter term is interpreted as noise produced by the uncertainty in the initial conditions. In our context, there is no noise contribution because the initial condition is fully resolved x˜ = 0. The closure term is exclusively memory. Equation (2.3) is often called a GLE. In equation (2.3), the dependency on φ˜ has been removed without introducing any approximation. However, equation (2.3) is not a ROM because of the inherent cost. It requires multiple evaluations of the exponential of a (large) (N − m) × (N − m) matrix multiplied by a (N − m) × m matrix. The cost of evolving this ROM in time should be comparable to that of ˆ computing A11 φ. To make equation (2.3) a genuine ROM, the memory term has to be approximated. Indeed, 22 ˆ which is called the memory the approximation depends on the integrand A12 e(t−s)A A21 φ(s), kernel. Let us assume that the entire set of eigenvalues of A22 is strictly negative. Then a relevant assumption to make is that the memory has finite support, i.e. there exists some memory length τ > 0 such that: t t 22 22 ˆ ds ≈ ˆ ds. A12 e(t−s)A A21 φ(s) A12 e(t−s)A A21 φ(s) 0

t−τ

This is highlighted in figure 1 for matrices (A11 , A12 , A21 , A22 ) = (−0.8, 1, 1, −1). The initial conditions are (ˆx, x˜ ) = (1, 0). The subgrid term A12 φ˜ at time t is equal to the area of the shaded 22 ˆ over the elapsed time range s ∈ [0, t]. In these snapshots, the yellow region A12 e(t−s)A A21 φ(s) memory does appear to be finite. In addition, the decay profile of the kernel, from its most recent (s = t) to its oldest (s = 0) contributions, suggest that the subgrid/memory term could be reasonably approximated using the most recent kernel contribution and a good guess of the memory length τ . In other words, if the kernel has a simple decay profile (exponential, for instance) one could consider a closure model of the form: ˜ = A12 φ(t)

t 0

ˆ ds ≈ A12 C(τ )A21 φ(t), ˆ A12 e(t−s)A A21 φ(s) 22

with C(τ ) being either a scalar or a (N − m) × (N − m) matrix, depending on the level of sophistication desired. The memory length τ does not have to be the same for all m components of the memory kernel, nor does it have to be constant throughout the run. If C(τ ) is simple enough, the evaluation of the memory approximation can be made reasonable enough for use as a ROM closure model. It is important to appreciate the fact that the above discussion on developing closure models based on memory was straightforward only because rewriting equation (2.1) as a GLE, with all outstanding terms explicitly known, was possible. As one can anticipate, however, for nonlinear systems, additional challenges arise. A first nonlinear ODE considered in this work is the Brusselator [31]:

and

⎫ d ⎪ φ1 = A − Bφ1 − φ1 + φ12 φ2 ⎪ ⎬ dt ⎪ d ⎪ ⎭ φ2 = Bφ1 − φ12 φ2 , dt

(2.4)

where φˆ = φ1 , φ˜ = φ2 , m = 1 and N = 2. A and B are constants. The Brusselator will serve as a toy problem to familiarize the reader with the concepts and the method presented in this work.

...................................................

0

4

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

d φˆ = A11 φˆ + dt

t

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

(b)

1

2

3

4

5 s

6

7

8

9 10

(c) 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0

0

1

2

3

4

5 s

6

7

8

9 10

0

1

2

3

4

5 s

6

7

8

9 10

(d)

0

1

2

3

4

5 s

6

7

8

9 10

1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0

ˆ to the memory at t is exactly Figure 1. Linear system: snapshots of the memory kernel. ∀ s ∈ [0, t], the contribution of φ(s) 12 A22 (t−s) 21 ˆ A φ(s). Memory appears finite, and the decay profile of the kernel suggests that the memory can be approximated A e ˆ and a measure of τ only. (a) t = 1.9 s, (b) t = 3.8 s, (c) t = 5.7 s and (d) t = 9.6 s. using the contribution of φ(t)

It is possible to write a GLE for the evolution of φˆ for nonlinear systems. For instance, it is possible to rewrite the subgrid term φˆ 2 φ˜ in the Brusselator as t ˆ 2 φ(t) ˜ = K1 (φ(s), ˆ φ(t) t − s) ds + F1 (x, t), 0

with F1 (x, t) dropping if x˜ = 0. Further details are provided in the next sections. The technique used to obtain the above GLE is based on projection operator methods [32] and was derived by Mori [18] and Zwanzig [22] for non-equilibrium, Hamiltonian systems. Chorin [14] extended it to general dynamical systems. In this work, important concepts of the M-Z formalism are discussed before presenting the derivation [14,16] of the GLE in the general case. In the case of nonlinear systems, very little is known about the memory kernel (contribution of ˆ only) and the noise (value at t = 0 only). More details are provided in the following sections. φ(t) In short, memory approximations are hard to devise in this case because of the outstanding challenge posed by an auxiliary orthogonal dynamics equation. The goal of the present work is to develop a numerical method to estimate the memory kernel for nonlinear systems.

(b) Liouville equivalence The Liouville operator associated with equation (1.1) is L=

N  i=1

Ri (x)

∂ . ∂xi

...................................................

contribution of fˆ(s) to memory fˆ A12f˜ = memory

0

5

1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

Consider the following linear PDE in phase space:

6 (2.5)

where g(x) is a function of the initial conditions. Equation (2.5) is the Liouville equation. Chorin showed that, for any regular function g(x), the solutions of equations (1.1) and (2.5) satisfy: u(x, t) = g(φ(x, t)).

(2.6)

We refer to this result as the Liouville equivalence. In the case g(x) = x, the solution to the Liouville equation is the same as that of equation (1.1). The Liouville equivalence demonstrates that the solutions u(x, t) of equation (2.5) for all functions g(x) are known once the solution φ(x, t) to equation (1.1) is known. This is an instance of the method of characteristics for advection equations. This equivalence is a fundamental step in deriving the GLE (§2d). In semi-group notation, the solution of the Liouville equation can be rewritten as: u(x, t) = etL u(x, 0) = etL g(x), where etL is the semi-group generated by L. The Liouville equivalence endows this semi-group with a remarkable composition property: etL g(x) = g(etL x).

(2.7)

This property is of central interest to this work. As a side note, the semi-group etL is also known as the Koopman semi-group [33,34] associated with the continuous-time dynamical system (1.1).

(c) Projection operators Projection operators [32] are central in non-equilibrium problems and the M-Z formalism. For purposes of introduction, the optimal prediction setting of Chorin [14] will be initially assumed. Consider the initial conditions x to be random variables drawn from some probability distribution μ(x) and a probability density function ρ(x). Denote Γ , the vector space in which x and φ(t) lie (RN in our case). For a function g(x) that operates on Γ , define the expected value of g by  g(x)ρ(x) dx. E[g] = Γ

Denote L2 the Hilbert space of functions on Γ endowed with the inner product (f , g) = E[fg]. In the M-Z formalism, a projector P that transforms functions in L2 into functions of the resolved component xˆ ∈ Rm , is used. There are several different projectors [14] for which the M-Z formalism has been considered: (i) In irreversible statistical mechanics, the linear projection is given by ∀ f ∈ L2 , (Pf )(ˆx) =

m 

a−1 ij (f , xi )xi ,

i,j=1

where a−1 ij are the entries of a matrix whose inverse has entries aij = (xi , xj ). (ii) In optimal prediction, the conditional expectation of f given xˆ , is given by  f (x)ρ(x) d˜x 2 , d˜x = dxm+1 . . . dxN . ∀ f ∈ L , (Pf )(ˆx) = E[f |ˆx] =  ρ(x) d˜x This is also the orthogonal projection onto the span of all functions of xˆ .

...................................................

u(x, 0) = g(x),

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

∂ u(x, t) = Lu(x, t); ∂t

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(iii) The conditional expectation can be approximated by a finite-rank projection onto the span of an orthonormal set of function hk (ˆx), k = 1 . . . K:

k=1

Define Q = I − P the complimentary projector, where I is the identity operator. In this work, we focus on under-resolved simulations of nonlinear time-dependent systems with fully resolved initial conditions (i.e. x˜ = 0). The projector to be used for that set-up is a truncation: (Pf )(ˆx) = f (ˆx, 0) = f (ˆx).

(2.8)

This projector can be seen as a conditional expectation in the limit where the probability density function is a delta function centred at x˜ = 0.

(d) Generalized Langevin equation The Liouville equation (2.5) for g(x) = xj is given by ∂ tL e xj = L etL xj = etL Lxj . ∂t The r.h.s. term of the above equation can decomposed using I = P + Q: ∂ tL e xj = etL PLxj + etL QLxj . ∂t

(2.9)

The first term etL PLxj is a function of φˆ only. The point of recasting the dynamical system (1.1) as a linear PDE is that, using the Dyson formula [35]: t etL = etQL + e(t−s)L PL esQL ds, 0

the second term etL QLxj in equation (2.9) that depends on the full solution φ ∈ RN can be rewritten as t etL QLxj = etQL QLxj + e(t−s)L PL esQL QLxj ds. 0

The first term Fj (x, t) = etQL QLxj models the influence of the initial conditions. It is often interpreted as the noise produced by the uncertainty in the initial conditions. The second term is the general expression for the memory, which involves the past history of resolved physics ˆ φ(s) only. Ultimately, the dynamical system (1.1) has been recast as a set of GLEs for the state components φj (x, t) = etL xj : ∂ tL e xj = etL PLxj + ∂t

t 0

e(t−s)L PL esQL QLxj ds + etQL QLxj ,

(2.10)

or, in a more compact form: ∂ ˆ t)) + φj (x, t) = Rj (φ(x, ∂t Rj (ˆx) = (PRj )(ˆx),

t 0

ˆ s), t − s) ds + Fj (x, t), Kj (φ(x,

Fj (x, t) = etQL QLxj ,

Kj (ˆx, t) = PLFj (x, t).

ˆ t) only. The second In the r.h.s., the first term is the Markovian contribution that depends on φ(x, term is the memory that will involve contributions from all past values of the resolved physics ˆ s), s ∈ [0, t] only. The relationship between the memory kernel Kj and the noise Fj is referred to φ(x, as a ‘fluctuation dissipation theorem’ [14–16,22,36] in the original context of statistical mechanics.

...................................................

K  (f , hk )hk (ˆx).

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

∀ f ∈ L2 , (Pf )(ˆx) =

7

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

Another consequence of the projection operator used is the form assumed by the integrand. In [18,19,22], the integrand can be separated into a kernel and a time history partly because of the projection operator used. Such a definition of the kernel and associated terminology is not ˆ x, s), t − s) as the kernel. possible in the current context. Thus, we refer to the entire term Kj (φ(ˆ As explained in §1a for the linear case, equation (2.12) is not a reduced-order model. This is because of the intrinsically high cost associated with the evaluation of the kernel (despite being a function of φˆ ∈ Rm only!). For instance, evaluation of the kernel in the linear case involved operating on matrices of size much larger than m. In the nonlinear case, evaluating the kernel requires the solution of the orthogonal dynamics equation.

(e) Orthogonal dynamics F(x, t) = ( Fj (x, t) )1≤j≤N ∈ RN is the solution of a linear PDE: ∂ F(x, t) = QLF(x, t); ∂t

F(x, 0) = QLx.

(2.13)

This equation is called the orthogonal dynamics because its solution lies in the orthogonal space of P at all instants. Givon et al. [29] studied the existence of solutions to the orthogonal dynamics for Hamiltonian systems. They proved the existence of classical solutions for finite-rank projections and the existence of weak solutions for the conditional expectation. Even though the truncation projector can be seen as a limiting case of the conditional expectation, non-Hamiltonian systems are considered in this work. Therefore, existence of solutions to the associated orthogonal dynamics cannot be supported using their work. It is possible to find a solution to the orthogonal dynamics when the dynamical system is linear. Consider equation (2.1) and define matrices A and A˜ ∈ RN×N as     A11 A12 0 A12 ˜ A= and A = . A21 A22 0 A22 The Liouville operator is given by

⎛ ⎞ N N   ∂ ⎝ Aij xj ⎠ . L= ∂xi i=1

We note that

j=1

  0 ˜ = Ax. QLx = R(x) − PR(x) = A x˜

Therefore, ∀ n ≥ 1, and it can be shown, by induction, that:

(QL)n x = A˜ n x, 

0 ∀ n ≥ 1, A = 0 ˜n

 A12 (A22 )n−1 . (A22 )n

...................................................

ˆ t) = φ(ˆ ˆ x, t). The corresponding equation is therefore: The authors choose to work with P φ(x, t ∂ ˆ x, t)) + Kj (φ(ˆ ˆ x, s), t − s) ds. (2.12) φj (ˆx, t) = Rj (φ(ˆ ∂t 0

8

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

The interpretation and implementation of each term is imposed by the projector P [36]. The ˆ t) is the target quantity the user wants to solve for. The projector is chosen such that P φ(x, evolution equation for the target is obtained by projecting the GLE: t ∂ (2.11) P etL xj = P etL PLxj + P e(t−s)L PL esQL QLxj ds. ∂t 0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

For any matrix B, the series

∞ n n n=1 (t /n!)(QL) x and:

n /n!)Bn

converges to etB . Therefore, so does the series

 ∞ n  t 0 (QL)n x = n! 0 n=1

   22 A12 etA A12 tA22 x = e x˜ 22 A22 A22 etA

⎛ ⎞  N N ∂ A12 A22 s  ⎝ x˜ LF(x, s) = Aij xj ⎠ e A22 ∂xi 

i=1

j=1

⎛ ⎞   N N 12   ∂ 22 A A12 A22 s ⎝ x˜ = = Aij xj ⎠ e e A s [A21 A22 ]x. A22 A22 ∂xi 



i=m+1

j=1

Applying the projector P, one gets   A12 A22 s 21 [A K(ˆx, s) = PLF(x, s) = e A22 Ultimately, we obtain

 22

A ]

xˆ

0(N−m)×1





 A12 A22 s 21 A xˆ . = e A22

  A12 A22 (t−s) 21 ˆ ˆ A φ(s). e K(φ(s), t − s) = A22

This is the exact expression for the memory kernel. The first m components correspond to what is in equation (2.3). As a side note, the above solution method did not assume analyticity of etQL , meaning that this semi-group could be described as etQL =

∞ n  t (QL)n . n!

n=0

In semi-group theory [37], such an expansion is convergent if the generator QL is a bounded operator on a Banach space. To find a solution for the linear case, we just showed that the series (tn /n!)(QL)n x was convergent.

3. A priori estimation of the memory Kernel This work seeks to further an approach Parish & Duraisamy [20] hinted at to solve the orthogonal dynamics and suggests a procedure to estimate the memory kernel in nonlinear dynamical systems. An equivalent form of equation (2.12) that is more fitting to the authors’ viewpoint is the following: t ˆ x, t)) = Kj (φ(ˆ ˆ x, s), t − s) ds, j = 1 . . . m. (3.1) wj = Rj (φ(ˆx, t)) − Rj (φ(ˆ 0

The l.h.s. term is the general expression for the subgrid terms, denoted wj , that need to be modelled in ROMs. The M-Z formalism shows that they can be exactly represented as memory integrals. As a reminder, the noise term dropped because we are within the setting the projector P imposes, which is fully resolved initial conditions. The goal of this work is not to derive closure models based on memory approximations as in [14,20,21,23,38,39]. No closure model has been derived or tested in ROM simulations. Full-order ˆ x, t) and the ‘exact’ simulations have been run to provide for the ‘exact’ ROM solution snapshots φ(ˆ subgrid terms wj . The focus of this work is in estimating the memory kernel. The goal is to extract the kernel as in §2a, but for nonlinear systems. The memory kernel is unknown for nonlinear systems. Equation (3.1) is therefore the only verification model available. In the numerical results section, the exact ROM snapshots and

9 ...................................................

is a solution. Then:

n=0 (t

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

F(x, t) =

∞

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

subgrid terms are first computed. Then the procedure is applied to estimate the memory kernel ˆ x, s), t − s) at the required points and time ranges. Kj (φ(ˆ

F(x, t) = etQL F(x, 0), where F(x, 0) = QLx = QR(x) = R(x) − PR(x) = R(x) − R(ˆx). Earlier, the composition property of the semi group etL was introduced: etL g(x) = g(etL x) = g(φ(x, t)). The main assumption of this work consists in granting the semi-group of the orthogonal dynamics etQL the composition property for g(x) = QLx. Defining φ Q (x, t) = etQL x ∈ RN , one has: F(x, t) = etQL F(x, 0) ≈ F(etQL x, 0) = F(φ Q (x, t), 0) = R(φ Q (x, t)) − R(φˆ Q (x, t)), ˆ By definition, φ Q (x, t) satisfies where φˆ Q is defined in the same manner as φ. ∂ Q φ = QLφ Q ; ∂t

φ Q (0) = x.

(3.2)

Using the expression φ Q (x, t) = etQL x and the fact that etQL and QL commute, one obtains: ∂ Q φ = etQL QLx = etQL QR(x) = etQL (R(x) − R(ˆx)). ∂t Using composition again, we finally obtain d Q φ (t) = R(φ Q (t)) − R(φˆ Q (t)); dt

φ Q (0) = x.

(3.3)

This is an ODE, where the r.h.s. term is F(x, t). Accordingly, Fj (x, t) can be computed by evolving in time equation (3.3) and taking the jth component of the r.h.s. Equation (3.3) results from granting the semi-group of the orthogonal dynamics a composition property. Let us introduce the Liouville operator LQ associated with that ODE: LQ =

N  i=1

 ∂ ∂ = (QRi )(x) . ∂xi ∂xi N

(Ri (x) − Ri (ˆx))

i=1

Using the Liouville equivalence for g(x) = F(x, 0) = QLx, we show that g(φ Q (x, t)) = F(x, t) satisfies ∂ F(x, t) = LQ F(x, t); ∂t

F(x, 0) = QLx.

(3.4)

Basically, the orthogonal dynamics equation has been replaced with an ersatz (or ‘pseudo orthogonal dynamics’) whose solution can be related to that of equation (3.3)—that we call pseudo orthogonal ODE—via the Liouville equivalence. Whether or not granting etQL a composition property is valid depends on how close are the solutions to the orthogonal dynamics equation (2.13) and its ersatz equation (3.4) are. For linear

...................................................

The solution of the orthogonal dynamics can be formally expressed as

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

(a) Pseudo orthogonal dynamics

10

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

systems, the pseudo orthogonal ODE is the following:

11 φ Q (0) = x.

The ersatz solution is given by the r.h.s., that is:     A12 tA22 A12 ˜ Q φ = e x˜ . F(x, s) = A22 A22 This is the same solution as the one derived in §2e. It has already been shown that this solution leads to the exact kernel. Granting etQL a composition property for the observable g(x) = QLx is therefore valid for linear systems. For nonlinear systems, however, solutions of the orthogonal dynamics are unknown. In the absence of reference nonlinear problems where the solution of the orthogonal dynamics is known, an indirect verification approach driven by equation (3.1) is required. The numerical results obtained in §4 suggest that there exists a class of nonlinear problems for which the assumption etQL g(x) = g(etQL x) for g(x) = QLx is worth considering.

(b) Computing the kernel The memory term for the jth resolved variable is given by t t t t ˆ − s), s) ds = Kj (φ(s), ˆ e(t−s)L PLFj (x, s) ds = e(t−s)L Kj (ˆx, s) ds = Kj (φ(t t − s) ds. 0

0

0

0

ˆ and t − s are known quantities and therefore the challenge lies in computing Both φ(s) Kj (ˆx, s) = PLFj (x, s). Without the application of P, one has LFj (x, s) =

N  i=1

Ri (x)

∂ Fj (x, s). ∂xi

(3.5)

These N sensitivities could be separately estimated, but it is more efficient to view LF as one ¯ single sensitivity along the unitary direction R(x) = R(x)/ R(x) : LF(x, s) = R(x)

N  i=1

¯ i (x) ∂ F(x, s) = R(x) (∇ ¯ F(x, s)) R R(x) ∂xi

= R(x) lim

ε−>0

¯ F(x + εR(x), s) − F(x, s) . ε

Applying P and using finite difference (Fréchet derivative), Kj (ˆx, s) ≈ R(ˆx)

¯ x), s) − Fj (ˆx, s) Fj (ˆx + εR(ˆ ε

= R(ˆx)

¯ x), s) Fj (ˆx + εR(ˆ ε

.

The term Fj (ˆx, s) dropped because Fj (ˆx, s) = P etQL x = 0 by definition of the orthogonal dynamics. In the numerical tests, the value of the sensitivity factor ε is taken small enough so that the converged value is attained, but not too small in order to avoid subtractive cancellation errors. Ultimately, the memory kernel is approximated as ˆ ˆ t − s) ≈ R(φ(s)) Kj (φ(s),

¯ φ(s)), ˆ + εR( ˆ Fj (φ(s) t − s) ε

,

s ∈ [0, t].

(3.6)

In other words, evaluating the memory kernel amounts to computing the sensitivity with respect to the initial conditions, of the solution to the orthogonal dynamics in a direction that ˆ is determined by the resolved component φ(s). Note that equation (3.6) and the latter statement are independent of how the orthogonal dynamics have been solved.

...................................................

d Q φ˜ = A22 φ˜ Q ; dt

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

d Q φˆ = A12 φ˜ Q ; dt

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

Once the kernel is computed, the memory follows. With a rectangle rule for instance, one can make the following approximation: Nt 

ˆ n ), t − sn )s. Kj (φ(s

n=1

(c) Summary and cost Denote T the time domain discretized into Nt equally spaced time instants T = {0, t1 , . . . , tNt }. A full-order simulation of equation (1.1) has been run, so that the ‘exact’ ROM solution snapshots ˆ Nt )} are available. Denote Cfull the cost associated with computing these Nt vector ˆ 1 ), . . . φ(t {φ(t values. Denote  tn ˆ Kj (φ(s), t − s) ds, Mn = 0

the memory term at t = tn , which is also the subgrid term at the same instant wj (tn ). Then the first five terms starting at n = 1 can be approximated with a rectangle rule in the form: ˆ 1 ), 0)], M1 = t[Kj (φ(t ˆ 1 ), t) + Kj (φ(t ˆ 2 ), 0)], M2 = t[Kj (φ(t ˆ 1 ), 2t) + Kj (φ(t ˆ 2 ), t) + Kj (φ(t ˆ 3 ), 0)], M3 = t[Kj (φ(t ˆ 1 ), 3t) + Kj (φ(t ˆ 2 ), 2t) + Kj (φ(t ˆ 3 ), t) + Kj (φ(t ˆ 4 ), 0)] M4 = t[Kj (φ(t ˆ 1 ), 4t) + Kj (φ(t ˆ 2 ), 3t) + Kj (φ(t ˆ 3 ), 2t) + Kj (φ(t ˆ 4 ), t) + Kj (φ(t ˆ 5 ), 0)]. M5 = t[Kj (φ(t

and

Thus, to compute M1 up to M5 , one must compute: — — — — —

ˆ 1 ), s) for s ∈ {0, t, 2t, 3t, 4t}. Kj (φ(t ˆ 2 ), s) for s ∈ {0, t, 2t, 3t}. Kj (φ(t ˆ 3 ), s) for s ∈ {0, t, 2t}. Kj (φ(t ˆ 4 ), s) for s ∈ {0, t}. Kj (φ(t ˆ 4 ), s) for s ∈ {0}. Kj (φ(t

ˆ n ), s) requires computing Fj (φ(t ˆ n) + According to equation (3.6), each kernel value Kj (φ(t ¯ ˆ n )), s). It is obtained—see §3a—by solving the pseudo orthogonal ODE: εR(φ(t ∂ Q φ (t) = R(φ Q (t)) − R(φˆ Q (t)); ∂t

¯ φ(t ˆ n ) + εR( ˆ n )), φ Q (0) = φ(t

(3.7)

and retrieving the first m components of the r.h.s. term at t = s. Let us assume that the costs of evolving equation (3.7) and equation (1.1) in time are the same. Then the cost CM−Z of computing the memory terms Mn at all Nt time instants scales as CM−Z ∝

Nt − 1 Cfull . 2

(3.8)

ˆ tn − s), the scaling that is obtained Since computing Mn requires computing n values of Kj (φ(s), should not be surprising. It is emphasized again that the purpose of the proposed method is to estimate the memory kernel a priori in nonlinear systems. Such scaling arguably limits the size of the problems that can be considered. Note however, ˆ n ), s) that the described procedure can be easily implemented in parallel. The kernel values Kj (φ(t at the required instants s can be computed separately for each n = 1 . . . Nt .

...................................................

0

ˆ Kj (φ(s), t − s) ds ≈

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

t

12

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

(b) w1 = memory estimated memory

2.5

1.0 0.5 0

5

10

15 t

20

25

30

0

5

10 t

15

20

Figure 2. Brusselator: the dotted line is w1 (t) = φ1 (t)2 φ2 (t) in equation (3.1) computed using full-order solution snapshots. The solid line is the same quantity but evaluated by integrating the kernel (r.h.s. of equation (3.1)) over time. (a) (A, B) = (1, 1.7) − dt = 1.5 × 10−2 s and (b) (A, B) = (1, 3) − dt = 10−2 s . (Online version in colour.)

Equation (3.8) gives the cost of computing the kernel values in the full range s ∈ [0, t]. If the size of the problem is such that CMZ is a limiting factor, one can still use the procedure to compute the ˆ t − s) in a truncated range s ∈ [t − τ , t] with τ taken such that the corresponding kernel Kj (φ(s), cost CMZ (τ ) is affordable. This can be seen as an indirect way of investigating finite memory (provided that the method is always accurate). For all the problems considered in this work, the kernel was estimated in the full range s ∈ [0, t] even though the observations that followed suggested that it was not needed.

4. Applications In this section, the effectiveness of the kernel estimation technique is assessed in three problems of increasing complexity—the Brusselator, the VBE and the K-S equation.

(a) Brusselator Two configurations of the Brusselator equation (2.4) are considered: — (A, B) = (1, 1.7), where the system is asymptotically stable. — (A, B) = (1, 3), where the system exhibits limit cycle oscillations (LCO). The Liouville operator for the Brusselator is given by L ≡ (A + x21 x2 − Bx1 − x1 )

∂ ∂ + (Bx1 − x21 x2 ) . ∂x1 ∂x2

The pseudo orthogonal ODE is given by d Q φ = (φ1Q )2 φ2Q ; dt 1

d Q φ = −(φ1Q )2 φ2Q ; dt 2

φ Q (0) = x,

(4.1)

and F1 (x, s) = φ1Q (x, s)2 φ2Q (x, s). The original ODE (2.4) and the pseudo orthogonal ODE (4.1) are solved using fourth-order Runge–Kutta time integration. The kernel K1 is then computed at the required points using equation (3.6) with ε = 10−8 . Figure 2a and b shows the numerical results for the stable and LCO configurations, respectively. There are some amplification and phase errors, but the results are good overall.

...................................................

1.5

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

2.0

0

13

18 16 14 12 10 8 6 4 2 0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

(b)

2.0

1.5

1.5

1.0

1.0

0.5

0.5 0

5

10

15 s

20

25

0

30

(c)

0

5

10

15 s

20

25

30

0

5

10

15 s

20

25

30

(d) 3.0

3.0

2.5

2.5

2.0

2.0

1.5

1.5

1.0

1.0

0.5

0.5

0 0

5

10

15 s

20

25

30

0

Figure 3. Brusselator (stable configuration): memory kernel snapshots. w1 (t) = φ1 (t)2 φ2 (t). The (estimated) contribution of φ1 (s) to the memory at t (yellow shaded area) is K1 (φ1 (s), t − s)ds for s ∈ [0, t]. (a) t = 4.28 s, (b) t = 8.55 s, (c) t = 17.83 s and (d) t = 30 s.

In figures 3 and 4, we provide snapshots of the estimated memory kernel K1 (φ1 (s), t − s) for (x1 , x2 ) = (1, 0). Figure 3a–d was computed for the stable configuration (A, B) = (1, 1.7), while figure 4a–d was computed for the unstable configuration (A, B) = (1, 3). In both configurations, the memory appears to have finite support.

(b) Viscous burgers equation The VBE in discrete Fourier space is considered: d ik uk + dt 2



up uq = −νk2 uk ,

k ∈ F ∪ G.

(4.2)

p+q=k (p,q)∈F∪G

The physical domain [0, 2π ] is discretized into 2N equidistant points {x0 , x1 , . . . , x2N−1 } with x0 = 0 and x2N−1 = 2π − x. The parameter ν is the viscosity. The relationship between the coefficients uk and the grid values u(xn ) is given by uk =

2N−1 

u(xn ) e2πik(n/N) ,

n=0

where F and G represent the sets of resolved and unresolved Fourier modes, respectively. Denote u = (uk )1≤k≤N the vector of Fourier coefficients and uˆ the same vector with all of G modes set to 0.

...................................................

2.0

2.5

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

contribution of f1(s) to memory w1 = memory estimated memory f1

2.5

0

14

3.0

3.0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

(b) contribution of f1(s) to memory w1 = memory estimated memory f1

0

2

4

6

8

10 12 14 16 18 20 s

(c)

20 18 16 14 12 10 8 6 4 2 0

15

0

2

4

6

8

10 12 14 16 18 20 s

0

2

4

6

8

10 12 14 16 18 20 s

(d)

0

2

4

6

8

20 18 16 14 12 10 8 6 4 2 0

10 12 14 16 18 20 s

Figure 4. Brusselator (LCO configuration): memory kernel snapshots. (a) t = 8 s, (b) t = 9 s, (c) t = 10 s and (d) t = 15 s.

Denote u0 = u(0). The initial conditions are given by [40] u(x, 0) =

m 

 2E(k) sin(kx + βk );

E(k) =

k=1

5−5/3

if 1 ≤ k ≤ 5

k−5/3

if k > 5

.

Note that the summation index stops at m, which is the number of resolved modes, in order to have fully resolved initial conditions. If the index went all the way down to N, the initial conditions would not fall within our scope. The subgrid terms for this problem are given by wk = −

ik 2





up uq −

p+q=k (p,q)∈F∪G



 up uq = P

p+q=k (p,q)∈F

t 0

ˆ Kk (u(s), t − s) ds,

The Liouville operator is ⎞

⎛ L≡

N−1 ⎜ k=1

⎜ ik ⎜− ⎝ 2

 p+q=k (p,q)∈F∪G

⎟ ∂ ⎟ u0p u0q − νk2 u0k ⎟ . ⎠ ∂u0k

Applying equation (3.6) to the VBE, we obtain ˆ ˆ t − s) ≈ R(u(s)) Kj (u(s),

¯ u(s)), ˆ + εR( ˆ Fj (u(s) t − s) ε

,

k ∈ F.

...................................................

20 18 16 14 12 10 8 6 4 2 0

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

20 18 16 14 12 10 8 6 4 2 0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

0.9

60

0.4

20

0.3

0.6

40 k

0.5

30

0.4

20

0.3

0.2 10

0.2 10

0.1 0

0.5

1.0 t

1.5

2.0

0

0.1 0

0.5

1.0 t

1.5

2.0

0

Figure 5. VBE: norm of (wk )1≤k≤64 represented as a continuous contour. (a) Integration over time of the estimated kernel and (b) using the full-order solution.

(a)

0.04 0.03 0.02 0.01 0 –0.01 –0.02 –0.03 –0.04 –0.05 –0.06

(b)

0.25 0.20

XF

0.15 0.10 0.05

0

0.5

1.0 t

1.5

0

2.0

0

0.5

1.0 t

1.5

2.0

Figure 6. VBE: numerical results (dotted line, using the full-order solution; solid line, integration of the estimated kernel) under the two metrics. (a) Mean (wk )1≤k≤Nc —real part and (b) contribution of the subgrid terms to the resolved energy decay. (Online version in colour.)

where F(u0 , s) is the r.h.s. term at time s of the corresponding pseudo orthogonal ODE: d Q ik u =− dt k 2

d Q ik uk = − dt 2 and uQ (0) = u0 .







Q uQ p uq

−

p+q=k (p,q)∈F∪G

 p+q=k (p,q)∈F∪G



 Q uQ p uq

p+q=k (p,q)∈F2 Q uQ p uq

−

 p+q=k (p,q)∈F2

,

k ∈ F,

 Q uQ p uq

− νk2 uQ k ,

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ k ∈ G⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

(4.3)

(i) Results The parameters are m = 64, N = 1048, tf = 2, ε = 10−8 , dt = 10−3 s and ν = 10−3 . A fourth-order Runge–Kutta integration scheme with 3/2 padding is used to solve equations (4.2) and (4.3).

...................................................

30

0.7

50

0.6 0.5

0.8

16

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

0.7

40

0.9

60

0.8

50

k

(b)

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

contribution of uˆ(s) to memory MZ memory exact memory

0

0

0.5

1.0 s

0.5

1.0 s

1.5

2.0

1.5

2.0

1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

(d) 1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

17

0

0.5

0

0.5

1.0 s

1.5

1.0 s

1.5

2.0

2.0

Figure 7. VBE: snapshots of the kernel (imaginary part) for cut-off mode 64. (a) t = 0.5, (b) t = 1.0, (c) t = 1.5 and (d) t = 2.0.

(a)

(b) 60 50 40

k

30 20 10 0

0.1

0.2

0.3 t

0.4

0.5

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0.21 0.20 0.19 0.18 t 0.17 0.16 0.15 0.14 0.13 10

20

30 k

40

50

60

Figure 8. VBE: (a) A first estimation of how rapidly the kernel Kj (ˆu, t) can be expected to decay with time t. (b) For each wavenumber k, an estimation of the time τ after which the decay profile goes below 1%.

A global picture of the results is given in figure 5 where contour plots of the norm of the subgrid terms (wk )1≤k≤64 computed using the full-order solution and reconstructed with our method are given. Two metrics are used to measure how well the subgrid terms were reproduced:

...................................................

(c)

(b)

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

(b)

18

10–1 t

0

0.05

0.10

0.15

0.20

0.25

10–2 102

t

m

Figure 9. VBE: impact of the ROM resolution m on the memory length. The higher the level of resolution, the faster the kernel decays. Panel (a) shows averaged cut-off decay profiles Km (ˆu, t) for different values of m. Panel (b) suggests a power-law relationship between m and the memory length.

0.40

0.05

0.35

0.04 0.03

0.30

0.02

0.25

0.01

t 0.20

0 –0.01

0.15

–0.02

0.10

–0.03

0.05 0

–0.04 1

2

3 x

4

5

6

–0.05

Figure 10. K-S: full-order solution in physical space.

— The first metric is a mean over all the resolved modes (figure 6a). The agreement between the estimation and the full-order results is excellent, although amplification errors are notable in the first third of the modes. This metric only quantifies whether the dominant memory terms have been well reproduced. For this fundamental problem in particular, the higher the wavenumber k, the bigger the magnitude of wk . — To obtain a more complete picture, the second metric quantifies the contribution of the subgrid terms to the resolved part of the energy decay rate. As        d d d 1 dE =− Re u∗k uk , uk u∗k + u∗k uk = − − dt F 2 dt dt dt k∈F

k∈F

the contribution from wk is defined as ΞF = − k∈F Re(u∗k wk ). The fact that the low wavenumber Fourier modes uk are more energetic than the high wavenumber modes compensates for the opposite trend in the subgrid terms wk , and makes ΞF a good indicator of whether the low wavenumber subgrid terms have been well reproduced

...................................................

m = 32 m = 64 m = 128 m = 192 m = 256 m = 384 m = 512

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a) 30

15 10 5 0

0.05

0.10

0.15

0.20 t

0.25

0.30

0.35

0.40

(b) 30 25 20 k 15 10 5 0

0.05

0.10

0.15

0.20 t

0.25

0.30

0.35

0.40

5000 4500 4000 3500 3000 2500 2000 1500 1000 500 0

Figure 11. K-S: norm of (wk )1≤k≤32 represented as a continuous contour. (a) Integration over time of the estimated kernel and (b) computed from the full-order solution.

(figure 6b). It appears that for this problem, the proposed method only enables a faithful reproduction of the high wavenumber subgrid terms. ˆ t − s) at four Figure 7 shows snapshots of the imaginary part of the memory kernel K64 (u(s), different time instants t. These figures suggest the finite support of the memory length and a decaying sinusoidal form of the kernel. To estimate the memory length, decay profiles for the kernels Kj have been built by averaging ˆ n ), s)|, s ∈ [0tn ], over all the trajectory points φ(t ˆ n ). This is shown and scaling the profiles |Kj (φ(t in figure 8a as a discrete contour plot. Figure 8b gathers the memory lengths associated with these m decay profiles. Figure 9a–b shows the effect of the size of the reduced-order model on the estimated memory length for w64 (k = 64 is the cut-off wavenumber). It is seen that as the number of resolved modes increases, the memory length decreases. For this specific problem, the relationship between the memory length and the resolution level amounts to a simple scaling law.

(c) Kuramoto–Sivashinsky equation The K-S equation [41,42] in discrete Fourier space is given by ik  d uk = (k2 − νk4 )uk − up uq , dt 2 p+q=k

k ∈ F ∪ G.

...................................................

20 k

19

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

25

5000 4500 4000 3500 3000 2500 2000 1500 1000 500 0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

500

20

400

real

100 0 –100 –200 –300 –400

0

0.05

0.10

0.15

0.20 t

0.25

0.30

0.35

0.40

Figure 12. K-S: numerical results under the first metric—mean (wk )1≤k≤Nc (real part) versus time. The dotted line is computed using the full-order solution. The solid line is obtained by integrating over time the estimated kernel. (Online version in colour.)

2.5

XF (×105)

2.0 1.5 1.0 0.5 0 0

0.05

0.10

0.15

0.20 t

0.25

0.30

0.35

0.40

Figure 13. K-S: numerical results under the second metric ΞF . Same legend as in figure 12. (Online version in colour.)

This equation has similarities with the VBE, but develops richer behaviour due to the presence of both second- and fourth-order spatial derivatives. The linear part k2 − νk4 is strictly positive √ for k < 1/ ν. The corresponding modes act as sources of energy for the entire system. Stability is ensured by the convection term that transfers the energy produced by the source at the large scales (low wavenumber) to the small scales (high wavenumber) where −νk4 dominates. The viscosity ν is chosen such that the system would display (figure 10) a state of persistent dynamical disorder. The initial conditions are the same as in the VBE. Since the K-S system is stiffer than the VBE, a fourth-order exponential time differencing scheme [43,44] is used instead of the Runge–Kutta scheme for robustness. The initial conditions and problem parameters are the same as for Burgers equation, except m = 32, tf = 0.4 s and dt = 10−4 s.

(i) Results The overall picture (contour of the subgrid terms) is given in figure 11. The mean of wk and ΞF are shown in figures 12 and 13, respectively. By both metrics, the method performs well. Figure 14 suggests a very short memory length compared with the Burgers equation. The estimated decay profiles and their corresponding memory lengths are provided in figure 15a,b.

...................................................

200

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

300

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

(a)

(b)

0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 s

(c)

1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

21

0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 s

(d)

0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 s

1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 s

Figure 14. K-S: snapshots of the kernel (real part) for mode 21. The estimated kernels decay very rapidly. The memory at time t is completely determined by the most recent kernel values Kj (ˆu(s), t − s), s ∈ [t − τ , t]. Note: The memory terms and the computed kernels have been scaled to fit together in the graphs. (a) t = 0.1, (b) t = 0.2, (c) t = 0.3 and (d) t = 0.4. (a) 30

(b)

25 20 k

15 10 5 0

0.005 0.010 0.015 0.020 t

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0.020 0.018 0.016 t

0.014 0.012 0.010 0.008 0.006 5

10

15

20

25

30

k

Figure 15. K-S: (a) A first estimation of the rate at which the kernel Kj (ˆu, s) decays with time. (b) For each wavenumber k, an estimation of the time τ after which the decay profile goes below 1% of its maximal amplitude. The presence of a finite support in time for the memory suggests that, in order to reconstruct the subgrid terms wj at time t from the M-Z perspective, one would only need to compute the ˆ t − s) for s ∈ [t − τ , t]. kernel Kj (φ(s),

...................................................

contribution of uˆ(s) to memory Re(w21) = memory estimated memory

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

1.0 0.8 0.6 0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1.0

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

5. Conclusion and perspectives

implemented, performed, and interpreted all the simulations in consultation with E.J.P. and K.D. A.G. wrote the manuscript. All authors gave final approval for publication. Competing interests. The authors have no competing interests. Funding. This research was funded by the AFOSR (grant no. FA9550-16-1-0309) under the project LES Modelling of Non-local effects using Statistical Coarse-graining (Tech. Monitor: Jean-Luc Cambier). Acknowledgements. A.G. thanks Dr Ghislain Haine for insightful conversations from which this work benefited.

References 1. Holmes P, Lumley JL, Berkooz G. 1996 Turbulence, coherent structures, dynamical systems and symmetry. Cambridge, UK: Cambridge University Press. 2. Sirovich L. 1987 Turbulence and the dynamics of coherent structures, part III: dynamics and scaling. Q. Appl. Math. 45, 583–590. (doi:10.1090/qam/910464) 3. Aubry N, Holmes P, Lumley J, Stone E. 1988 The dynamics of coherent structures in the wall region of a turbulent boundary layer. J. Fluid Mech. 192, 115–173. (doi:10.1017/S002211 2088001818) 4. Berkooz G, Holmes P, Lumley J. 1993 The proper orthogonal decomposition in the analysis of turbulent flows. Annu. Rev. Fluid. Mech. 25, 539–575. (doi:10.1146/annurev.fl. 25.010193.002543)

...................................................

Data accessibility. This article has no additional data. Authors’ contributions. A.G. developed the procedure starting from the foundations E.J.P. and K.D. laid. A.G.

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

The M-Z formalism provides a mathematically consistent framework for the construction of reduced-order models of dynamical systems. Using projection operators in phase space, the M-Z approach recasts a high-order nonlinear ODE system into generalized Langevin equation (GLE) for the resolved physics. In the GLE, the contribution from the unresolved physics is exactly represented in the form of a non-local memory integral that involves the past history of the resolved physics and a noise term that drops in the setting of the projector used. In this work, we presented an a priori procedure to estimate the memory kernel for the ˆ t − s) amounts to the truncation projector. In this context, computing the memory kernel Kj (φ(s), evaluation of the sensitivity with respect to the initial conditions of the solution to the orthogonal dynamics a direction imposed by the ROM trajectory. Instead of solving the orthogonal dynamics equation—which is a high-dimensional linear PDE—we propose to work with a simplified equation. This simplification is a result of the assumption that the semi-group of the orthogonal dynamics is a composition operator for the observable g(x) = QLx. The ersatz has the advantage of being tractable by evolving in time a characteristics equation that we call the pseudo orthogonal ODE. It has to be emphasized that the method and analysis proposed in this work is limited to the truncation projector only. The procedure is exact in the linear case where the kernel is known analytically. Based on numerical results for the Brusselator, the VBE and the K-S equation, we conjecture that the assumption made on the semi-group of the orthogonal dynamics is appropriate for a certain class of nonlinear problems. For Burgers, the dominating, near cut-off subgrid terms were well reproduced while the low-wavenumber were underestimated. For K-S, the method performed better overall. For the Burgers equation, a simple power-law scaling between the size of the resolved set and the memory length was observed. For K-S, the estimated memory length was shorter than for the Burgers equation. Further work will examine the validity of our ersatz approach for more complex problems. This work attempts to estimate the entire integrand of the memory term. The use of finiterank projections [14] to approximate the projection operator will allow for the decomposition of the memory integrand into a product of a time-dependent term and the state of the system. While such a decomposition may provide more insight into the modelling of the kernel, a large number of basis functions may be required, demanding the development of sparse approximation techniques.

22

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

23 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

5. Willcox K, Peraire J. 2002 Balanced model reduction via the proper orthogonal decomposition. AIAA J. 40, 2323–2330. (doi:10.2514/2.1570) 6. Rowley CW. 2005 Model reduction for fluids using balanced proper orthogonal decomposition. Int. J. Bifurcation Chaos 15, 997–1013. (doi:10.1142/S0218127405012429) 7. Gugercin S, Antoulas AC. 2004 A survey of model reduction by balanced truncation and some new results. Int. J. Control 77, 748–766. (doi:10.1080/00207170410001713448) 8. Moore G. 1981 Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE. Trans. Automat. Contr. 26, 17–32. (doi:10.1109/TAC.1981.1102574) 9. Veroy K, Patera AT. 2005 Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-bases a posteriori error bounds. J. Numer. Methods Fluids 47, 773–788. (doi:10.1002/fld.867) 10. Rozza G. 2011 Reduced basis approximation and error bounds for potential flows in parametrized geometries. Commun. Comput. Phys. 9, 1–48. (doi:10.4208/cicp.100310.260710a) 11. Beattie CA, Gugercin S. 2005 Krylov-based model reduction of second-order systems with proportional damping. In Proc. of the 44th IEEE Conf. on Decision and Control, and the European Control Conference 2005, Sevilla, Spain, pp. 2278–2283. Piscataway, NJ: IEEE. (doi:10.1109/CDC.2005.1582501) 12. Smagorinsky J. 1963 General circulation experiments with the primitive equations. Mon. Weather Rev. 91, 99–164. (doi:10.1175/1520-0493(1963)0912.3.CO;2) 13. Deardorff J. 1970 A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. J. Fluid. Mech. 41, 453–480. (doi:10.1017/S0022112070000691) 14. Chorin A, Hald O, Kupferman R. 2002 Optimal prediction with memory. Phys. D 166, 239–257. (doi:10.1016/S0167-2789(02)00446-3) 15. Chorin A, Stinis P. 2005 Problem reduction, renormalization, and memory. Commun. Appl. Math. Comput. Sci. 1, 1–27. (doi:10.2140/camcos.2006.1.1) 16. Chorin AJ, Hald O. 2009 Stochastic tools for mathematics and science. Berlin, Germany: Springer. 17. Chorin AJ, Hald O, Kupferman R. 2000 Optimal prediction and the Mori-Zwanzig representation of irreversible processes. Proc. Natl Acad. Sci. USA 97, 2968–2973. (doi:10.1073/ pnas.97.7.2968) 18. Mori H. 1965 Transport, collective motion and Brownian motion. Prog. Theor. Phys. 33, 423–450. (doi:10.1143/PTP.33.423) 19. Zwanzig R. 1973 Nonlinear generalized Langevin equations. J. Stat. Phys. 9, 215–220. (doi:10.1007/BF01008729) 20. Parish E, Duraisamy K. 2016 Reduced order modeling of turbulent flows using statistical coarse-graining. In 53rd AIAA Aerospace Sciences Meeting, 13–17 June 2016, Washington, DC, USA. Reston, VA: AIAA. (doi:10.2514/6.2016-3640) 21. Parish EJ, Duraisamy K. 2017 Non-Markovian closure models for large eddy simulations using the Mori-Zwanzig formalism. Phy. Rev. Fluids 2, 014604. (doi:10.1103/Phys RevFluids.2.014604) 22. Zwanzig R. 2001 Nonequilibrium statistical mechanics. New York, NY: Oxford University Press. 23. Bernstein D. 2007 Optimal prediction of Burgers’s equation. Multis. Model. Simul. 6, 27–52. (doi:10.1137/060651720) 24. Stinis P. 2007 Higher order Mori-Zwanzig models for the Euler equations. SIAM Multiscale Model. Simul. 6, 741–760. (doi:10.1137/06066504X) 25. Venturi D, Cho H, Karniadakis GE. 2016 Mori-Zwanzig approach to uncertainty quantification. In Handbook of uncertainty quantification, pp. 1037–1073. Berlin, Germany: Springer. 26. Venturi D, Karniadakis GE. 2014 Convolutionless Nakajima–Zwanzig equations for stochastic analysis in nonlinear dynamical systems. Proc. R. Soc. A 470, 20130754. (doi:10.1098/rspa. 2013.0754) 27. Cho H, Venturi D, Karniadakis GE. 2014 Statistical analysis and simulation of random shocks in stochastic Burgers equation. Proc. R. Soc. A 470, 20140080. (doi:10.1098/rspa.2014.0080) 28. Stinis P. 2015 Renormalized Mori-Zwanzig reduced models for systems without scale separation. Proc. R. Soc. A 471, 20140446. (doi:10.1098/rspa.2014.0446) 29. Givon D, Hald OH, Kupferman R. 2005 Existence proof for orthogonal dynamics and the Mori-Zwanzig formalism. Israel J. Math. 145, 221–241. (doi:10.1007/BF02786691) 30. Darve E, Solomon J, Amirali K. 2009 Computing generalized Langevin equations and generalized Fokker–Planck equations. Proc. Natl Acad. Sci. USA 106, 10 884–10 889. (doi:10.1073/pnas.0902633106)

Downloaded from http://rspa.royalsocietypublishing.org/ on September 30, 2017

24 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 473: 20170385

31. Prigogine I. 1980 From being to becoming: time and complexity in the physical sciences. New York, NY: W.H. Freeman and Company. 32. Hynes JT, Deutch JM. 1975 Non-equilibrium problems—projection operator techniques. In Physical chemistry—an advanced treatise, vol. XIB, ch. 11, pp. 739–836. New York, NY: Academic Press. 33. Mezic I. 2005 Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dyn. 41, 309–325. (doi:10.1007/s11071-005-2824-x) 34. MarkoBudii M, Mohr R, Mezic I. 2012 Applied Koopmanism. Chaos 22, 047510. (doi:10.1063/1.4772195) 35. Evans D, Morriss G. 1990 Statistical mechanics of nonequilibrium liquids. London, UK: Academic Press. 36. Givon D, Kupferman R, Stuart A. 2004 Extracting macroscopic dynamics: model problems and algorithms. Nonlinearity 17, R55–R127. (doi:10.1088/0951-7715/17/6/R01) 37. Pazy A. 1983 Semigroups of linear operators and applications to partial differential equations. Applied Mathematical Sciences, vol. 44. New York, NY: Springer. 38. Stinis P. 2006 Stochastic optimal prediction for the Kuramoto–Sivashinsky equation. Multiscale Model. Simul. 2, 580–612. (doi:10.1137/030600424) 39. Parish EJ, Duraisamy K. 2017 A dynamic subgrid scale model for large eddy simulations based on the Mori-Zwanzig formalism. J. Comp. Phys. 349, 154–175. (doi:10.1016/j.jcp. 2017.07.053) 40. Li Y, Wang Z. 2015 A priori and a posteriori evaluations of subgrid stress models with the Burgers’ equation. In 53rd AIAA Aerospace Sciences Meeting, Kissimmee, FL, USA. Reston, VA: AIAA. (doi:10.2514/6.2015-1283) 41. Hyman JM, Nicolaenko B. 1986 The Kuramoto-Sivanshinsky equation: a bridge between partial differential equations and dynamical system. Phys. D 18, 113–126. (doi:10.1016/01672789(86)90166-1) 42. Nicolaenko B, Scheurer B, Temam B. 1985 Some global dynamical properties of the Kuramoto-Sivashinsky equation: nonlinear stability and attractors. Phys. D 16, 155–183. (doi:10.1016/0167-2789(85)90056-9) 43. Cox SM, Matthews PC. 2002 Exponential time differencing for stiff systems. J. Comput. Phys. 176, 430–455. (doi:10.1006/jcph.2002.6995) 44. Kassam A, Trefethen L. 2005 Fourth-order time-stepping for stiff PDEs. J. Sci. Comput. 26, 1214–1233. (doi:10.1137/S1064827502410633)