A Unified Framework for Multiscale Modeling using the Mori-Zwanzig Formalism and the Variational Multiscale Method
Eric J. Parish Karthik Duraisamy
math.NA











parish@umich.edu
kdur@umich.edu
Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI 48109, USA














We describe a paradigm for multiscale modeling that combines the Mori-Zwanzig (MZ) formalism of Statistical Mechanics with the  Variational Multiscale (VMS) method. The MZ-VMS approach leverages both VMS scale-separation projectors as well as phase-space projectors to provide a systematic modeling approach that is applicable to  non-linear partial differential equations. Spectral as well as continuous and discontinuous finite element methods are considered. The framework leads to a formally closed equation in which the effect of the unresolved scales on the resolved scales is non-local in time and appears as a convolution or memory integral. The resulting non-Markovian system is used as a starting point for model development.  We discover that unresolved scales lead to memory effects that are driven by an orthogonal projection of the coarse-scale residual and inter-element jumps. It is further shown that an MZ-based finite memory model is a variant of the well-known adjoint-stabilization method. For hyperbolic equations, this stabilization is shown to have the form of an artificial viscosity term. We further establish connections between the memory kernel and approximate Riemann solvers. It is demonstrated that, in the case of one-dimensional linear advection, the assumption of a finite memory and a  linear quadrature leads to a closure term that is formally equivalent to an upwind flux correction.










§ INTRODUCTION



Advances in computational hardware and numerical algorithms have allowed for the pursuit of numerical simulations of increasingly complex systems. Direct computations of practical multiscale problems, however, have remained prohibitively expensive due to the presence of a wide range of length and time scales. Incorrectly accounting for unresolved scales can lead to inaccurate solutions, and in many cases, unstable numerical schemes. The pursuit of efficient and robust numerical methods is thus a pacing item in computational physics. The Variational Multiscale (VMS) approach of Hughes <cit.> is an elegant mathematical procedure to derive stabilization/closure schemes for numerical simulations of multiscale problems. The VMS procedure is centered around a sum decomposition of the solution u in terms of resolved scales ũ and unresolved or fine-scales u^'. The key challenge in VMS is to obtain a representation for u^' in terms of ũ. This is typically done by virtue of approximated Green's functions. Given the approximation to the fine-scales, the impact of the fine-scales on the coarse-scales can then be numerically computed. In Hughes' pioneering work, it was demonstrated that classic stabilization techniques, such adjoint stabilization and streamline upwind Petrov-Galerkin, can be derived through specific approximations to the Green's function and manipulations of the multiscale equations. The Variational Multiscale Method has since gained significant attention, and reviews may be found in <cit.>. The VMS procedure by itself, however, does not eliminate the closure problem. Representation of the fine-scale state in terms of the coarse-scale state, by virtue of a Green's function or otherwise, is required. For the non-linear multiscale problems commonly encountered in science and engineering, the development of such relations has proved challenging and the VMS procedure is typically extended in an adhoc 
fashion. 

A less publicized but conceptually similar approach to VMS is the optimal prediction framework developed by Chorin and co-workers <cit.>. Chorin's optimal prediction framework, which is a reformulation of the Mori-Zwanzig (MZ) formalism of statistical mechanics <cit.>, can be viewed as a  model order reduction strategy for ordinary differential equations. The framework centers around the use of projection operators that separate the phase space of an ordinary differential equation into resolved and unresolved manifolds. Using this decomposition, a high-dimensional non-linear Markovian dynamical system can be recast into an equivalent, lower-dimensional non-Markovian system. This formulation is mathematically exact. In the lower-dimensional system, which is commonly referred to as the generalized Langevin equation (GLE), the effect of the unresolved scales on the resolved scales is non-local in time and appears as a convolution integral. This term is typically referred to as memory. The memory term emerging from the MZ formalism is analogous to the Green's function emerging from the VMS procedure. While the memory integral is intractable in general non-linear problems, it serves as a mathematically rigorous starting point  for the construction of closure models. 

In addition to the pioneering work of Chorin and collaborators, a significant body of research on Mori-Zwanzig approaches exists in the literature.  The MZ formalism has been examined extensively in fields such as molecular dynamics <cit.> and uncertainty quantification <cit.>. In the context of solution of partial differential equations, the main body of work has examined approximations to the memory term and their application to the semi-discrete systems emerging from Fourier-Galerkin spatial discretizations.  Stinis and coworkers <cit.> have developed several models for approximating the memory, including finite memory and renormalized models, and examined their performance on the Burgers equation and the Euler equations. In the same spirit, the work of Zhu, Dominy, and Venturi <cit.> examines additional approximations to the memory as well as error estimates for various MZ models. The present authors have examined the performance of the t-model and several of Stinis' models in the context of the Navier-Stokes equations <cit.>. We have further investigated the development of a parameter-free MZ-based model <cit.> that combines ideas from the MZ community with the Germano identity.

The Mori-Zwanzig procedure, however, has received minimal attention as a practical tool for closure modeling in the context of the numerical solution of partial differential equations. The vast majority of the work has been undertaken on either idealized systems or problems that pertain to the statistical mechanics community. Further, MZ-based approaches have not  been applied to numerical simulations of partial differential equations using practical discretization techniques such as finite elements. In fluid dynamics and Large Eddy Simulation, MZ-related modeling has been pursued on the semi-discrete systems emerging from spatial discretizations using Fourier-Galerkin spectral methods. MZ-based models have been investigated in the context of Fourier-Galerkin simulations of Burgers equation <cit.>, the Euler equations <cit.>, and the Navier–Stokes equations <cit.>. This body of work has demonstrated the potential of using the MZ procedure for model development. The Fourier-Galerkin spectral method, however, is limited to canonical problems and periodic domains. The extension of the Mori-Zwanzig procedure to more general numerical methods has remained unclear.
This is largely due to the challenges in formulating the MZ framework in a systematic manner for the ordinary differential equations that emerge from the numerical discretization of partial differential equations. These challenges stem from the fact that the key concept in MZ is a tractable decomposition of the discrete unknowns into a coarse-scale resolved set and a fine-scale unresolved set. While this separation is clear in the context of Fourier-Galerkin methods, it is less straightforward for more general numerical schemes. Another road-block to the further development of MZ methods is the lack of clarity regarding connections to other popular techniques. This work addresses the aforementioned issues. 

 

The main goal of this work is to describe a closure modeling approach that combines the Mori-Zwanzig formalism and the Variational Multiscale Method. We refer to this framework as MZ-VMS. We also establish novel connections with several existing stabilization methods including orthogonal subscales, adjoint stabilization, and upwind fluxes. Spectral methods as well as continuous/discontinuous finite element methods are considered. Similar to the VMS procedure, the framework relies on an a priori decomposition of the solution space. The resulting equations for ũ and u^' are then considered at the semi-discrete level. The Mori-Zwanzig formalism is used to eliminate the dependence of the coarse-scales on the fine-scales and serves as a starting point for model development. Figure <ref> provides an overview for how the MZ-VMS framework can be used for model development. 




The organization of this paper will be as follows: In Section <ref> we will provide a review of the variational multiscale approach. In Section <ref>, we will outline the development of the Mori-Zwanzig formalism for the VMS equations in the spectral method case. Section <ref> presents new insights into the memory term and will discuss several MZ-based models. Section <ref> will extend the derivation to the finite element method. Section <ref> will discuss connections between MZ-based modeling approaches, stabilization methods, artificial viscosity, and approximate Riemann solvers. Conclusions and perspectives will be provided in Section <ref>.



§ THE VARIATIONAL MULTISCALE METHOD

We start by reviewing the variational multiscale method. We refer the interested reader to <cit.> for a more complete description.
The variational multiscale method is typically introduced by considering a homogeneous Dirichlet problem,

    L u = f     in    Ω,


    u = 0     on    Γ,

where Ω⊂ R^d. The variational counterpart to Eq. <ref> is now formulated. We restrict ourselves to the Galerkin setting for simplicity. Let V≡ H_0^1 ( Ω ) represent the space of the solution and weighting functions. The variational problem is defined as follows: find u ∈V such that ∀ w ∈V,

    (w, Lu)= (w,f),

where (·, ·) is the L^2 inner product.
The variational multiscale method utilizes a decomposition of the solution space into a coarse-scale resolved space Ṽ⊂V and a fine-scale unresolved space V^'⊂V. In VMS, the solution space is expressed as a sum decomposition,

    V = Ṽ⊕V^'.


Let Π̃ be the linear projector onto the coarse-scale space, Π̃ : V→Ṽ. Various choices exist for the projector Π̃ <cit.>, and here we exclusively consider Π̃ to be the L^2 projector,

    (w̃,Π̃ u) = (w̃,u)     ∀w̃∈Ṽ, u ∈V.

The fine-scale space, V^', becomes the orthogonal complement of Ṽ in V such that,

    (w̃,Π^'u) = 0    ∀w̃∈Ṽ , u ∈V,

where Π^' = 𝐈 - Π̃.
With this decomposition, the solution can be represented as

    u = (Π̃ + Π^')u =  ũ + u^',

and the same for w. For purposes of clarity, it is assumed that ũ and u^' are homogeneous on Γ.
By virtue of the linear independence of the fine and coarse trial spaces, Eq. <ref> separates into two sub-problems,

    (w̃,Lũ) +  (w̃,Lu^') = (w̃,f),     ∀w̃∈Ṽ,


    (w^',Lũ) +  (w^',Lu^') = (w^',f)     ∀w^'∈V^'.

Equations <ref> and <ref> may be equivalently written as

    (w̃,Lũ) + (L^*w̃,u^') = (w̃,f),     ∀w̃∈Ṽ,


    (w^', Lũ) +  (w^',Lu^') = (w^',f)     ∀w^'∈V^',

where we have used the integration-by-parts formula,

    ( w̃ , L u^') = (L^* w̃, u^'),     ∀w̃∈Ṽ, u^'∈V^'.







The operator L^* is the adjoint of L. The objective is to analytically solve for the fine-scale solution, u^', and inject it into the coarse-scale equation. This can be achieved by solving the Green's function problem corresponding to the Euler-Lagrange equations of the fine-scale equation. The final expression for the fine-scale state is 

    u^'(y) = - ∫_Ωg^'(x,y) (L ũ - f )(x) dΩ_x,

where g^' is the fine-scale Green's function satisfying

    Π^'L^* g^'(x,y) = Π^'δ (x-y)     ∀  x ∈Ω,


    g^'(x,y) = 0      on     Γ.








       
             [t]0.95
The  equation for the resolved, “coarse" scales can be obtained as

    (w̃,Lũ)  + ( L^* w̃,  M^'( L ũ -f)  ) = (w̃,f),

where

    ( L^* w̃,  M^'( L ũ -f)  ) = - ∫_Ω  (L^* w̃)(y) ∫_Ωg^'(x,y)(Lũ - f)(x) dΩ_x dΩ_y.


             
          
       


In practice, the fine-scale Green's function is approximated. A common model assumes the fine-scale Green's function to be

    g^'(x,y ) ≈τ(x) δ (x - y),

where τ is a stabilization parameter. This approximation leads to the following expression for the fine-scale solution,

    u^'(y) ≈ -τ(L ũ - f ).

The VMS coarse-scale equation then appears as follows

    (w̃,L ũ)   - (  L^* w̃, τ(Lũ - f) ) = (w̃,f).

The approximation given by Eq. <ref> is an adjoint-type stabilization technique and was proposed <cit.> before the VMS procedure. One of the significant accomplishments of the VMS framework is that it elucidates the mathematical origins of such stabilization techniques.

Remarks


  * In practice the Green's function is  approximated. The type of approximation to the Green's function defines various VMS-based models.

  * The extension of the Green's function to non-linear problems is not straight-forward.

  * The VMS procedure does not eliminate the modeling problem, and some method must be devised to approximate u^'. This issue is amplified in non-linear and unsteady problems. It is for this purpose that we develop the Mori-Zwanzig formalism for VMS.





 §.§ Hierarchical Refinement and Static Condensation

The fine-scale Green's function can be better understood through hierarchical refinement. The expression of the state variable in the trial space is

    u (x,t)= ∑_j = 0^∞ w_j(x) a_j(t),

where a_j are the modal coefficients for with we seek to solve. Compactly, 

    u(x,t) = 𝐰^T(x) 𝐚(t).

We consider the hierarchical decomposition of the state variable,

    u (x,t)= ũ(x,t) + u^'(x,t) = ∑_j = 0^N w_j(x) a_j(t) + ∑_j = N^∞ w_j(x) a_j(t),

where the first N basis functions are in Ṽ and the remaining basis functions are in V^'. Compactly,

    u(x,t) =  ũ(x,t) + u^'(x,t)  = 𝐰̃^T ã + ŵ^T 𝐚^'.

Hierarchical decompositions may be obtained, for example, with Legendre polynomials or Fourier series. Note that such a decomposition constrains ũ∈Ṽ and u^'∈V^'.
With this definition, one can analytically solve the fine-scale state u^' in terms of the coarse-scales,

    𝐮̂ = 𝐊̂^-1(𝐰^',L 𝐮̃ - f),

where 

    𝐊̂ = (𝐰^',L 𝐰^'^T).

The coarse-scale equation becomes

    (w̃, L ũ -f)  +  (w̃, L u^' ) = 0     ∀w̃∈Ṽ.





§ FORMULATION OF THE MORI-ZWANZIG FORMALISM FOR VMS (SPECTRAL METHODS)

We now formulate the Mori-Zwanzig formalism for the globally smooth case within the context of the VMS method. The MZ procedure works with the semi-discrete system that emerges after a spatial discretization. Similar to VMS, the discrete unknowns are partitioned into a coarse (resolved) and fine (unresolved) scales. The dependence of the coarse-scale equation on the fine-scales is removed by obtaining an analytic solution to the fine-scale solution. In the linear case, this can be achieved through integrating factors, while for the non-linear case MZ makes use of Liouville equations and Duhamel's principle. The formal removal of the fine-scale state does not lead to a reduction in computational complexity, but instead serves as a starting point for model development. We first illustrate the philosophy in the linear case before proceeding to the non-linear problem.


 §.§ Linear initial value problem: integrating factors

 Consider the linear initial value problem,

    ∂ u/∂ t + Lu = f       x ∈Ω, t ∈ (0,T),

subject to the boundary and initial conditions,

    u(x,t)    = 0       x ∈Γ,t ∈ (0,T), 
    
         u(x,0)    = u_0       x ∈Ω.

Note that L is again a differential operator. Letting V≡ H_0^1(Ω) again denote the test and trial space, the Galerkin formulation at the semi-discrete level leads to the weighted residual problem,

    (w, u_t ) + (w , Lu -f) = 0     ∀ w ∈V.

Separating the test/trial space via the sum decomposition leads to two sub-problems,

    (w̃ , ũ_t ) + (w̃, L ũ  )  + (w̃ , u^'_t ) + (w̃, Lu^'  ) = (w̃, f)      ∀w̃∈Ṽ,


    (w^' , ũ_t ) + (w^', L ũ )   + (w^' , u^'_t ) + (w^', L u^'  ) = (w^', f)      ∀w^'∈V^'.

We assume u_0  ∈Ṽ such that

    ũ(x,0) = ũ_0,


    u^'(x,0) = 0.

Consider the set of discrete equations arising from Eqns. <ref> and <ref>. Let the vector of trial/test functions be denoted by 𝐰. The expression of the state variable in the trial space is

    u (x,t)= ∑_j = 0^∞ w_j(x) a_j(t),

where a_j are the modal coefficients to be solved for. Compactly, 

    u(x,t) = 𝐰^T(x) 𝐚(t).

The key step in the formulation of the Mori-Zwanzig framework is the restriction that the coarse-scale space and fine-scale space are L^2 orthogonal such that

    (w̃,w^') = 0     ∀w̃∈Ṽ,w^'∈V^'.

Constraining the coarse-scale and fine-scale spaces to be L^2 orthogonal leads to hierarchical discretizations in the sense that the fine-scale stencil does not impact the coarse-scale stencil. With the definition of the projector in Eq. <ref>, the fine-scale space is implicitly L^2 orthogonal such that Eq. <ref> is satisfied.  

Equations <ref> and <ref> correspond to an infinite dimensional set of ordinary differential equations. The ODE set arising from Eq. <ref> is finite dimensional, while that arising from Eq. <ref> is infinite dimensional. With the orthogonality constraint, these ODE systems are written as

    (𝐰̃ , ũ_t ) + (𝐰̃, Lũ  )  + (𝐰̃,Lu^' ) = (𝐰̃, f),


    (𝐰^' , u^'_t ) + (𝐰^', Lũ  )  + (𝐰^', Lu^' ) = (𝐰^', f).



Similar to the VMS approach, the objective is to analytically solve for the fine-scale state, u^'. This can be achieved through the use of integrating factors. For simplicity, we restrict ourselves to the case of orthonormal basis functions, although no such simplification is necessary. One obtains

    u^'(t) = -∫_0^t  𝐰^'^T e^ -s(𝐰^', L 𝐰^'^T) (𝐰^',L ũ(t-s) - f) ds ,

where e^ s(𝐰^', L 𝐰^'^T) is the matrix exponential. Injecting the fine-scale solution into the coarse-scale equation yields,

    (𝐰̃ , ũ_t ) + (𝐰̃, L ũ)  - (𝐰̃,  ∫_0^t L 𝐰^'^T e^ -s(𝐰^', L 𝐰^'^T) (𝐰^',L ũ(t-s) - f) ds  ) = (𝐰̃, f ).



       
             [t]0.95
The convolution integral is referred to as memory.
Equation <ref> can be written in a more transparent form by recognizing that the memory contains terms similar to a projection. Define

    Π^'^Q(t) g = 𝐰^'^T e^ -t(𝐰^', L 𝐰^'^T) (𝐰^',g).

With this notation,

    (𝐰̃ , ũ_t )  + (𝐰̃, L ũ)  - (𝐰̃,  ∫_0^t L  Π^'^Q(s) [ L ũ(t-s) - f ] ds  ) = (𝐰̃, f ).



             
          
       



Remarks:
 

  * Equation <ref> is an exact, closed equation for the coarse-scales.

  * The formal removal of the fine-scales leads to a coarse-scale problem that is non-local in time.

  * Similar to the VMS approach, the memory is driven by the residual of the coarse-scale equation. Additionally, the projection associated with Eq. <ref> constrains the residual to be in the space of the fine-scales. This is similar to the VMS approach, where the Green's function problem is defined on the fine-scales. 

  * Even though Equation <ref> is written for the coarse scales, no reduction in computational cost has been achieved as we have simply transformed the system into an integro-differential equation.

  * Integrating factors are an instance of Duhamel's principle for ordinary differential equations. Duhamel's principle will be used extensively in what follows.

  * Equation <ref> should be viewed as a starting point for model development in complex problems.





 §.§ Extension to non-linear problems: The Mori-Zwanzig Formalism

In the linear case, we showed that it is possible to obtain an analytic expression for the fine-scale state through the use of integrating factors. The method of integrating factors is, unfortunately, not generalizable to non-linear problems. The MZ procedure addresses non-linearity by casting a non-linear ordinary differential equation as an equivalent linear partial differential equation in phase space. The fine-scale state is then analytically solved in terms of evolution operators in phase space. This framework is detailed in what follows.

We consider the initial-value problem

    ∂ u/∂ t + R(u) = f           x ∈Ω, t ∈ (0,T),

where R is a non-linear differential operator. The governing equations are subject to  boundary and initial conditions

    u(x,t)    = 0       x ∈Γ,t ∈ (0,T), 
    
         u(x,0)    = u_0       x ∈Ω.

 With the trial and test space V≡ H_0^1, the discrete multiscale equations for the Galerkin method can be written as

    (𝐰̃, ũ_t) + (𝐰̃, R( ũ ) )  + (𝐰̃, R(u ) -  R( ũ )) = (𝐰̃, f),


    (𝐰^', u^'_t) + (𝐰^', R( ũ ) )  + (𝐰^', R(u ) -  R( ũ )) = (𝐰^', f),

with 𝐚(t=0) = 𝐚_0. We again consider orthonormal basis functions for simplicity, although this is not required (eg. see Section <ref>). We denote the Hilbert space in which 𝐚_0 resides as H. 
The objective is again to analytically solve for the fine-scale state. 


  §.§.§ Transformation to Phase Space and the Liouville Equation

The Green's function used in the traditional VMS approach as well as the integrating factor approach described in the previous section relies on the principle of superposition and is limited to linear systems. The MZ procedure addresses non-linearity by casting the original semi-discrete Galerkin system as a partial differential equation that exists in phase space. It is worth emphasizing that, although the Mori-Zwanzig formalism has its roots in the work of Mori and Zwanzig, the work of Chorin and collaborators <cit.> is a significant revamp of the formalism and extends the concept to general systems of ordinary differential equations. The following discussion is inspired by  <cit.> and we refer the reader to both <cit.> and <cit.> for clarification on any of the following points.

The non-linear ODE system defined by Eqns. <ref> and <ref> can be cast as the following linear partial differential equation <cit.>,

    ∂ v/∂ t = Lv,

where L is the Liouville operator and is given by,

    L = ∑_j=0^∞ (w_j, f - R(u_0))∂/∂ a_0j.

Equation <ref> is referred to as the Liouville equation and is an exact statement of the original dynamics.
The solution to Eq. <ref> is given by <cit.>,

    v(𝐚_0,t) = g(𝐚(𝐚_0,t)).

Semi-group notation is now used, in which the solution to the Liouville equation may be expressed as,

    v(𝐚_0,t) = e^t L  g(𝐚(𝐚_0,0)) .

The proper interpretation of Eq. <ref> is that  v(𝐚_0,t), is given by the solution of the partial differential equation defined by the evolution operator e^t L (in this case Eq. <ref>), with initial conditions v(𝐚_0,0) = g(𝐚(𝐚_0,0)).
By virtue of Eq. <ref>, it is seen that the evolution operator e^t L is a Koopman operator <cit.>, and we can alternatively write the solution as

    v(𝐚_0,t) =  g( e^t L𝐚(𝐚_0,0)).

The implications of e^t L are significant. It demonstrates that, given the trajectories 𝐚(𝐚_0,t), the solution v is known for any observable g. This is analogous to the method of characteristics for hyperbolic systems. Noting that L and e^t L commute, Eq. <ref> may be written in the semi-group notation as,

    ∂ v/∂ t = e^t LL v(𝐚_0,0).

A set of equations for the resolved modes can be obtained by taking g(𝐚_0,t) = 𝐚̃_0,

    ∂/∂ t  e^t Lã_0= e^t LLã_0.



  §.§.§ Projection Operators and the Liouville Equation

We proceed by decomposing the space of initial conditions, H, into a resolved and unresolved subspace,

    H = H̃⊕H^'.

The associated projection operators are defined as 𝒫: H→H̃ and 𝒬 = I - 𝒫. We consider a projection operator that is appropriate for the deterministic initial conditions considered here (Eqns. <ref> and <ref>),

    Pf(𝐚̃_0,𝐚^'_0) = ∫_H f(𝐚̃_0,𝐚^'_0) δ(𝐚^'_0) d 𝐚^'_0,

which leads to

    Pf(𝐚̃_0,𝐚^'_0) = f(𝐚̃_0,0).



Other projections, such as conditional expectations are possible  <cit.>, but will not be pursued in the present work. It is important to emphasize that the projectors P and Q operate on H and are fundamentally different from the L^2 projectors Π̃ and Π^' (see Figure <ref>). With the projection operators, the Liouville equation can be split as,

    ∂/∂ t  e^t Lã_0 = e^t LPLã_0 +  e^t LQLã_0.

The objective now is to remove the dependence of the right hand side of Eq. <ref> on the unresolved scales, 𝐚^' (i.e. QL𝐚̃). To demonstrate how this may be achieved, consider the partial differential operator governed by the semigroup e^t L written as,

    ∂/∂ t   -    L = 0.

We will refer to Eq. <ref> as the homogeneous problem. Consider now the inhomogeneous problem with forcing PL,

    ∂/∂ t   -    L = -PL.

Making use of the identity I = P  + Q, the inhomogeneous problem can be written as

    ∂/∂ t   - QL  =  0.

 Eq. <ref> is referred to in the literature as the orthogonal dynamics operator, and can be conceptualized as a Liouville operator with forcing. The evolution operator given by the orthogonal dynamics is e^t QL. Here, we can leverage the  linearity of the partial differential operators and make use of superposition. The solution to the orthogonal dynamics equation can be expressed in terms of solutions to the homogeneous Liouville equation through Duhamel's principle (in operator form),

    e^t ℒ = e^t 𝒬ℒ + ∫_0^t e^(t - s)ℒ𝒫ℒ e^s 𝒬ℒ ds.

Inserting Eq. <ref> into Eq. <ref>, the generalized Langevin equation is obtained,

    ∂/∂ t  e^t Lã_0=   e^tLPLã_0_Markovian +   e^tQLQLã_0_Noise + 
     ∫_0^t e^(t - s)ℒ𝒫ℒ e^s 𝒬ℒQLã_0 ds_Memory.

The system described in Eq. <ref> is precise and not an approximation to the original ODE system. For notational purposes, define

    F_j(𝐚_0,t) = e^tQLQL𝐚_0,      K_j(𝐚_0,t) = PLF_j(𝐚_0,t).

We refer to K(𝐚_0,t) as the memory kernel. It can be shown that solutions to the orthogonal dynamics equation are in the null space of P, meaning PF_j(𝐚_0,t) = 0. 
By the definition of the initial conditions (Eq. <ref>), the noise-term is zero and we obtain,

    ∂/∂ t    e^t Lã_0 = e^tLPLã_0+   ∫_0^t e^(t - s)ℒ𝒫ℒ e^s 𝒬ℒQLã_0 ds.



       
             [t]0.95
Equation <ref> can be written in a more transparent form,

    (𝐰̃,ũ_t)  + (𝐰̃, R(ũ) )  -    ∫_0^t K(ã(t-s),s) ds = (𝐰̃, f ),

where K_j(𝐚_0,t)  = PL e^tQLQL𝐚_0. Note that the time derivative is represented as a partial derivative due to the Liouville operators embedded in the memory. 

             
          
       


Remarks


  * Equation <ref> is precisely a Galerkin discretization of Eq. <ref> with the addition of a memory term originating from scale separation. 

  * Equation <ref> is a non-local closed equation for the coarse-scales.

  * Evaluation of the memory term is not tractable as it involves the solution of the evolution operator, e^t QL. This is referred to as the orthogonal dynamics and is discussed in the following section. Instead, Eq. <ref> is viewed as a starting point for the construction of closure models.

  * The formulation for the smooth case is limited to spectral methods on canonical domains. In Section <ref>, we extend the discussion to the case of finite elements.





§ THE MEMORY TERM: INSIGHT AND MODELING

The MZ-VMS procedure itself does not provide a reduction in computational complexity as it has replaced the fine-scale state with a memory term. This memory term relies on solutions to the orthogonal dynamics equation, which  is intractable in the general case. The MZ-VMS procedure instead provides an exact representation of the fine-scale state in terms of the coarse-scales. This is used as a starting point for model development. In this section, we expand our discussion of the orthogonal dynamics equation and discuss several modeling strategies. In particular we discuss the value of the memory kernel at s=0 and demonstrate that many existing MZ models are residual-based closures.




 §.§ The orthogonal dynamics

In the Variational Multiscale method, the fine-scale state is parameterized in terms of the coarse-scale state by virtue of a fine-scale Green's function. The Mori-Zwanzig procedure instead uses Duhamel's principle to relate the solution of the orthogonal dynamics equation to the Liouville equation. This allows for the elimination of fine-scales. The evolution operator of the orthogonal dynamics is given by e^t QL. It is important to recognize that, while the evolution operator e^t L is a Koopman operator, no such result exists for e^t QL in the general non-linear case. As a consequence, evaluating terms evolved by e^t QL requires one to directly solve the orthogonal dynamics equation. This is, in general, intractable. 

To help clarify the interaction of the memory term and the orthogonal dynamics, Figure <ref> depicts the memory term in s-t space. In Figure <ref>, the evolution of the solution in time is denoted by the solid blue line at s = 0. To evaluate the memory, solutions to the orthogonal dynamics equation, F(ã(t),s), must be evolved in psuedo-time s using initial conditions that depend on the solution at time t. This is depicted by the dashed red lines in Figure <ref>. This leads to a three-dimensional surface in s-t space, as seen in Figure <ref>. Evaluation of the memory integral then requires a path integration backwards in time along the dashed-lines in Figure <ref>, yielding the shaded yellow region in Figure <ref>.

A quantity that is of particular interest is the memory kernel evaluated at s=0, which is denoted by the solid blue line.  This term, K(ã(t),0), drives the memory term and is typically leveraged to develop closure models <cit.> within the MZ setting.  A clear derivation for the smooth case is presented in Appendix <ref> and one finds the important result,
















       
             [t]0.95

    K(ã(t),0)     =
    [t]
            ∫_Ω∫_Ω (𝐰̃R')(x)  Π^'(x,y) (R(ũ) - f)(y)  dΩ_y dΩ_x,

where R' = ∂R/∂ũ. 

             
          
       


Remarks


  * Equation <ref> shows that the memory is driven by an orthogonal projection of the coarse-scale residual. If this residual is zero, no information is added to the memory. Further, if the coarse-scale residual is fully resolved, no information is added to the memory.

  * The form of the memory at s=0 displays many similarities to Eq. <ref>. This will be further discussed in Section <ref>.








 §.§ Models for the Memory

The construction of an appropriate surrogate to the memory term requires an understanding of the structure of the orthogonal dynamics equation. Due to the challenges associated with the solution of very high-dimensional partial differential equations, to the knowledge of the authors no direct attempt has been undertaken to solve the orthogonal dynamics equation. The most general attempt to extract the memory term and orthogonal dynamics is presented in <cit.>, where Hermite polynomials and Volterra integral equations are used to approximate the memory (and hence orthogonal dynamics). This procedure was shown to provide a reasonably accurate representation of the memory for a low-dimensional dynamical system. The procedure, however, is intractable for high-dimensional problems. This fact is exemplified in the work of Bernstein <cit.>, where  the methodology is applied to the Burgers equation. Another notable attempt to extract the memory is the work of Gouasmi et al. <cit.>, where it is assumed that the semi-group emerging from the orthogonal dynamics equation is a composition operator. This allows the orthogonal dynamics to by solved by virtue of an auxiliary set of ordinary differential equations. This method was shown to be exact for linear systems. It further provided reasonable results for mildly non-linear problems and suggested the presence of finite memory effects. The success of the method, however, is problem-dependent and its accuracy is challenging to assess, from a theoretical standpoint.

Despite the complexity and minimal understanding of the orthogonal dynamics, various surrogate models for the memory exist. These models are typically based on series expansions or geometrical arguments and have been applied to problems of varying complexity. While a comprehensive review of all MZ-based models is outside the scope of this paper, the following subsections will outline a variety of MZ-based models and discuss their mathematical underpinnings.



  §.§.§ The t-model and renormalized variants

The t-model devised by Chorin <cit.> approximates the memory using a zeroth order Taylor expansion of the memory about s=0,

    ∫_0^t K(ã(t-s),s) ds≈ t K(ã(t),0),

where K(ã(t),0) = e^t LPLQLã_0. The t-model may be alternatively derived via a rectangular quadrature rule starting from s=0, or by approximating the orthogonal dynamics evolution operator with e^t QL≈ e^t L. The t-model has traditionally been viewed as a long memory model <cit.>, and has been applied with varying degrees of success to Fourier-Galerkin solutions of the Burgers equation, Euler equations, and incompressible Navier-Stokes equations.

More recently, it has been recognized that the t-model requires an additional scaling to maintain accuracy. A class of renormalized models derived by Stinis <cit.> approximates the memory with a series expansion of the form

    ∫_0^t K(ã(t-s),s) ds ≈∑_j=1^N C_j (-1)^j+1t^j/j!Pe^t L (PL)^jQL𝐚_0k.

 Stinis <cit.> also presents a methodology to determine the coefficients C_j using certain norms of the full solution. 



  §.§.§ The τ-model

The τ-model is motivated by the idea that memory has a finite support in time and makes the approximation,

    ∫_0^t K(ã(t-s),s) ds ≈τ K(ã(t),0),

where τ is the memory length. A dynamic procedure to compute τ was presented by Parish and Duraisamy <cit.>. The methodology leverages the Germano identity and assumes scale similarity to construct an energy transfer constraint between two-levels of coarse-graining.



  §.§.§ Finite Memory Models

A hierarchy of models that again leverage the idea of finite memory was constructed by Stinis <cit.>. The idea is to repeatedly differentiate the memory with respect to time to obtain a set of hierachical equations. This hierarchy is then truncated. The first order model form that emerges from this procedure is

    ∫_0^t K(ã(t-s),s) ds ≈M(t),

where M(t) obeys the equation

    dM/dt = -2 /τM + 2 e^t LPLQLã_0 ,

with τ again being a memory length. 


  §.§.§ Complete Memory Models

A class of models has recently been developed by Price and Stinis <cit.> that assumes an almost commutativity of of PL and QL to obtain an expansion of the memory term that is demonstrated to be accurate for long time intervals. The second order expansion is given by

    ∫_0^t K(ã(t-s),s) ds ≈ t P e^t LPLQLã_0 - t^2/2P e^t L[ PLPLQLã_0 - PLQLQLã_0 ] .



  §.§.§ Faber Approximations

As a final example, Zhu and Venturi <cit.> consider approximations to the memory that are based on a Faber series expansion of the orthogonal dynamics,

    e^t QL =  ∑_j=0^∞a_j(t) F_j(QL),

where F_j is the j_th order Faber polynomial and a_j are the basis coefficients. In addition, <cit.> provides a cogent summary of several additional expansion techniques used to approximate the memory.






 §.§ Discussion

A number of researchers have examined the application of Mori-Zwanzig-based models to the Navier–Stokes equation in a  spectral Fourier-Galerkin setting <cit.>. The discussion to this point presents the general basis for these methods in the context of the Variational Multiscale Method. The demonstration that the memory is driven by an orthogonal projection of the residual is a hitherto unappreciated fact that provides  new insight into the mechanics of MZ-based models.
In addition, a more subtle but equally important result emerged from the previous discussion. 
Since the fine-scale space is infinite dimensional, the projection of the residual onto the fine-scale space may be computed as 

    Π^' (f - R( ũ )) = ( f - R( ũ ) ) - Π̃ (f - R( ũ )).

The computation of the model form only requires that one can represent the coarse-scale residual on the numerical grid.  Fine-scale basis functions do not have to be processed. The fact that this term is a projection is significant. The traditional wisdom in MZ is that one has to define an a priori finite-dimensional set of ODEs for the fine-scale equation. The resulting MZ model is then a function of the (somewhat arbitrary) definition of the finite-size fine-scale equation as well as the fine-scale test/trial functions. Here we have shown that, in the general setting where the fine-scale equation is infinite dimensional, the memory at s=0 can be computed using only coarse-scale information. 

The discussion to this point has been limited by the assumptions of smoothness, which is only appropriate for methods utilizing globally smooth basis functions. The Fourier-Galerkin spectral method is one such example. These methods are limited to canonical domains and are not sufficiently robust for the complex flows and domains commonly encountered in science and engineering. In Section <ref>, we extend the discussion to the case of finite elements.




§ MZ-VMS FORMULATION FOR THE FINITE ELEMENT METHOD

In this section, we develop the formulation of the MZ-VMS framework for finite element methods. The following discussion will encompass both continuous and discontinuous Galerkin methods. As before, denote the test and trial space as V = Ṽ⊕V^', where V≡ L^2(Ω). We have relaxed the H^1(Ω) restriction to allow for discontinuous solutions.
Define T_h to be the decomposition of the domain Ω into a set of non-overlapping elements, T, over Ω_k with boundaries Γ_k. We seek solutions for the coarse-scales in the finite element space

    Ṽ = { w ∈ L^2(Ω) :  v|_T∈ P^k(T),  ∀ T ∈T_h },

where P^k is the space of polynomials up to degree k.
The above definition allows for discontinuities between elements. The continuous Galerkin formulation makes use of a subspace of Eq. <ref> that enforces continuity between elements. 
Some notation is beneficial before proceeding. Define,

    (·,·)_Ω = ∑_k ∈T (·,·)_Ω_k,      (·,·)_Γ = ∑_k ∈T (·,·)_Γ_k.

The FEM weak formulation of Eq. <ref> becomes

    (w, u_t )_Ω + (w, R(u) )_Ω+ ( w , b(u))_Γ  =  (w, f )_Ω    ∀ w ∈V,

where b(u) is a boundary operator that arises via integration-by-parts. In the context of discontinuous Galerkin, one may consider Eq. <ref> to be the “DG strong form". The finite element method proceeds to approximate the solution u in a finite dimensional subspace, Ṽ⊂V. With V^'Ṽ, we obtain the multiscale weak formulation,

    (w̃, ũ_t )_Ω +  (w̃, R(ũ) )_Ω + (w̃, R(u) - R(ũ) )_Ω + ( w̃ , b(ũ))_Γ +  ( w̃ , b(u) - b(ũ))_Γ =  (w̃, f )_Ω    ∀w̃∈Ṽ,


    (w^', ũ_t )_Ω +  (w^', R(ũ) )_Ω + (w^', R(u) - R(ũ) )_Ω + ( w^' , b(ũ))_Γ +  ( w^', b(u) - b(ũ))_Γ =  (w^', f )_Ω    ∀w^'∈V^'.

To illustrate the Mori-Zwanzig procedure, we consider the discrete system resulting from Eqns. <ref> and <ref>. As the basis functions in FEM are typically not orthonormal, we retain the mass matrices in what follows. Equations <ref> and <ref> become,

    d ã/dt =  𝐌̃^-1[  -(𝐰̃, R(ũ) )_Ω - (𝐰̃, R(u) - R(ũ) )_Ω +  (𝐰̃, f )_Ω - ( 𝐰̃ , b(ũ))_Γ -  ( 𝐰̃ , b(u)- b(ũ))_Γ],


    d 𝐚^'/dt =  𝐌^'^-1[  -(𝐰^', R(ũ) )_Ω - (𝐰^', R(u)- R(ũ) )_Ω +  (𝐰^', f )_Ω - ( 𝐰^' , b(ũ))_Γ -  ( 𝐰^' , b(u) - b(ũ))_Γ],

where the mass matrices are

    𝐌̃ = (𝐰̃, 𝐰̃^T)_Ω,     𝐌^' = (𝐰^', 𝐰^'^T)_Ω.

Note that there is no coupling between the mass matrices as a  result of  L^2 orthogonality of the coarse and fine-scales. Compactly, we can express the entire system as

    d 𝐚/dt =  𝐌^-1[  -(𝐰, R(u) )_Ω - ( 𝐰 , b(u))_Γ +  (𝐰, f )_Ω],

with 𝐌 = (𝐰, 𝐰^T)_Ω.
Through the Mori-Zwanzig procedure, we can integrate out the fine-scale variable and express Eq. <ref> as

    (𝐰̃, ũ_t )_Ω +  (𝐰̃, R(ũ) )_Ω + ( 𝐰̃ , b(ũ))_Γ =  (𝐰̃, f )_Ω  +  𝐌̃∫_0^t K(𝐚̃(t-s),s) ds.

The memory kernel is given by

    K_j(𝐚_0,t) = PLe^t QLQLa_0j,

where 

    L≡∑M_ij^-1[  -(w_j, R(u_0) )_Ω- ( w_j , b(u_0))_Γ  +  (w_j, f )_Ω] ∂/∂ a_0i .

The summation is over all global unknowns. The memory kernel at s=0 is again of interest. The derivation closely follows that of the smooth case and one finds,


       
             [t]0.95

    𝐌̃K(ã(t),0)     =
    [t]
            ∫_Ω∫_Ω (𝐰̃R')(x)  Π^'(x,y) (R(ũ) - f)(y)  dΩ_y dΩ_x  
       + ∫_Ω∫_Γ (𝐰̃R')(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΩ_x 
       + ∫_Γ∫_Ω (𝐰̃ b')(x)   Π^'(x,y) (R(ũ) - f)(y) dΓ_y dΩ_x 
       + ∫_Γ∫_Γ (𝐰̃  b')(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΓ_x,

where again

    R' = ∂R/∂ũ,      b' = ∂ b/∂ũ.

             
          
       



Remarks


  * Compared to the smooth case, the finite element method gives rise to additional interactions between the coarse and fine-scales. 

  * The coarse-scale equation, Eq. <ref>, is again non-local in time.

  * The memory kernel at s=0 contains both volumetric and surface integrals. The form of this term is again similar to that  obtained by Hughes <cit.> for the “rough" FEM case. This will again be discussed in Section <ref>.






§ CONNECTIONS OF MZ-VMS WITH EXISTING CONCEPTS

In Section <ref>, it was seen that all models utilize the term  K(ã(t),0), which is written equivalently as e^t LPLQLã_0. This value drives the memory term and has been discussed throughout the previous sections. In this section we investigate the mechanics of this term in detail and draw parallels between memory effects, artificial viscosity, and upwinding.

We consider the memory kernel at s=0 for the FEM case. Recall that this term appears as

    𝐌̃K(ã(t),0)     = 
    [t]
            ∫_Ω∫_Ω (𝐰̃R')(x)  Π^'(x,y) (R(ũ) - f)(y)  dΩ_y dΩ_x  
       + ∫_Ω∫_Γ (𝐰̃R')(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΩ_x 
       + ∫_Γ∫_Ω (𝐰̃ b')(x)   Π^'(x,y) (R(ũ) - f)(y) dΓ_y dΩ_x 
       + ∫_Γ∫_Γ (𝐰̃  b')(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΓ_x,

where R' = ∂R/∂ũ and b' = ∂  b /∂ũ. The term Π^' is the L^2 projection onto the fine-scales,

    Π^'(x,y) =  𝐰^'^T(x)𝐌^'^-1𝐰^'(y).

In the globally smooth case, the terms involving the boundary operators drop and only the first line of Eq. <ref> is retained. 

Remarks


  * The memory is driven by the residual of the coarse-scales projected onto V^'. When the residual of the coarse-scales is negligible, no additional information is added to the memory. Models such as the t and τ-model are inactive. Further, if the coarse-scale residual is non-zero but exists only in Ṽ, models such as the t and τ-model are again inactive and no information is added to the memory. 

  * The orthogonal projection,  Π^', can be conceptualized as a mechanism to  enforce the approximation to constrain the fine-scale state to  the correct trial space, V^'. A significant body of work on orthogonal subgrid-scale models in the context of the Variational Multiscale Method has been undertaken by Codina <cit.>.

  * The boundary terms in the FEM formulation give rise to surface integrals. As will be shown later, these surface integrals can, in turn, give rise to jump operators between elements and can add artificial diffusion to the system.

  * It is seen that the orthogonal projector, Π^'(x,y), can be viewed as an approximation to the fine-scale Green's function. This will be discussed in the next section.




 §.§ The τ-model and an orthogonal approximation to the fine-scale Green's function

Approximating the memory with the τ-model gives rise to the following closed equations for the coarse-scales,

    (𝐰̃, ũ_t) + (𝐰̃, R(ũ) )  -   τ𝐌̃ K(ã(t),0) = (𝐰̃, f ).

To provide insight into this formulation, let R be a linear operator, R = L. The steady-state solution achieved by the τ-model is

    (𝐰̃, Lũ )  -   τ𝐌̃ K(ã(t),0) = (𝐰̃, f ).

Integrating by parts, the memory at s=0 may be expressed as

    𝐌̃K(ã(t),0)     = 
    [t]
            ∫_Ω∫_Ω (L^*𝐰̃)(x)  Π^'(x,y) (L ũ - f)(y)  dΩ_y dΩ_x  
       + ∫_Ω∫_Γ (L^*𝐰̃)(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΩ_x 
       + ∫_Γ∫_Ω (b^*𝐰̃ )(x)   Π^'(x,y) (Lũ - f)(y) dΓ_y dΩ_x 
       + ∫_Γ∫_Γ (b^* 𝐰̃)(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΓ_x,

where b^* is the corresponding boundary operator produced through the integration by parts.
The memory term has the same form as that obtained by Hughes <cit.>. In the steady case, the orthogonal projection can be viewed as an approximation to the fine-scale Green's function,

    g^'(x,y) ≈τ(x,y) Π^'(x,y).


The τ-model can be  further viewed as a variant of the adjoint stabilization method, where the only difference is the inclusion of the orthogonal projector. To provide insight into this approximation, we consider the fine-scale Green's functions for the advection diffusion operator in one dimension,

    L = c ∂/∂ x - ν∂^2/∂ x^2.

The fine-scale Green's function is examined for the continuous Galerkin method with first order finite elements using the approach of Hughes and Sangalli <cit.>. We consider the case where x ∈ (0,1) with 16 elements and physical parameters ν = 0.001, c = 1.
Figure <ref> compares the orthogonal approximation to the Green's function (with τ=1) to the exact fine-scale Green's function. It  is seen that the orthogonal Green's function leads to an approximation that modulates the local residual. In contrast to the true Green's function, it is seen that the approximate orthogonal Green's function is highly localized. Further, it should be noted that the approximate Green's function does not include any physics related to the associated governing equations.









 §.§ Residual-Based Artificial Viscosity

To further clarify the role of the e^t LPLQLã_0 term, we consider the hyperbolic conservation law,

    ∂𝐮/∂ t +  ∇·𝐅(𝐮) = 0     in    Ω.

The semi-discrete system obtained through the FEM discretization is,

    ∫_Ω𝐰𝐮_t dΩ + ∫_Ω𝐰∇·𝐅(𝐮) dΩ +  ∫_Γ𝐰𝐛(𝐮,𝐧) dΓ = 0,

where again 𝐛 is a boundary operator, and 𝐧 is the normal vector at element interfaces. Application of the MZ-VMS procedure leads to

    ∫_Ω𝐰̃𝐮_t dΩ + ∫_Ω𝐰̃∇·𝐅(ũ) dΩ + ∫_Γ𝐰̃𝐛(ũ,𝐧) dΓ  =  𝐌̃∫_0^t K(ã(t-s),s)ds.



       
             [t]0.95
The value of the memory at s=0 can be expressed as,
 
    𝐌̃K(ã(t),0) =  ∫_Ω𝐰̃∇·𝐅'(𝐪) dΩ +  ∫_Γ𝐰̃𝐛'(𝐪, 𝐧) dΓ,

where 𝐪 is given by,

    ∫_Ω𝐰^'𝐪 dΩ =  ∫_Ω𝐰^'∇·𝐅(ũ) dΩ +  ∫_Γ𝐰^'𝐛(ũ,𝐧)  dΓ.

The term 𝐅' = ∂𝐅/∂ũ is the flux Jacobian and 𝐛' is the numerical flux function linearized about 𝐮̃. The resulting coarse-scale equation for the τ-model is

    ∫_Ω𝐰̃𝐮_t dΩ + ∫_Ω𝐰̃∇·( 𝐅(ũ) - τ𝐅'(𝐪))dΩ + ∫_Γ𝐰̃( 𝐛(ũ,𝐧) - τ𝐛'(𝐪,𝐧) ) dΓ = 0.

             
          
       



Consider now Eq. <ref> augmented with an artificial viscosity term that is proportional to the orthogonal projection of the divergence of the flux,


    ∂𝐮/∂ t + ∇·𝐅 = τ∇·𝐅'( Π^'∇·𝐅).

A standard discretization technique for this second order equation is to split it into two first order equations <cit.>,

    ∂𝐮/∂ t + ∇·(𝐅(𝐮) - 𝐅'(𝐪)) = 0,

with

    𝐪 =Π^'∇·𝐅(𝐮).

Assuming that the boundary operators are handled analogously, the discretization of Eq. <ref> and Eq. <ref> through finite elements leads to precisely Eqns. <ref> and <ref>.

Remarks


  * For a hyperbolic conservation law, the memory is driven by a non-linear term that acts as a type of non-linear artificial viscosity. 

  * The magnitude of the artificial dissipation is proportional to the projection of the flux onto the fine-scales. If the flux term is fully resolved, no information is added to the memory.

  * Due to the appearance of the orthogonal projector, it is difficult to comment on the sign of the artificial viscosity. While proofs exist showing that the term e^t LPLQL𝐚̃_0 is globally dissipative in certain settings <cit.>, no such result is readily apparent in the general case.






 §.§ Global Dissipation of the Memory Kernel

An elegant proof regarding the stability of the t-model was delivered by Hald and Stinis in <cit.>. This proof is presented here with a slight extension.. We examine the impact of the memory kernel at s=0 (which is closely related to the τ-model) on the global energy, J. Define J as,

    J = 1/2∫_Ω u^2 dΩ = 1/2∫_Ω( ∑_j a_j w_j(𝐱)  )^2 dΩ.

To derive an equation for the evolution of the energy, define the inner product,

    ( 𝐚,𝐛)_J =   ∫_Ω (𝐰^T 𝐚) (𝐰^T 𝐛) d Ω.

 The inner product in Eq. <ref> produces the volume integral of the physical space functions that 𝐚 and 𝐛 are representing. Note that for an orthonormal basis,
  
    ∫ w_j (𝐱) w_i(𝐱) d𝐱 = c_ijδ_ij.

 With w_k being orthogonal, the inner product in Eq. <ref> can be written as
 
    ( 𝐚,𝐛)_J =   ∫_Ω (𝐰^T 𝐚) (𝐰^T 𝐛) d Ω =  ∑_j c_j a_j b_j.

Since the inner product can be written as a simple sum, we have the following important property

    (𝐚 + 𝐛, 𝐜 + 𝐝)_J = (𝐚, 𝐜)_J  + (𝐛, 𝐝)_J.

We now proceed to demonstrate the energy transfer properties of the τ-model. 
We can write an equation for the rate of change of the energy by taking the inner product with 𝐚,

    ( 𝐚 , d 𝐚/dt)_J = ( 𝐚̃ , 𝐟( 𝐚̃, 𝐚^' ) )_J +  ( 𝐚̂ , 𝐠( 𝐚̃, 𝐚^' ) )_J,

or

    d J/dt  = ( 𝐚̃ , 𝐟( 𝐚̃, 𝐚^' ) )_J +  ( 𝐚̂ , 𝐠( 𝐚̃, 𝐚^' ) )_J,

In writing Eq. <ref>, we have made use of the property in Eq. <ref>. Consider the case where the original numerical method conserves the quantity J, 

    d J/dt = (𝐚̃ ^i , 𝐟^i( 𝐚̃, 𝐚^' ) ) +  ( 𝐚̂ ^i , 𝐠^i( 𝐚̃, 𝐚^' ) )
    = 0.

Setting the unresolved variables to zero shows

    ( 𝐚̃^i , 𝐟^i( 𝐚̃, 𝐚^' ) ) = 0.

Similarly, let ã = ϵ𝐡, 

    ( 𝐚̃^i , 𝐟^i( 𝐚̃, ϵ𝐡) ) +  ( ϵ𝐡^i, 𝐠^i( 𝐚̃,  ϵ𝐡) )
    = 0.

Expanding in a Taylor series gives

    ( 𝐚̃^i , 𝐟^i( 𝐚̃, 0 ) )  +  ϵ( 𝐚̃^i , ∂𝐟^i/∂𝐚^'𝐡) + ϵ( 𝐡^i, 𝐠^i( 𝐚̃,  0) )
    + O(ϵ^2) = 0.

Dividing by ϵ, letting ϵ→ 0, and setting 𝐡 = 𝐠( 𝐚̃,  0) gives

    ( 𝐚̃^i , ∂𝐟^i/∂𝐚^'𝐠 ( 𝐚̃,  0)  ) = - ( 𝐠^i( 𝐚̃,  0) , 𝐠^i( 𝐚̃,  0) )

Noting that the memory term at s=0 can be written for the coarse-scale equation as,

    K(ã(t),0) = ∂𝐟^i/∂𝐚^'𝐠 ( 𝐚̃,  0),

One obtains

    ( 𝐚̃^i ,K(ã(t),0)  ) = - ( 𝐠^i( 𝐚̃,  0) , 𝐠^i( 𝐚̃,  0) )  ≤ 0 .

It is seen that, globally, the memory term at s=0 dissipates energy. This dissipation rate is proportional to the norm of the right hand side of the fine-scale equation evaluated at u^'=0. This is best thought of as a Leonard stress type term.

    d 𝐚̃^i/∂ t = 𝐟^i( 𝐚̃, 0 )  + τ∂𝐟^i/∂𝐚̂ g(ã,0) .

Taking the inner product of Eq. <ref> with ã^i and inserting Eq. <ref>, we have for the t-model and τ-model,

    d J̃/∂ t =  - τ( 𝐠^i( 𝐚̃,  0) , 𝐠^i( 𝐚̃,  0) ) ≤ 0.

Equation <ref> demonstrates that, for a numerical method conserving energy for the i_th variable, the MZ model will dissipate energy for that variable. 




 §.§ Relationship to the Upwind Flux

The relationship between artificial viscosity and upwinding is a well-recognized aspect in numerical methods. In the previous section, it was demonstrated that the memory is driven by a term that resembles artificial viscosity. This term can be explicitly linked to upwinding. To demonstrate this, we consider the τ-model applied to the discontinuous Galerkin discretization of the linear advection equation. The complete derivation of what follows is provided in Appendix <ref>.

The linear advection equation is given by,

    ∂ u/∂ t +  ∂ f/∂ x = 0,

with f = cu.
The discontinuous Galerkin discretization of Eq. <ref> leads to the following weak formulation on the k_th element,

    ∫_Ω_k𝐰_k ∂ u/∂ t d Ω - ∫_Ω_k∂𝐰_k/∂ x f dΩ = -∫_Γ_k𝐰_k f^*(u,n) dΓ,

where f^* is a numerical flux function that provides coupling between the elements. Through the MZ-VMS framework, the τ-model leads to the following coarse-scale equation, 

    ∫_Ω_k𝐰̃_k ∂ũ/∂ t d Ω - ∫_Ω_k∂𝐰̃_k/∂ x f(ũ )dΩ   = -∫_Ω_k𝐰̃_k f^*(ũ,𝐧) dΓ  + τ𝐌̃K_k(ã(t),0).

The memory term depends on the numerical flux function. In the case where the numerical flux function is a central flux, the memory at s=0 (Eq. <ref>) is found to be

    𝐌̃K_k(ã(t),0)  = 
    S_1 𝐰̃_k^R  c^2/2(   u_k+1^L - u_k^R )
    - S_1 𝐰̃_k^L  c^2/2( u_k^L - u_k-1^R  ),

where u_k^L and u_k^R  are the values of u on the left right edges of the k_th element, respectively (see Figure <ref>).

The constant S_1 is a scalar whose value is given in Appendix <ref>. Noting that the natural time scale for linear advection is |c|^-1, a possible selection for the memory length in the τ model is

    τ = 1/|c| S_1,

which leads to


       
             [t]0.95

    τ𝐌̃K_k(ã(t),0)  = 
    𝐰̃_k^R  |c|/2( u_k+1^L - u_k^R  )
    -𝐰̃_k^L  |c|/2(u_k^L  - u_k-1^R  ).

It is seen that in this setting, the τ-model provides precisely the upwind correction for linear advection. 

             
          
       























Remarks


  * The above exercise demonstrates the role of upwind fluxes as implicit subgrid-scale models. Notably, it was formally proved that upwind fluxes can be derived through the approximation of memory effects.

  * In the non-linear case, the first three terms in Eq. <ref> will not drop. Instead, the memory term is driven by both volumetric and surface integrals. 

  * The last line of Eq. <ref> will manifest itself in the form of a flux function in the non-linear case. The linearized operator b^* will contain flux Jacobians, and similarities may be drawn between this portion of the memory term and approximate Riemann solvers.






 §.§ Discussion

The previous section outlined the development of the MZ-VMS framework to the discontinuous Galerkin method. The addition of flux terms in the DG approach led to additional interactions between the coarse and fine-scales. The resulting memory term was seen to be driven by a non-linear artificial viscosity term. It is important to note that the treatment of the viscous and inviscid fluxes in this setting are identical. It was further shown that, in the case of a central flux in one-dimension, a portion of this term resembles a numerical flux function. This gives rise to the following important point: an MZ-VMS model should be constructed off of a non-dissipative flux function. Indeed, it has been recognized that flux schemes such as the Roe flux are, themselves, implicit subgrid-scale models <cit.>. Constructing an MZ-VMS subgrid-scale model from a scheme that already contains an implicit subgrid-scale model appears to be unwise.

Many physical conservation laws are known to obey the conservation of an auxiliary quantity known as entropy. We believe that a numerical method that obeys this conservation law is the most robust starting point for the development of an MZ-VMS model.  Fourier-Galerkin methods for the incompressible Navier-Stokes equations, where most MZ research has been performed, can shown to be entropy (which is also the kinetic energy) preserving. For discontinuous Galerkin methods, the development of entropy-conservative schemes for square entropy functions is straightforward. The τ-model constructed on such an entropy conserving scheme can be shown to be globally dissipative via a straightforward extension of the proof provided in <cit.>. The development of entropy-conservative DG methods for non-square methods is less straightforward, as is the development of the MZ-VMS idea.

 


§ NUMERICAL EXAMPLE

We consider a simple numerical example to demonstrate the potential of a VMS-MZ model. The example considered here is the viscous Burgers equation in one-dimension,

    ∂ u/∂ t + 1/2∂ u^2/∂ x = ν∂^2 u/∂ x^2.

The viscous Burgers equation is a one-dimensional model for the Navier-Stokes equations and has been a testing ground for many subgrid models. We consider high-order discontinuous Galerkin simulations of Eq. <ref> on the domain x ∈ [0,2π). We consider the initial condition

    u(x) = U_0^* ∑_i=1^k_c√(2 E(k_i))sin(k_i x + β_i),

where E(k) = 5^-5/3 if 1≤ k ≤5 and E(k) = k^-5/3 for k > 5. Eq. <ref> initializes an energy spectrum with a -5/3 slope for k > 5. The phase angle β is a random number in the domain [-π,π]. A constant seed value is used in all of the simulations. The simulations considered take k_c = 32 and are performed with a viscosity of ν = 0.005.



 §.§ Construction of the τ-model

The discontinuous Galerkin weak form for Burgers equation with the τ-model is

    (𝐰̃_k, ũ_t )_Ω_k -  ( ∂/∂ x𝐰̃_k, 1/2ũ^2 - τũΠ^' q )_Ω_k=  -( 𝐰̃_k ,  F̃^*(ũ) - τA^*( Π^' q))_Γ_k

with

    (𝐰,q)_Ω_k =   ( ∂/∂ x𝐰_k , 1/2ũ^2)_Ω_k  -  ( 𝐰_k ,    F^*(ũ))_Γ_k.

The numerical flux function is taken to be an entropy-conserving flux <cit.>,

    F^*(u^+,u^-) = 1/6(u^+^2 + u^+ u^- + u^-^2).

The linearized flux function is

    A^*(u'^+,u'^-) = 1/6(2u^+ u'^+  +  u^- u'^+ + u^+ u'^- + 2 u^-u'^-).

Note that the auxiliary variable q only needs to be solve for q ∈V^'.


 §.§ Results

Simulations are performed using an eighth-order discontinuous Galerkin method with 8 elements. We consider simulations that use an entropy conserving flux (which we denote as no model), with the entropy conserving flux plus the τ-model, and with the Roe flux. Figure <ref> shows the resulting profiles for the kinetic energy, dissipation rate, physical space profiles, and energy spectra. Simulations performed with no implicit subgrid-scale model (i.e. with the entropy conserving flux) are seen to not dissipate enough kinetic energy. This leads to the formation of unphysical oscillations and a build up of energy at high frequencies. The Roe flux is capable of acting as a subgrid-scale model, but is also seen to be slightly under-dissipative. This is evident in the physical space velocity profiles in several locations, with the most notable being around x=1.5. The τ-model provides further improved predictions. The evolution of kinetic energy is reasonably captured and the physical space velocity profiles are improved upon from the Roe flux.










§ NUMERICAL EXAMPLE

Discontinuous Galerkin simulations of the one-dimensional Euler equations are considered to demonstrate the optimal prediction-based models. The Euler equations are given by

    ∂ρ/∂ t + ∂/∂ x ( ρ u) = 0,


    ∂ρ u/∂ t + ∂/∂ x (ρ u^2 - p) = 0,


    ∂ρ E/∂ t + ∂/∂ x (u(E + p)) = 0.

The simulation of the Shu-Osher problem is considered. The Shu-Osher problem consists of a Mach 3 shock interacting with a sinusoidal density wave. The accurate simulation of the flow requires methods that can handle discontinuities without dissipating the small-scale density fluctuations. The problem setup is given by,

    ρ = {[        3.857143          x < -4; 1 + 0.2 sin(5x)          x ≥ -4;                 ]     
    u = {[ 2.629369   x < -4;        0   x ≥ -4;          ]    
    p= {[ 10.333333    x < -4;         1    x ≥ -4;           ].

High-order methods are known to perform poorly for shocks. Artificial "ringing" across a shock front with high order polynomials is a significant issue and often leads to failure of a numerical method. This has motivated a significant body of work for shock capturing methods for discontinuous Galerkin <cit.>. Stability issues are also known to occur for smooth flows when the simulation is under-resolved. In this case smooth features become aliased and spurious oscillations can occur. This often leads to the failure of the numerical method as well. Here we investigate the capabilities of the τ-model as a general method to handle both shocks and otherwise smooth subgrid phenomena. The Roe flux is used to model the inter-element flux function arising from the MZ procedure.

Simulations using a total of 400 degrees of freedom at polynomial orders p=0,1, and 3 are considered. At this resolution, baseline discontinuous Galerkin simulations fail for the p=1 and p=3 cases. The inclusion of the MZ-τ-model stabilizes the simulation. The resulting solution predictions at t=2.0 are shown in Figure <ref>. The p=1 and p=3 simulations contain spurious oscillations around the shock fronts, but accurately capture the location of the shock-fronts as well as the sinusoidal post-shock density waves. 




§ CONCLUSIONS



This work outlined the development of a multiscale modeling framework that combines the Mori-Zwanzig formalism of statistical mechanics with the Variational Multiscale method. We refer to this framework as MZ-VMS. Although the MZ formalism has previously been recognized as a powerful tool for model reduction, the formulation of MZ-based approaches for the numerical simulation of partial differential equations using practical discretization techniques such as finite elements has remained unclear. The use of MZ-VMS is established as a practical tool for model development within the spectral and finite element setting and elegant connections were developed with existing modeling strategies. 

The MZ-VMS approach leverages an exact transformation of the discretized equations into phase-space, allowing for a systematic modeling approach that is applicable to both unsteady linear and non-linear systems. The method results in an exact coarse-scale equation where the effect of the unresolved scales on the resolved scales are represented in the form of a memory integral. The exact evaluation of the memory term is, however, not tractable as it depends on the orthogonal dynamics, which can only be computed by solving a very high dimensional PDE in phase space. Various strategies for modeling the memory were discussed. 

We demonstrated the following key results:





    
  * Similar to the VMS Green's function, the memory term is driven by both an orthogonal projection of the coarse-scale residual and jumps at element interfaces.  This discovery lends  insight into MZ-based models and provides the first links between MZ-based methods and existing stabilization techniques.

    
  * In the general setting where the fine-scale equation is infinite dimensional, the memory at s=0 can be computed using only coarse-scale information. This is a critical observation for model development.

    
  * In the steady linear case, the MZ-τ-model leads to a residual-based closure that may be viewed as a variant of the adjoint stabilization method. 

    
  * In the case of hyperbolic conservation laws, there are inherent links between the memory term, artificial viscosity, and upwinding. 

    
  * In particular, we showed that the MZ-τ-model can recover an upwind flux for discontinuous Galerkin discretizations of linear advection. 

    
  * In the non-linear case, the MZ-τ-model will lead to flux functions that contain flux Jacobians, and similarities may be drawn between this portion of the memory term and approximate Riemann solvers. 



The MZ-VMS framework can serve as a powerful tool for multiscale modeling. Models derived from the MZ-VMS procedure are typically obtained through mathematically consistent arguments (rather than phenomenological arguments) and can readily be applied to complex multiphysics systems. As an example, we have obtained promising results in the context of Large Eddy Simulations of turbulent flows <cit.>. The framework can also be naturally extended to applications in projection-based reduced order modeling and error estimation. 






§ ACKNOWLEDGMENTS
 


This research was funded by the AFOSR under the project LES Modeling of Non-local effects using Statistical Coarse-graining (Tech.  Monitor:  Jean-Luc Cambier). 





§ DERIVATION OF THE MEMORY AT S=0 (SPECTRAL METHOD CASE)

A variety of Mori-Zwanzig-based models require evaluating the memory term at s=0, 

    K(ã(t),0) = e^t LPLQLã_0.

Evaluation of PLQLã_0 is an exercise in algebra, and we provide the step-by-step process for the smooth case where jump terms between the elements are not considered. We consider the non-linear equation,

    ∂ u/∂ t + R(u) = f.

The evaluation of PLQLã_0 is as follows:


  * First compute e^t LL𝐚̃_0. This is simply the right hand side of the coarse-scale equation,

    e^t LL𝐚̃_0 =  (𝐰̃,  f - R( ũ )  ) -(𝐰̃,  R(u) - R( ũ )  ).


  * Next compute e^t LQL𝐚̃_0. By definition of how we've written the right hand side,

    e^t LQL𝐚̃_0 =  -(𝐰̃,  R(u) - R( ũ )  ),

or in terms of the basis coefficients

    e^t LQL𝐚̃_0 =  -(𝐰̃,  R(𝐰^T 𝐚) - R( 𝐰̃^T 𝐚̃ )  ).


  * Compute e^t LLQL𝐚̃_0. This is simplified by recognizing that the Liouville operator is the Frechet derivative in the direction of the right hand side. For some function scalar function g,

    L g(u(𝐚_0)) = ∂ g/∂𝐚_0L𝐚̃_0.

The linearization in the Frechet derivative is with respect to the modal variables, 𝐚_0. For the mapping from u = 𝐰^T 𝐚, the chain rule gives

    ∂ g/∂𝐚_0L𝐚̃_0 = ∂ g/∂ u_0𝐰^T L𝐚̃_0.

The Liouville operator can be applied to QLã_0 by linearizing R with respect to u, and evaluating the resulting model at R(u),

    e^t LLQL𝐚̃_0 =  -(𝐰̃,  R'(𝐰^T  [  (𝐰̃,  f - R( ũ )  ) -(𝐰̃,  R(u) - R( ũ )  )  ]  )  
    
     - R'(𝐰̃^T  [    (𝐰̃,  f - R( ũ )  ) -(𝐰̃,  R(u) - R( ũ )  ) ]    ).

 
  * Compute e^t LPLQL𝐚̃_0. This sets all terms involving 𝐚^' to zero and we obtain

    e^t LPLQL𝐚̃_0 =  -(𝐰̃,  R'(𝐰^T  (𝐰,f - R'( ũ ) )  )  - R'(𝐰̃^T   (𝐰̃,f - R( ũ ) )  )     ).


  * Equation <ref> may be simplified by recognizing that terms appear as projections,

    e^t LPLQL𝐚̃_0 =  -(𝐰̃,  R'( Π (f - R( ũ ) )   )  - R'( Π̃ (f - R( ũ ) )   )  )     ).


  * With Π^' = Π - Π̃ and noting that R' is linear, 

    e^t LPLQL𝐚̃_0 =  -(𝐰̃,  R'( Π^' (f - R( ũ ) )   )     ).






§ DERIVATION OF THE Τ-MODEL FOR LINEAR ADVECTION

This Appendix details the derivation for the τ-model for the discontinuous Galerkin discretization of the linear advection equation in one-dimension using Legendre polynomials. Note that the Legendre polynomials form an orthogonal family of polynomials. The linear advection equation is given by

    ∂ u/∂ t +  ∂ f/∂ x = 0,

with f = cu.
The discontinuous Galerkin discretization of Eq. <ref> leads to the following weak formulation,

    ∫_Ω_k w ∂ u/∂ t d Ω - ∫_Ω_k∂ w/∂ x f(u) dΩ = -∫_Γ_k w f^*(u,𝐧) dΓ,

where f^* is a numerical flux function that provides coupling between the elements.
Separating scales, we can express Eq. <ref> by

    ∫_Ω_kw̃∂ũ/∂ t d Ω - ∫_Ω_k∂w̃/∂ x f(ũ )dΩ - ∫_Ω_k∂w̃/∂ x f( u^') dΩ  = -∫_Γ_kw̃ f^*(ũ,𝐧) dΓ  -∫_Γ_kw̃ f^*(u^',𝐧) dΓ.


    ∫_Ω_kw^'∂u^'/∂ t d Ω + ∫_Ω_kw^'∂/∂ xf( ũ) dΩ + ∫_Ω_kw^'∂/∂ xf( u^') dΩ  = ∫_Γ_kw^' (f - f^*)(ũ,𝐧) dΓ∫_Γ_kw^'(f -  f^*)(u^',𝐧) dΓ.

Note that we have expressed the fine-scale equation, Eq. <ref>, in the DG strong form. Through the MZ-VMS framework, we can express the coarse-scale equation as

    ∫_Ω_kw̃∂ũ/∂ t d Ω - ∫_Ω_k∂w̃/∂ x f(ũ )dΩ   = -∫_Γ_kw̃ f^*(ũ,n) dΓ  + 𝐌̃K(ã(t),0) ,

where M̃ is the coarse-scale mass matrix.
From Eq. <ref> we recall that the memory term at s=0 may be expressed as

    𝐌̃K(ã(t),0)     = 
    [t]
            ∫_Ω∫_Ω (L^*𝐰̃)(x)  Π^'(x,y) (L ũ )(y)  dΩ_y dΩ_x  
       + ∫_Ω∫_Γ (L^*𝐰̃)(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΩ_x 
       + ∫_Γ∫_Ω (b^*𝐰̃ )(x)   Π^'(x,y) (Lũ )(y) dΓ_y dΩ_x 
       + ∫_Γ∫_Γ (b^* 𝐰̃)(x)  Π^'(x,y) b(ũ(y))  dΓ_y dΓ_x,

where L = c ∂/∂ x and the adjoint operator is L^* = -  c ∂/∂ x. The boundary operators on element k correspond to the flux functions,

    ∫_Γ_k b_k(u) dΓ  = -[ f(u_k^R) - f^*(u_k^R,u_k+1^L ) ] + [ f(u_k^L) - f^*(u_k-1^R,u_k^L )]

where f^* is a numerical flux function and u_k^L and u_k^R correspond to the value of u on the left and right boundary of the k_th element (see Figure <ref>). Similarly we have

    ∫_Γ_k b_k^*(u) dΓ  = f^*(u_k^R,u_k+1^L )  -  f^*(u_k-1^R,u_k^L ).

For notational purposes we denote

    Δ f^R_k = f(u^R_k)  - f^*(u_k^R,u_k+1^L ) ,     Δ f^L_k = f(u^L_k)  -  f^*(u_k-1^R,u_k^L )


Due to orthogonality, one can show that all terms in Eq. <ref> involving volume integrals drop and the memory term on the k_th element simplifies to,

    𝐌̃K_k(ã(t),0)  = ∫_Γ_k∫_Γ (b_k^* 𝐰̃_𝐤)(x)  𝐰^'^T(x) 𝐌^'^-1𝐰^'(y) b_k(ũ(y))  dΓ_y dΓ_x.

Evaluating the integral over Γ_y,

    𝐌̃K_k(ã(t),0)     = 
    [t]
        -∫_Γ_k (b_k^* 𝐰̃_k)(x)  𝐰^'_k(x)^T 𝐌^'^-1[ 𝐰^'_k^R Δ f_k^R  - 𝐰^'_k^L Δ f_k^L ]  d Γ_x.

Evaluation of the integral of Γ_x then yields the expression

    𝐌̃K_k(ã(t),0)     = 
    [t]
    -𝐰̃_k^R  f^* (      𝐰^'_k^R^T𝐌^'^-1[ 𝐰^'_k^R Δ f_k^R  - 𝐰^'_k^L Δ f_k^L ]  , 
       𝐰^'_k+1^L^T𝐌^'^-1[ 𝐰^'_k+1^R Δ f_k+1^R  - 𝐰^'_k+1^L Δ f_k+1^L   ) 
    
      + 𝐰̃_k^L  f^* (      𝐰^'_k-1^R^T𝐌^'^-1[ 𝐰^'_k-1^R Δ f_k-1^R  - 𝐰^'_k-1^L Δ f_k-1^L ]  , 
       𝐰^'_k^L^T𝐌^'^-1[ 𝐰^'_k^R Δ f_k^R - 𝐰^'_k^L Δ f_k^L ]   ) 
      .

This can be expressed compactly as

    𝐌̃K_k(ã(t),0)     = 
    [t]
    -𝐰̃_k^R  f^* (      S_1 Δ f_k^R - S_2 Δ f_k^L  , S_2 Δ f_k+1^R -S_1 Δ f_k+1^L  ) 
    
     + 𝐰̃_k^L  f^* (      S_1 Δ f_k-1^R -S_2 Δ f_k-1^L  , S_2 Δ f_k^R  - S_1 Δ f_k^L ) 
      .

where S_1 and S_2 are scalars given by,

    S_1 = 𝐰^'^R^T𝐌^'^-1𝐰^'^R ,      S_2 = 𝐰^'^R^T𝐌^'^-1𝐰^'^L .

Note that for the Legendre polynomials we have

    𝐰^'^L^T𝐌^'^-1𝐰^'^L  = 𝐰^'^R^T𝐌^'^-1𝐰^'^R = S_1.

Further, for Legendre polynomials one will note that the ratio S_2/S_1≈1/N, where N is the number of basis functions. In the limit that the fine-scales have infinite support (N →∞), we may drop the terms involving S_2 to have

    𝐌̃K_k(ã(t),0)  = 
    -𝐰̃_k^R  f^* (   S_1 Δ f_k^R , -S_1 Δ f_k+1^L ) 
      + 𝐰̃_k^L  f^* (   S_1 Δ f_k-1^R  , 
       - S_1 Δ f_k^L   )

It is seen that the memory is dependent on the numerical flux function. Here, we examine the case where the numerical flux is a central flux,

    f^*(u_k^R,u_k+1^L) = c/2( u_k^R + u_k+1^L).

Note that for the central flux we have

    Δ f_k^R = c/2(u_k^R  -  u_k+1^L  ),     Δ f_k^L =  c/2( u_k^L - u_k-1^R  ) .

Further,

    f^*(Δ f_k^R, - Δ f_k+1^L ) = c^2/2(u_k^R -  u_k+1^L ),      f^*(Δ f_k-1^R, - Δ f_k^L ) = c^2/2(u_k-1^R -  u_k^L ).

We see that Eq. <ref> reduces to 

    𝐌̃K_k(ã(t),0)  = 
    -S_1 𝐰̃_k^R  c^2/2(  u_k^R - u_k+1^L )
    + S_1 𝐰̃_k^L  c^2/2(u_k-1^R -  u_k^L  ).

Noting that the natural time scale for linear advection is |c|^-1, a possible selection for the memory length in the MZ-τ-model is

    τ = 1/|c| S_1,

which leads to

    τ𝐌̃K_k(ã(t),0)  = 
    𝐰̃_k^R  |c|/2( u_k+1^L - u_k^R  )
    -𝐰̃_k^L  |c|/2(u_k^L  - u_k-1^R  ).

Equation <ref> is precisely the upwind correction for linear advection. 






10


hughes0
Hughes, T. J., Multiscale phonomena: Green's functions, the
  Dirichlet-to-Neumann formulation, subgridscale models, bubbles and the
  origins of stabilized methods, Computer methods in applied mechanics
  and engineering, Vol. 127, 1995, pp. 387–401.

hughes1
Hughes, T. J., Feijoo, G., Mazzei, L., and Qunicy, J., The variational
  multiscale method - a paradigm for computational mechanics, Computer
  methods in applied mechanics and engineering, Vol. 166, 1998,
  pp. 173–189.

hughes2004multiscale
Hughes, T. J., Scovazzi, G., and Franca, L. P., Multiscale and
  stabilized methods, Encyclopedia of computational mechanics, 2004.

codina2017variational
Codina, R., Badia, S., Baiges, J., and Principe, J., Variational
  multiscale methods in computational fluid dynamics, Encyclopedia of
  computational mechanics, E. Stein, R. de Borst, and TJ Hughes, Eds., Wiley
  Online Library, 2017.

ChorinOptimalPrediction
Chorin, A. J., Hald, O., and Kupferman, R., Optimal prediction and the
  Mori-Zwanzig representation of irreversible processes, Proc. Natl Acad.
  Sci., Vol. 97, No. (doi:10.1073/pnas.97.7.2968), 2000, pp. 2968–2973.

HaldPredictFromData
Chorin, A. J., Hald, O. H., and Kupferman, R., Prediction from Partial
  Data, Renormalization, and Averaging, Journal of Scientific
  Computing, Vol. 28, No. 2-3, September 2006, pp. 245–261.

GivonOrthogonal
Givon, D., Hald, O. H., and Kupferman, R., Existence Proof for
  Orthogonal Dynamics and the Mori-Zwanzig Formalism, Israel Journal of
  Mathematics, Vol. 145, No. 1, 2005, pp. 221–241.

StatisticalMechanics
Chorin, A. J. and Hald, O. H., Stochastic Tools for Mathematics and
  Science, Springer, 2013.

MoriTransport
Mori, H., Transport, collective motion, and Brownian motion, Prog. Theoret. Phys, Vol. 33, No. 3, 1965, pp. 423–455.

ZwanzigTransport
Zwanzig, R., Memory Effects in Irreversible Thermodynamics, Phys. Rev., Vol. 124, No. 4, 1961, pp. 983.

Darve_CGMZ
Lee, H. S., Ahn, S.-H., and Darve, E. F., Building a Coarse-Grained
  Model Based on the Mori-Zwanzig Formalism, MRS Procedings,
  Vol. 1753, 2015.

Karniadakis_CGMZ
Li, Z., Bian, X., Li, X., and Karniadakis, G. E., Incorporation of
  memory effects in coarse-grained modeling via the Mori-Zwanzig formalism,
  The Journal of Chemical Physics, Vol. 143, No. 24, 2015.

stinisEuler
Hald, O. H. and Stinis, P., Optimal prediction and the rate of decay
  for solutions of the Euler equations in two and three dimensions, Proceedings of the National Academy of Sciences, Vol. 104, No. 16, 2007,
  pp. 6527–6532.

stinisHighOrderEuler
Stinis, P., Higher order Mori-Zwanzig models for the Euler equations,
  Multiscale Model. Simul., Vol. 6, No. 3, 2007, pp. 741–760.

Stinis-rMZ
Stinis, P., Renormalized reduced models for singular PDEs, Commun. Appl. Math. Comput. Sci., Vol. 8, No. 1, 2013, pp. 39–66.

stinis_finitememory
Stinis, P., Mori-Zwanzig reduced models for uncertainty
  quantification I: Parametric uncertainty, arXiv:1211.4285, 2012.

PriceMZ
Price, J. and Stinis, P., Renormalized Reduced Order Models with
  Memory for Long Time Prediction, arXiv:1707.01955v1, 2017.

zhu_mz_faber
Zhu, Y. and Venturi, D., Faber approximation to the Mori-Zwanzig
  equation, arXiv:1708.03806v2, 2017.

zhu_mz_error
Zhu, Y., Dominy, J. M., and Venturi, D., Rigorous error estimates for
  the memory integral in the Mori-Zwanzig formulation, arXiv:1708.02235v2, 2017.

parishAIAA2016
Parish, E. and Duraisamy, K., Reduced Order Modeling of Turbulent
  Flows Using Statistical Coarse-graining, 46th AIAA Fluid Dynamics
  Conference, AIAA AVIATION Forum, Washington, D.C., June 2016.

ParishMZ1
Parish, E. J. and Duraisamy, K., Non-Markovian Closure Models for
  Large Eddy Simulation using the Mori-Zwanzig Formalism, Physical Review
  Fluids, Vol. 2, No. 1, 2017.

Parish_Dtau2
Parish, E. J. and Duraisamy, K., A Dynamic Subgrid Scale Model for
  Large Eddy Simulations Based on the Mori-Zwanzig Formalism, Journal of
  Computational Physics, Vol. 349, No. 15, 2017, pp. 154–175.

bernsteinBurgers
Bernstein, D., Optimal Prediction of Burgers's Equation, Multiscale Model. Simul., Vol. 6, No. 1, 2007, pp. 27–52.

RenormalizedMZ
Stinis, P., Renormalized Mori-Zwanzig-reduced models for systems
  without scale separation, Proceedings of The Royal Society A,
  Vol. 471, No. 20140446, 2015.

GouasmiMZ1
Gouasmi, A., Parish, E. J., and Duraisamy, K., A priori estimation of
  memory effects in reduced-order models of nonlinear systems using the
  Mori-Zwanzig formalism, Proceedings of The Royal Society A,
  Vol. 473, No. 20170385, 2017.

ChandyFrankelLES
Chandy, A. J. and Frankel, S. H., The t-model as a Large Eddy
  Simulation model for the Navier-Stokes Equations, Multiscale Model.
  Simul., Vol. 8, No. 2, 2009, pp. 445–462.

collis_vms
Collis, S., Monitoring unresolved scales in multiscale turbulence
  modeling. Physics of Fluids, Vol. 13, No. 6, 2001, pp. 1800–1806.

hughes_sangalli
Hughes, T. and Sangalli, G., Variational Multiscale Analysis: The
  Fine-Scale Green's Function, Projection, Optimization, Localization, and
  Stabilized Methods, SIAM Journal on Numerical Analysis, Vol. 45,
  No. 2, 2007, pp. 539–557.

adjointFEM
Franca, L., Frey, S., and Hughes, T., Stabilized finite element
  methods: 1. Application to the advective-diffusive model, Computer
  methods in applied mechanics and engineering, Vol. 95, No. 2, 1992,
  pp. 253–276.

Koopman
Koopman, B., Hamiltonian systems and transformations in Hilbert
  space, Proceedings of the National Academy of Sciences of the US A,
  Vol. 17, No. 5, 1931, pp. 315–318.

chorin_book
Chorin, A. and Hald, O., Stochastic Tools in Mathematics and Science,
  Springer-Verlag, 2005.

codina_oss
Codina, R., Stabilization of incompressibility and convection through
  orthogonal sub-scales in finite element methods, Computer methods in
  applied mechanics and engineering, Vol. 190, No. 1579-1599, 2000.

codina_oss0
Codina, R., Stabilized finite element approximation of transient
  incompressible flows using orthogonal subscales, Computer methods in
  applied mechanics and engineering, Vol. 191, No. 4295-4321, 2002.

codina_oss1
Guasch, O. and Codina, R., Statistical behavior of the orthogonal
  subgrid scale stabilization terms in the finite element large eddy simulation
  of turbulent flows, Computer methods in applied mechanics and
  engineering, Vol. 261-262, No. 154-166, 2013.

codina_ossburgers
Bayona, C., Baiges, J., and Codina, R., Variational multi-scale
  approximation of the one-dimensional forced Burgers equation: the role of
  orthogonal sub-grid scales in turbulence modeling, International
  Journal for Numerical Methods in Fluids, 2017.

br1
Bassi, F. and Rebay, S., A High Order Accurate Discontinuous Finite
  Element Method for the Numerical Solution of the Compressible Navier-Stokes
  Equations, Journal of Computational Physics, Vol. 131, 1997,
  pp. 267–279.





