Scale-selective dissipation in energy-conserving finite element schemes for two-dimensional turbulence
aa r X i v : . [ m a t h . NA ] A p r Scale-selective dissipation in energy-conserving finite elementschemes for two-dimensional turbulence
Andrea Natale and Colin J. Cotter
Department of Mathematics, Imperial College London, London, SW7 2AZ, UK
November 5, 2018
Abstract
We analyse the multiscale properties of energy-conserving upwind-stabilised finite elementdiscretisations of the two-dimensional incompressible Euler equations. We focus our attentionon two particular methods: the Lie derivative discretisation introduced in Natale and Cotter(2016a) and the Streamline Upwind/Petrov-Galerkin (SUPG) discretisation of the vorticity ad-vection equation. Such discretisations provide control on enstrophy by modelling different typesof scale interactions. We quantify the performance of the schemes in reproducing the non-localenergy backscatter that characterises two-dimensional turbulent flows.
Numerical simulations are usually unable to fully represent all the dynamically relevant scales inatmospheric flows (Williamson, 2007; Shutts, 2005; Thuburn et al. , 2014). For such unresolvedflows, different discretisations should be compared not only in terms of accuracy and order of con-vergence, but also on the way the unresolved scales affect the computation. More specifically, asboth kinetic energy and enstrophy are statistically relevant in determining large scale flow evolu-tion (Majda and Wang, 2006), their dynamics within the resolved scales should be appropriatelycaptured by the numerical scheme in use.At low resolution, energy conservation is generally considered an appropriate requirement fornumerical schemes modelling atmospheric flows (see, e.g., Shutts, 2005). Enstrophy conservation, onthe other hand, induces a piling up of enstrophy at small scales, a phenomenon known as “spectralblocking”. Such a phenomenon is well known, and it has been studied extensively in the contextof two-dimensional turbulence. The classical picture is described in Maltrud and Vallis (1993),where the different behaviour of energy and enstrophy transfer are studied in terms of triads ofwavenumbers. In particular, one sees that in the enstrophy intertial range the most significanttriads interactions are very nonlocal, with two of the three wavenumbers being much higher thanthe third. The large wavenumbers are dominant when considering enstrophy transfer, which canthen be considered to be local. Energy transfer, on the other hand, also has a non negligiblenonlocal contribution due to the smallest wavenumber, which causes an upscale transfer known asbackscatter.The general approach in numerical schemes is to focus on enstrophy dissipation, to avoid itspiling up at small scales. This however ignores the energy backscatter and yields simulations thatare too dissipative. One possible solution is achieved by combining some form of artificial hyper-viscous dissipation with an “energy fixer” (Shutts, 2005), which reinserts the dissipated energy inthe form of random perturbations or specific large-scale velocity patterns.1huburn et al. (2014) proposed an energy fixer modelled on the instantaneous vorticity dis-tribution and applied it to various finite difference and finite volume schemes for the barotropicvorticity equation, i.e. the two-dimensional incompressible Euler equations, which they consideredas a model problem. Moreover, they compared the behaviour of the various schemes in terms ofenergy and enstrophy tendencies in spectral space. In particular, they quantified how well the dis-cretisations would reproduce the energy and enstrophy tendencies due to small scales in a truncatedDNS simulation.In this paper, we perform a similar analysis, but apply it specifically to finite element upwind-stabilised and energy-conserving discretisations. The reason for our interest is that upwind stabili-sation in finite element methods can be regarded as a form of scale-selective dissipation, providing away of modelling the interaction of scales in the simulation. Depending on the scheme, such an inter-action can be either between resolved and quasi-resolved scales, or between resolved and unresolvedscales (as a proper subgrid model). This property has been used extensively in designing discreti-sations for fluid simulation (see, e.g., Hughes, 1995; Burman, 2007; Becker et al. , 2011), but to theauthors’ knowledge it has never been analysed in the context of energy conserving simulations. Inrecent years, finite element methods have been developed and considered for numerical weather pre-diction and climate modelling (Fournier et al. , 2004; Thomas and Loft, 2005; Dennis et al. , 2011;Kelly and Giraldo, 2012; Giraldo et al. , 2013; Brdar et al. , 2013; Bao et al. , 2015; Marras et al. ,2015). The use of stabilised finite element methods in atmospheric simulations has been investi-gated in (Marras et al. , 2013). Compatible finite element methods have been developed in orderto extend conservation and stability properties from C-grid staggered finite difference methods tothe finite element setting (Cotter and Shipton, 2012; McRae and Cotter, 2014; Natale and Cotter,2016b); these methods provide the context for this paper.Here, we focus our attention on two algorithms that model different scales interactions while pre-serving kinetic energy. The first is the discretisation developed in Natale and Cotter (2016a), whichconsists of a finite element H (div)-conforming scheme with an energy-conserving upwind stabilisa-tion. Such a discretisation was derived using Hamilton’s principle of stationary action applied to anappropriately defined Lagrangian, representing the total kinetic energy. Then, the scheme conservesthe Lagrangian by construction and possesses the attractive feature of a variational derivation, withpossible extensions to different Hamiltonian models. Furthermore, we will observe that the conser-vative upwind stabilisation implicitly induces an energy transfer from small quasi-resolved scales tofully resolved scales, avoiding in this way spectral blocking without dissipating energy. Our secondfocus will be the SUPG stabilised discretisation of the vorticity/stream function formulation of theEuler equations introduced in Tezduyar et al. (1988); Tezduyar (1989). The SUPG stabilisation canbe interpreted as a subgrid model in which the unresolved scales are proportional to the residual ateach time step (Hughes, 1995). Moreover, as the stabilisation is applied to the vorticity advectionequation, only enstrophy is dissipated while kinetic energy is conserved exactly.This paper is organised as follows. In Section 2 we introduce the finite element framework andwe describe the different schemes analysed in the paper. In Section 3 we discuss the multiscale inter-pretation of the upwind stabilisation. In Section 4 we present some numerical tests demonstratingthe behaviour of the schemes in terms of spectral energy and enstrophy tendencies. Conclusions arepresented in Section 5. In this section we introduce some basic finite elements tools (Section 2.1) and describe the methodsstudied in this paper (Section 2.2). In particular, as we treat mixed formulations in velocity/pressure2s well as stream function/vorticity form, we will need to introduce finite element spaces for all thesevariables.
Consider a two-dimensional domain Ω. We denote by T a triangulation of Ω, i.e. a decomposition ofΩ in triangular elements T ⊂ Ω. Let E be the set of element edges of the triangulation T . We fix anorientation for all edges e ∈ E by specifying the unit normal vector n e . An orientation on the edge e ∈ E defines a positive and negative side on e , so that any finite element scalar or vector function f on Ω has two possible restrictions on e , denoted f + and f − , with f + := f | T + and T + ∈ T beingthe element with outward normal n e . Therefore, we can define on each e ∈ E the jump operator[ f ] := f + − f − and the average operator { f } := ( f + + f − ) / T ∈ T , wedenote by P r ( T ) the linear space of polynomials on T of degree up to r . Then the stream functionand vorticity space is given by V = CG r +1 := { ψ : ψ | T ∈P r +1 ( T ) ∀ T ∈ T ,ψ + | e = ψ − | e ∀ e ∈ E} , (1)which is the Continuous Galerkin finite element space of degree r + 1. Clearly, functions in V arecontinuous on Ω, i.e. V ⊂ C (Ω). More precisely, as we will employ periodic boundary conditions,the stream function space will be the space ¯ V of functions in V with zero average on Ω.We now define the velocity and pressure spaces. These are defined in order to produce a stablemixed finite element discretisation, see Boffi et al. (2013) for details. The velocity space V is givenby V = BDM r := { u : u | T ∈ ( P r ( T )) ∀ T ∈ T , u + | e · n e = u − | e · n e ∀ e ∈ E} , (2)which is the Brezzi-Douglas-Marini finite element space of degree r (Boffi et al. , 2013). Note thata vector field in V is not completely continuous, as only the normal components at the edgesare continuous. The subspace of V of continuous velocity fields is denoted by V l , i.e. V l := V ∩ ( C (Ω)) , and it will be used to model the large scale component of the velocity. The quasi-resolved small scale velocity space V s is the orthogonal complement of V l in V , with respect tothe L inner product. In other words, we have V = V l ⊕ V s .The pressure space V is given by V = DG r − := { p : p | T ∈ P r − ( T ) ∀ T ∈ T } , (3)which is the Discontinuous Galerkin finite element space of degree r − et al. , 2013). Moreprecisely, the pressure space will be the space ¯ V of functions in V with zero average on Ω.For any r , the couple ( V , ¯ V ) can be used to produce a stable mixed finite element formulationand enforce the divergence-free constraint exactly on the velocity field. In practice, this constraintsthe velocity to the space ∇ ⊥ V ⊆ V , where the superscript ⊥ denotes a clockwise rotation of π/ V and ¯ V are possible, however, for simplicity, we willrestrict ourself to the spaces mentioned in this section.Finally, for any scalar fields f and g on Ω, we define ( f, g ) Ω := R Ω f · g d x and k f k Ω := p ( f, f ) Ω .The same notation is used for vector and tensor fields, and for integrals over edges or elements.3 .2 Upwind methods for the two-dimensional Euler equations We now introduce two categories of upwind-stabilised finite element schemes for the two-dimensionalEuler equations. The first is based on the velocity/pressure formulation and the second on the streamfunction/vorticity formulation. We start with the former. The Euler equations can be written interms of velocity and pressure as follows, (cid:26) ∂ t u + ( u · ∇ ) u + ∇ p = 0 , ∇ · u = 0 , (4)where u and p are the velocity and pressure fields respectively, and the domain Ω is assumed tobe doubly-periodic. Using the divergence-free condition the first equation in the system (4) can berewritten as ∂ t u + ∇ · ( u ⊗ u ) + ∇ p = 0 , (5)which we refer to as flux form. Alternatively, using standard vector calculus identities, we can alsowrite it as ∂ t u + u ⊥ ( ∇ ⊥ · u ) + ∇ P = 0 , (6)where P := p + k u k /
2, and k · k denotes the Euclidean norm. We refer to Equation (6) asLie derivative formulation, as the advection term in this formulation may be interpreted as the Liederivative of the velocity one-form (Natale and Cotter, 2016a). Then, the general form of the mixedfinite element discretisation of the system in (4) requires us to find ( u , ˜ p ) ∈ V × ¯ V such that (cid:26) ( ∂ t u , v ) Ω + a ( u ; u , v ) + s ( u ; u , v ) − (˜ p, ∇ · v ) Ω = 0 , ( ∇ · u , q ) Ω = 0 , (7)for all ( v , q ) ∈ V × ¯ V , where ˜ p can represent either P or p . For our choice of finite element spaces,the system (7) ensures that ∇ · u = 0 is satisfied pointwise at all times. The forms a ( · ; · , · ) and s ( · ; · , · ) define the discretisation of the advection term, and in particular the second defines theupwind stabilisation. We consider here two particular choices. Flux form discretisation
This is the standard discretisation used in DG schemes for the Navier-Stokes equations (see, e.g., Cockburn et al. , 2004; Guzm´an et al. , 2015, for an application to theEuler equations). Note that such discretisation is not energy-conserving, and therefore we will onlyuse it as a reference due to its similar structure to the Lie derivative discretisation, introduced below.The flux form scheme is derived by writing the advection term in flux form as in Equation (5) andintegrating by parts. It is defined by a ( u ; u , v ) := − X T ∈T ( u , ( u · ∇ ) v ) T + X e ∈E ( u · n e { u } , [ v ]) e , (8) s ( u ; u , v ) := X e ∈E ( c e u · n e [ u ] , [ v ]) e , (9)where c e is a function defined on the edges, and dependent on u . Generally, c e = α u · n e / (2 | u · n e | ),where α > α = 1, the two terms can be combinedby replacing the average in the advection form a with the upwind value of u , i.e. u + if u · n e > u − otherwise. For 0 < α <
1, the average term becomes a skewed average.4ote that, if c e = α u · n e / (2 | u · n e | ) then s ( u ; u , u ) ≥
0, therefore it represents a dissipationterm in the discretisation. Denoting by h the maximum element diameter in the triangulation T , we define τ h := h/ sup Ω ( k u k ) to be the local time scale associated to the mesh. Then, thedissipation induced by the upwinding acts on a time scale bounded from above by τ h /α . Thisfollows immediately from the following string of inequalities, s ( u ; u , u ) ≤ α Ω ( k u k ) X e ∈E ([ u ] , [ u ]) e , ≤ Cα sup Ω ( k u k ) h − k u k = Cατ − h k u k , (10)where C > h . Lie derivative discretisation
We now consider the discretisation introduced in Natale and Cotter(2016a). It can be derived by writing the advection term in Lie derivative form as in Equation (6)and integrating by parts. It is defined by a ( u ; u , v ) := X T ∈T ( u ⊥ , ∇ ( u ⊥ · v )) T − X e ∈E ( { u ⊥ } · n e , [ u ⊥ · v ]) e , (11) s ( u ; u , v ) := − X e ∈E ( c e [ u ⊥ ] · n e , [ u ⊥ · v ]) e . (12)Note that in this case s ( u ; u , u ) = 0 independently of the choice of the function c e . Such a propertydirectly leads to conservation of energy. This can be verified by setting v = u and q = ˜ p in Equa-tion (7), and noting that we also have a ( u ; u , u ) = 0. However, interpreting the stabilisation termas we did for the flux form discretisation is not straightforward. Natale and Cotter (2016a) noticedthat for smooth advecting velocity such a discretisation coincides with the Eulerian Lie derivativediscretisation proposed in Heumann and Hiptmair (2011) for the linear advection-diffusion prob-lem. In Section 3 we elaborate on this observation, and we give a multiscale interpretation to thestabilisation term.Numerical tests show that both the flux form and the Lie derivative scheme converge withorder r + 1, in terms of velocity L error, when standard upwinding is employed, i.e. α = 1,see Guzm´an et al. (2015) and Natale and Cotter (2016a). A priori convergence estimates withsuboptimal convergence rate r can be found in the same references.We now consider the stream function/vorticity formulation of the two-dimensional Euler equa-tions which is the basis for the SUPG discretisation. This is given by (cid:26) ∂ t ω + u · ∇ ω = 0 , ∆ ψ = ω, (13)where u := ∇ ⊥ ψ . The general form of the mixed finite element discretisation of the system in (4)is given by : find ( ψ, ω ) ∈ ¯ V × V such that (cid:26) ( ∂ t ω, ϕ ) Ω + ˜ a ( u ; ω, ϕ ) + ˜ s ( u ; ω, ϕ ) = 0 , ( ∇ ψ, ∇ φ ) Ω = − ( ω, φ ) Ω , (14)for all ( ϕ, φ ) ∈ V × ¯ V . 5 UPG discretisation
The SUPG discretisation is based on the standard Galerkin discretisationfor vorticity advection, supplemented with a stabilisation term which is taken to be proportionalto the residual of the vorticity advection equation (Tezduyar et al. , 1988; Tezduyar, 1989), i.e. wechoose ˜ a ( u ; ω, ϕ ) := ( u · ∇ ω, ϕ ) Ω , (15)˜ s ( u ; ω, ϕ ) := ( τ s R ( ω ) , u · ∇ ϕ ) Ω , (16)where R ( ω ) := ∂ t ω + u · ∇ ω and where τ s ≥ a ( u ; ω, ψ ) = ˜ s ( u ; ω, ψ ) = 0 and setting ϕ = ψ in Equation (14). The stabilisation termis trivially consistent being proportional to the residual. Moreover it induces streamline dissipationon the enstrophy density. The factor τ s defines the time scale at which the dissipation acts and itis referred to as “intrinsic time scale” (Hughes, 1995). Generally speaking, it can be taken to beproportional to the global value τ h . However, defining τ s = τ h may lead to over-dissipative solutions,especially for higher order elements (Almeida and Silva, 1997; Codina et al. , 1992). For the linearadvection problem with high order elements, Almeida and Silva (1997) suggest the value τ s = βh T ξ k u k (17)where h T is the characteristic length of the element T ∈ T , depending on the shape and the localvelocity direction of the element, ξ = 1 for linear elements and ξ = 1 / β is a constant; in Almeida and Silva (1997) β = 1, but we will also consider differentvalues in our numerical tests.We note that the stream function/vorticity SUPG scheme can be re-interpreted as a mixedvelocity-pressure scheme, since if ψ ∈ V , then ∇ ⊥ ψ spans the divergence-free subspace of V , andin the absence of boundaries, ω ∈ V can be obtained by solving( γ, ω ) Ω = − ( ∇ ⊥ γ, u ) Ω , (18)for all test functions γ in V . Then we obtain an equivalent formulation that is an approximationof the incompressible Euler equations in Lie derivative form, with a ( u ; u , v ) =( v , u ⊥ ω ) Ω , (19) s ( u ; u , v ) =( v , − u ⊥ τ s R ( ω )) Ω . (20)This formulation allows us to extend the energy-conserving SUPG approach to compatible finite ele-ment method discretisations of the shallow water equations in velocity-height formulation, followingMcRae and Cotter (2014). The equivalent discretisation is( ∂ t u , v ) Ω + ( v , F ⊥ q ) Ω +( v , − F ⊥ τ s ( q t + u · ∇ q )) Ω − (cid:18) ∇ · v , gh + 12 | u | (cid:19) Ω = 0 , (21)( ∂ t h + ∇ · F , φ ) Ω = 0 , (22)( F , w ) Ω − ( u h, w ) Ω = 0 , (23)( γ, qh ) Ω + ( ∇ ⊥ γ, u ) Ω − ( γ, f ) Ω = 0 , (24)for all test functions ( v , φ, w , γ ) ∈ V × V × V × V , having introduced potential vorticity q ∈ V andmass flux F ∈ V , and where g is the acceleration due to gravity and f is the Coriolis parameter.6 Multiscale interpretation
We now describe the multiscale character of the Lie derivative (Section 3.1) and the SUPG discreti-sations (Section 3.2). We show that the Lie derivative scheme models the interaction of resolved andquasi-resolved scales within the velocity space V , whereas the SUPG scheme models the interactionof the resolved and unresolved scales in the vorticity advection equation. Consider the Lie derivative scheme and recall the decomposition V = V l ⊕ V s defined in Section 2.1. Let u = u l + u s , where u l and u s are the L projections of u onto V l and V s respectively. Then,taking v = u l in the expression for s in Equation (12), we obtain s ( u ; u , u l ) = − X e ∈E ( c e [ u ⊥ ] · n e , [ u ⊥ ] · u l ) e = − X e ∈E ( c e u l · n e [ u ] , [ u ]) e , (25)where we used the fact that u l is continuous and standard vector calculus identities. Similarly,taking v = u s yields s ( u ; u , u s ) = X e ∈E ( c e u l · n e [ u ] , [ u ]) e . (26)Due to the L orthogonality of u l and u s , by taking v = u l and v = u s in Equation (4) weobtain the evolution equations for the large scale and small scale kinetic energy respectively, i.e. E l := k u l k / E s := k u s k /
2. These are given by d E l d t + a ( u ; u , u l ) − ( p, ∇ · u l ) Ω = X e ∈E ( c e u l · n e [ u ] , [ u ]) e d E s d t + a ( u ; u , u s ) − ( p, ∇ · u s ) Ω = − X e ∈E ( c e u l · n e [ u ] , [ u ]) e (27)Then, if we set c e = α u l · n e / (2 | u l · n e | ) (by analogy to the usual choice in the flux form discreti-sation), we introduce an artificial energy transfer from small to large scales. This is because, forsuch a choice, s ( u ; u , u l ) ≤ s ( u ; u , u s ) ≥
0. Moreover, the total energy is conserved since s ( u ; u , u l ) + s ( u ; u , u s ) = s ( u ; u , u ) = 0.Finally, note that once the problem is reformulated as in Equation (27), the considerationsof the previous section on the dissipation time scale apply without major changes. In particular,the dissipation time scale in the small scale equation is bounded from above by τ h /α , with τ h := h/ (sup Ω k u l k ). Consider now the SUPG scheme. We follow the variational multiscale (VMS) framework, and inparticular Codina (2011), to give a simplified interpretation of the SUPG stabilisation as a subgrid7cale model. We assume that the exact vorticity ω ∈ V , where V is a sufficiently regular space, sothat the following variational formulation, (cid:26) ( ∂ t ω, ϕ ) Ω + ( u · ∇ ω, ϕ ) Ω = 0 , u = ∇ ⊥ (∆ − ω ) , (28)holds for all ϕ ∈ V . Assume that V = V ⊕ V u , where V u contains the information from theunresolved scales and it is not necessarily orthogonal in L to V . Then we can decompose ω =¯ ω + ω u , and by linearity u = ¯ u + u u . Using these decompositions, the vorticity equation becomes ( ∂ t ¯ ω + ∂ t ω u + ¯ u · ∇ ¯ ω + ¯ u · ∇ ω u , ¯ ϕ ) Ω +( u u · ∇ ¯ ω + u u · ∇ ω u , ¯ ϕ ) Ω = 0 , ( ∂ t ¯ ω + ∂ t ω u + ¯ u · ∇ ¯ ω + ¯ u · ∇ ω u , ϕ u ) Ω +( u u · ∇ ¯ ω + u u · ∇ ω u , ϕ u ) Ω = 0 , (29)for all ( ¯ ϕ, ϕ u ) ∈ V × V u . We now assume that the dynamics of the resolved and unresolved scalesis driven only by the resolved velocity component. This assumption can be justified heuristicallyas in Laval et al. (1999). If we also neglect the unsteady part of the unresolved component of thevorticity, we are left with the following system (cid:26) ( ∂ t ¯ ω, ¯ ϕ ) Ω + ( ¯ u · ∇ ¯ ω, ¯ ϕ ) Ω − ( ω u , ¯ u · ∇ ¯ ϕ ) Ω = 0( ∂ t ¯ ω, ϕ u ) Ω + ( ¯ u · ∇ ¯ ω + ¯ u · ∇ ω u , ϕ u ) Ω = 0 (30)for all ( ¯ ϕ, ϕ u ) ∈ V × V u . Finally, we approximate the unresolved scales equation by taking ( ¯ u ·∇ ω u , ϕ u ) ≈ τ − s ( ω u , ϕ u ). Then, a possible solution of the unresolved scales equation is ω u = τ s R (¯ ω ),with R (¯ ω ) = ∂ t ¯ ω + ¯ u · ∇ ¯ ω . The resulting scheme coincides with the SUPG scheme, with globalintrinsic time scale τ s . In fact, this derivation shows that for our problem the VMS and SUPGmethods are equivalent. This is due to the absence of the dissipation term in the equations studiedhere. We remark that different versions of the above derivation have been used to devise similaralgorithms. In Hughes (1995), for example, a different analysis leads to the definition of a local timescale, which originates from the Green’s function associated to the unresolved scales equation. InCodina and Blasco (2002) the time dependent term in the unresolved scales equation is retained,leading to a scheme in which the unresolved scales are dynamically tracked in time. These schemesstill conserve energy as they share the same structure, however, in this paper we limit ourself tostudy the simpler SUPG scheme.It should be noted that the SUPG scheme models enstrophy transfer between different scalesrather than energy transfer, as it was the case for the Lie derivative scheme. This is only possibledue to the fact that we have an explicit equation for the vorticity evolution. Then, the diffu-sion introduced by the stabilisation term only affects enstrophy, and does not compromise energyconservation. In this section we present some numerical results for the discretisations discussed in the previoussections. We start by reporting on their order of convergence using a manufactured solution test(Section 4.1). Then, we perform a forced turbulence simulation measuring energy and enstrophytendencies in spectral space, in order to quantify the multiscale behaviour of the schemes (Section 4.2and 4.3). We use as a benchmark test case the forced turbulence test performed in Thuburn et al. (2014), and we can compare our results to their reference solution and to the other schemes tested8 h SUPG Lie derivativeerror order error order1 7.40e-1 2.68e-2 – 4.12e-2 –3.70e-1 6.69e-3 2.00 5.38e-3 2.192.47-1 2.97e-3 2.00 9.00e-3 2.101.85e-1 1.67e-3 2.00 3.84e-3 2.092 7.40e-1 1.37e-3 – 3.03e-3 –3.70e-1 1.71e-4 3.01 3.20e-4 3.252.47-1 5.05e-5 3.00 8.94e-5 3.141.85e-1 2.13e-5 3.00 3.63e-5 3.13Table 1: Comparison between the SUPG and Lie derivative scheme in terms of the error k u − u h k Ω and the order of convergence, for the solution in Equation 31, and for α = β = 1.therein. Note that as in Thuburn et al. (2014) for all the tests of this section we regard the equationsof motion as non-dimensional, since they can always be brought in this form via an appropriate timerescaling. All computations were performed using the Firedrake software suite (Rathgeber et al. ,2015), which allows for symbolic implementation of finite element problems of mixed type. We alsoused additional resources for the implementation from the following references (Balay et al. , 2015,1997; Dalcin et al. , 2011; Hendrickson and Leland, 1995; Amestoy et al. , 2001, 2006). The method of manufactured solutions provides a useful way to investigate convergence rates whenanalytic solutions are not available. This consists in choosing an arbitrary function and adding anappropriate forcing to the system of equations we want to solve, so that this function is an exactsolution of the modified system. We let Ω = [0 , × [0 , x, y ), andpick as manufactured solution the one generated by the following stream function, ψ ( t, x, y ) = sin( πx ) sin( πy ) sin( πy − t ) e − t/σ , (31)with σ = 100, t ∈ [0 ,
1] and boundary conditions ψ = 0 on ∂ Ω. The L errors in velocity andvorticity are reported in Table 1 and 2 respectively. For the Lie derivative scheme we set α = 1,whereas for the SUPG scheme the stabilisation coefficient is given by Equation (17) with β = 1.Both schemes are integrated in time using the implicit midpoint integration rule with ∆ t = 10 − .Note that for the Lie derivative scheme the vorticity needs to be computed weakly at each time stepsince it is not used directly by the algorithm. The schemes behave similarly in terms of velocity,where r + 1 order of convergence is observed. However, the SUPG scheme yields lower error andconvergence rate greater than r + 1 for the vorticity, whereas the Lie derivative scheme only givesapproximately order r convergence. We now describe the reference turbulent solution considered in Thuburn et al. (2014). We start bymodifying the governing equation by adding a forcing and a dissipation term as follows: ∂ t ω + u · ∇ ω = f − ω/τ , (32)9 h SUPG Lie derivativeerror order error order1 7.40e-1 1.04e-1 – 1.58 –3.70e-1 2.16e-2 2.27 8.05e-1 0.972.47-1 8.62e-3 2.26 5.51e-1 0.931.85e-1 4.48e-3 2.28 4.21e-1 0.942 7.40e-1 3.50e-3 – 1.34e-1 –3.70e-1 2.79e-4 3.65 2.69e-2 2.322.47-1 6.74e-5 3.50 1.10e-2 2.211.85e-1 2.51e-5 3.43 5.79e-3 2.22Table 2: Comparison between the SUPG and Lie derivative scheme in terms of the error k ω − ω h k Ω and the order of convergence, for the solution in Equation 31, and for α = β = 1.on a periodic unit square domain Ω, with Cartesian coordinates ( x, y ). The fixed forcing is f =0 . πx ) and acts on the wave number k = 16. The dissipation term is the standard Ekmandrag with time scale τ = 100, which does not introduce selective decay (Majda and Wang, 2006).The initial conditions are ω | t =0 = sin(8 πx ) sin(8 πy )+0 . πx ) cos(6 πy )+0 . πx ) cos(4 πy )+0 .
01 sin(2 πy ) + 0 .
02 sin(2 πx ) , (33)after which the flow rapidly becomes chaotic.Thuburn et al. (2014) analysed the stationary state in terms of energy and enstrophy tendenciesin spectral space, with the aim of quantifying the effect of small scales (above a certain threshold) onlarge scales. In order to compute the reference values, they advanced the solution of Equation (32)until t = 200 using a spectral method with high resolution, computing the stream function ψ andthe vorticity ω . Then, for t ∈ [200 , J = ∇ ⊥ ψ · ∇ ω and then theenergy and enstrophy time derivative at each time step, using the following formulas˙ E ( k x , k y ) = Re ( ˆ ψ ∗ ˆ JN ) , ˙ Z ( k x , k y ) = Re ( ˆ ω ∗ ˆ JN ) , (34)where N is the number of grid points, ˆ ψ and ˆ J are the Fourier transforms of ψ and J (truncatedat the maximum retained wavenumber), and the superscript ∗ denotes complex conjugation. Theresults were integrated over angle in spectral space to give ˙ E ( k ) and ˙ Z ( k ), where k := q k x + k y .The same procedure was repeated truncating all the Fourier transforms with k ≤ k T , for severalvalues of k T , yielding ˙ E T ( k ) and ˙ Z T ( k ). Finally, the energy and enstrophy tendencies due to scales k > k T were computed as follows:˙ E SG ( k ) = ˙ E ( k ) − ˙ E T ( k ) , ˙ Z SG ( k ) = ˙ Z ( k ) − ˙ Z T ( k ) . (35)The results of such calculations are illustrated in Figure 1. The main features are a net non-local energy transfer from small to large scales and enstrophy dissipation concentrated close tothe truncation wavenumber. Reproducing this behaviour on the resolved scale of a numericalsimulation is very challenging. In the following, we review the behaviour of the schemes presentedin the previous section. 10
10 20 30 40 50 60 − T = 42Wavenumber × − T = 42Wavenumber × . − − T = 85Wavenumber × − − T = 85Wavenumber × . − − T = 170Wavenumber × − − T = 170Wavenumber × . (a)(b)(c) Figure 1: Energy and enstrophy tendencies for the reference solution, with k T = 42 (top), k T = 85(middle), k T = 170 (bottom) (Reproduced from Thuburn et al. (2014)).11 .3 Finite element turbulence tests We now repeat the experiment described in the last section, but using the methods defined inSection 2.2. In practice, we substitute ˙ E ( k ) and ˙ Z ( k ) in Equation (35) with those computed usingthe numerical Jacobian implied by the discretisations. Instead of restarting each simulation forthe initial condition specified by a DNS solution at t = 200 as in Thuburn et al. (2014), we startmonitoring the several schemes at a later time t = 300, so that each simulation can be consideredto have reached its own statistical equilibrium. In particular, we verified this by checking the timeevolution of the energy moving average. For all simulations we used a structured right triangulargrid and we denote by h the horizontal and vertical spacing of the nodes. Moreover, they are allintegrated in time using the implicit midpoint integration rule with ∆ t = 2 · − until t = 350.We compare the results for the flux form scheme, the Lie derivative scheme and the SUPG schemewhen using the same function spaces. The results for the flux form scheme are mainly reportedfor comparison with the Lie derivative scheme, since these schemes share a similar structure andthe same convergence rate. In all tests, the stabilisation factor for the SUPG scheme is given byEquation (17) with h T = h/ √
2, as suggested by Codina et al. (1992) for a symmetrical grid of righttriangles.In Figure 2 results are shown for r = 1, h = 1 / α = 1 and β = 1. The mesh resolution issuch that the maximum resolved wavenumber (computed using the grid associated to the streamfunction degrees of freedom) is approximately k = 85. For this setting we show the effect of thenumerical scheme on the integration for the cut-off wavenumbers tested in Thuburn et al. (2014).We start by analysing the energy behaviour (left column in Figure 2). At the lowest cut-offwavenumber k T = 42, the positive energy peak is overestimated for the Lie derivative scheme butit is completely absent when using the flux form discretisation. Both schemes are too dissipative atintermediate scales, with the flux form discretisation being more dissipative at the forcing scale. For k T = 85 the behaviour of the two schemes is substantially different. The flux form scheme appearsto subtract energy over the whole spectrum, and especially at the most energetic wavenumbers(since energy is inserted at k = 16). The Lie derivative scheme also subtracts energy over alarger wavenumber interval compared to the reference test in Figure 1, but energy is reinserted atlow frequencies, as expected from the discussion in the previous section. However, note that thepositive peak is considerably larger than the one in the reference solution.The SUPG scheme shows a much better agreement to the reference test, at equal resolution.Note, in particular, that for k T = 85 we see an energy dissipation peak close to the cut-off wavenum-ber, which was absent for the other schemes tested. This is probably due to the fact that the vorticityis explicitly used in the computation of the Jacobian, leading to a lower error in the computationof the energy and enstrophy tendencies. Such a conjecture is corroborated by the order of con-vergence results given in Section 4.1. Nonetheless, most of the dissipation is still acting on lowerwavenumbers, ranging almost down to the forcing scale.As for the enstrophy behaviour (right column in Figure 2), the Lie derivative and the flux formschemes show similar behaviour for k T = 42, where enstrophy is subtracted mainly at wavenumbersclose to the cut-off wavenumber as in Figure 1. For k T = 85 however, even if the negative peakat the highest wavenumber is still represented, this is much smaller than in the reference solution,and enstrophy appears to be dissipated comparatively more over the rest of the spectrum for bothschemes. Again, the SUPG scheme produces a much better agreement to the reference solution,and it is able to capture the sharp dissipation enstrophy peak also for k T = 85.In Figure 3 the same test is repeated at higher resolution h = 1 /
256 and same polynomial degree r = 1, with maximum retained wavenumber approximately equal to k = 170. For k T = 42, theresolved scales dominate the energy and enstrophy tendencies for all the schemes, so that we can12learly see the non-local energy backscatter correctly represented. However, the energy conservingschemes show again a larger positive peak at low frequencies. For k T = 85, we can see that theLie derivative scheme produces a positive energy peak, which is closer to its reference value, ifcompared to the coarser resolution test in Figure 2. The peak is absent in the results for flux formdiscretisation. Nonetheless, the plot for k T = 170 confirms that the energy transfer is still stronglylocal as dissipation is mostly acting close to the forcing scale. The SUPG scheme, on the otherhand, does not overestimate the energy peak for k T = 85, and produces a much more nonlocalenergy transfer, as it can be verified from the plots for k T = 170.Next, in Figure 4, we performed another test with large mesh size, i.e. we se set h = 1 / r = 2. For this case the number of degrees of freedom is thesame as the test case in Figure 2. The behaviour, however, is substantially better for the Liederivative scheme, with energy being dissipated over larger wavenumbers and the positive peakbeing comparable to the higher resolution tests in Figure 3. The improvement is more modest forthe SUPG scheme and does not compare to the increased accuracy achieved with mesh refinement.It should be noted that for both the Lie derivative and the SUPG schemes the upscale energytransfer is often accompanied by a spurious enstrophy injection at large scales. This effect isparticularly pronounced in the low resolution test in Figure 2. In order to obtain an indication ofthe consequences of such a behaviour we compare the vorticity fields for the different schemes inphysical space at fixed time t = 300, see Figures 6, 7 and 8. We observe that, especially at lowresolution ( h = 1 / r = 1), the vorticity field for the Lie derivative scheme appears smootherthan the one obtained using the flux form scheme. Their behaviour however is comparable, whereasthe SUPG scheme provides a more accurate representation of the vorticity field, with smaller andmore dense vortices. This reflects the higher accuracy of the SUPG discretisation, and it can alsobe linked to the unphysical enstrophy dissipation at intermediate scales which affects both the Liederivative and the flux form scheme. We also note that, when higher polynomial order is used( h = 1 / r = 2), the SUPG and Lie derivative scheme yield overshoots in the vorticity intensity.The effect is stronger for the SUPG scheme, although the overall behaviour of the solution is stillwell-captured.To complete the analysis, we examine the dependence of the results from the stabilisationparameters α and β , which can be used to regulate the level of dissipation. In Figure 5 we consider α, β ∈ { . , . , . } , h = 1 /
128 and r = 1. Remarkably, reducing α in the Lie derivative schemedoes not affect significantly the positive peak in the energy tendency plot. On the other hand,enstrophy dissipation moves to smaller scales. This, however, is due to the appearance of grid scaleoscillations, which can be noticed in the vorticity distribution at the final time t = 350 already for α = 0 .
5. Therefore, such behaviour is a signal of instability rather than improved accuracy. Similarconsiderations hold for the SUPG scheme. Note that for the SUPG scheme both the energy andenstrophy tendencies plots are not much affected by reducing β . However, once again, at the finaltime t = 350 and for β = 0 .
25 we observed gridscale oscillations appear in the vorticity field.
We now remove dissipation and forcing from the equations and compare the energy and enstrophyevolultion for the different schemes. We set as initial condition the one specified in Equation (33)and evolve the solution until t = 100, with α = 1 and β = 1. The results of this test are collectedin Figure 9.We observe that the Lie derivative and SUPG schemes conserve energy up to machine precisionas expected. The flux form scheme dissipates energy, although this behaviour is attenuated withrefinement. The enstrophy evolution for the Lie derivative and flux form scheme are very similar,13nd the largest differences can be seen at the lowest resolution test, i.e. for h = 1 / r = 1,where the Lie derivative scheme is slightly less dissipative. On the other hand, the SUPG schemedissipates substantially less enstrophy than the other two schemes. Note, in particular, the artificialdrop in enstrophy at the beginning of the simulation for the Lie derivative and flux form schemewhen r = 1, which is absent from the SUPG results. We analysed the Lie derivative finite element discretisation of the incompressible Euler equationsintroduced in Natale and Cotter (2016a) and the SUPG discretisation of the vorticity advectionequation, in terms of energy and enstrophy tendencies in a forced turbulence test case. The mainpoints of this paper can be summarised as follows: • The Lie derivative upwind discretisation in Natale and Cotter (2016a) can be reformulatedto model in a simple way energy backscatter from small to large scales, by an appropriatedecomposition of the velocity finite element space. The SUPG scheme, on the other hand,can be interpreted as a way to model enstrophy transfer towards unresolved scales. • Energy conservation in both schemes does not rely on an a priori choice of a vorticity per-turbation pattern, as is usually the case for methods based on energy fixers (Thuburn et al. ,2014). • The Lie derivative scheme can be generalised to different models. This is because, on onehand, it is not exclusively a two-dimensional method; on the other, as the method descendsfrom a Lagrangian formulation, one can derive similar algorithms by appropriately definingdifferent energy functionals (Natale and Cotter, 2016a). • Due to the nature of upwind stabilisation, the energy transfer is still strongly local for bothschemes, however the behaviour is comparable to that obtained with standard methods com-bined with energy fixers (Thuburn et al. , 2014). • The SUPG scheme shows a similar qualitative behaviour to the Lie derivative discretisationbut substantially better agreement with the reference solutions, probably due to higher orderrepresentation of the Jacobian. 14igure 2: Energy and enstrophy tendencies for the Lie derivative (a), the flux form (b) and theSUPG (c) scheme (with h = 1 / r = 1, α = 1, β = 1) on selected spectral intervals k ≤ k T .15igure 3: Energy and enstrophy tendencies for the Lie derivative (a), the flux form (b) and theSUPG (c) scheme (with h = 1 / r = 1, α = 1, β = 1) on selected spectral intervals k ≤ k T .16igure 4: Energy and enstrophy tendencies for the Lie derivative (a), the flux form (b) and theSUPG (c) scheme (with h = 1 / r = 2, α = 1, β = 1) on selected spectral intervals k ≤ k T .17igure 5: Energy and enstrophy tendencies for the Lie derivative and SUPG scheme (with h = 1 / r = 1) for different values of the stabilisation parameters α and β , on the spectral interval k ≤ a) Flux form (b) Lie derivative (c) SUPG Figure 6: Vorticity field at t = 300, for the tested schemes (with h = 1 / r = 1, α = 1, β = 1). (a) Flux form (b) Lie derivative (c) SUPG Figure 7: Vorticity field at t = 300, for the tested schemes (with h = 1 / r = 1, α = 1, β = 1). (a) Flux form (b) Lie derivative (c) SUPG Figure 8: Vorticity field at t = 300, for the tested schemes (with h = 1 / r = 2, α = 1, β = 1).19igure 9: Energy and enstrophy history for the Lie derivative (a), the flux form (b) and the SUPG(c) scheme (with α = 1, β = 1). 20 eferences Almeida, R. C. and Silva, R. S. (1997). A stable Petrov-Galerkin method for convection-dominatedproblems.
Computer Methods in Applied Mechanics and Engineering , SIAM Journal on Matrix Analysis andApplications , Parallel Computing , Modern Software Tools in Scientific Computing , pages 163–202. Birkh¨auserPress.Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout,V., Gropp, W. D., Kaushik, D., Knepley, M. G., McInnes, L. C., Rupp, K., Smith, B. F., Zampini,S., and Zhang, H. (2015). PETSc users manual. Technical Report ANL-95/11 - Revision 3.6,Argonne National Laboratory.Bao, L., Kl¨ofkorn, R., and Nair, R. D. (2015). Horizontally explicit and vertically implicit (HEVI)time discretization scheme for a discontinuous Galerkin nonhydrostatic model.
Monthly WeatherReview , ComptesRendus Mathematique , Mixed finite element methods and applications . Springerseries in computational mathematics. Springer, Berlin, Heidelberg.Brdar, S., Baldauf, M., Dedner, A., and Kl¨ofkorn, R. (2013). Comparison of dynamical cores forNWP models: comparison of COSMO and Dune.
Theoretical and Computational Fluid Dynamics , Computer Methods in Applied Mechanics andEngineering , Mathematics of Computation , pages 1067–1095.Codina, R. (2011).
Finite Element Approximation of the Convection-Diffusion Equation: Subgrid-Scale Spaces, Local Instabilities and Anisotropic Space-Time Discretizations , pages 85–97.Springer Berlin Heidelberg, Berlin, Heidelberg.Codina, R. and Blasco, J. (2002). Analysis of a stabilized finite element approximation of thetransient convection-diffusion-reaction equation using orthogonal subscales.
Computing and Vi-sualization in Science , (3), 167–174.Codina, R., O˜nate, E., and Cervera, M. (1992). The intrinsic time for the streamline upwind/Petrov-Galerkin formulation using quadratic elements. Computer Methods in Applied Mechanics andEngineering , Journalof Computational Physics , Advances in Water Resources , International Journal of High Performance ComputingApplications , page 1094342011428142.Fournier, A., Taylor, M. A., and Tribbia, J. J. (2004). The spectral element atmosphere model(SEAM): High-resolution parallel computation and localized resolution of regional dynamics.
Monthly Weather Review , SIAM Journal onScientific Computing , Preprint .Hendrickson, B. and Leland, R. (1995). A multilevel algorithm for partitioning graphs. In
Super-computing ’95: Proceedings of the 1995 ACM/IEEE Conference on Supercomputing (CDROM) ,page 28, New York. ACM Press.Heumann, H. and Hiptmair, R. (2011). Eulerian and semi-Lagrangian methods for advection-diffusion for differential forms.
Discrete and Continuous Dynamical Systems , ComputerMethods in Applied Mechanics and Engineering , Journal ofComputational Physics , Phys. Rev. Lett. ,
3, 4061–4064.Majda, A. and Wang, X. (2006).
Non-Linear dynamics and statistical theories for basic geophysicalflows . Cambridge University Press.Maltrud, M. E. and Vallis, G. K. (1993). Energy and enstrophy transfer in numerical simulations oftwo-dimensional turbulence.
Physics of Fluids A: Fluid Dynamics (1989-1993) , (7), 1760–1775.Marras, S., Moragues, M., V´azquez, M., Jorba, O., and Houzeaux, G. (2013). Simulations of moistconvection by a variational multiscale stabilized finite element method. Journal of ComputationalPhysics ,
52, 195–218. 22arras, S., Kelly, J. F., Moragues, M., M¨uller, A., Kopera, M. A., V´azquez, M., Giraldo, F. X.,Houzeaux, G., and Jorba, O. (2015). A review of element-based Galerkin methods for numericalweather prediction: Finite elements, spectral elements, and discontinuous Galerkin.
Archives ofComputational Methods in Engineering , pages 1–50.McRae, A. T. T. (2015).
Compatible finite element methods for atmospheric dynamical cores . Ph.D.thesis, Imperial College London.McRae, A. T. T. and Cotter, C. J. (2014). Energy-and enstrophy-conserving schemes for the shallow-water equations, based on mimetic finite elements.
Quarterly Journal of the Royal MeteorologicalSociety , Submitted to IMA Journal of Numerical Analysis .Natale, A. and Cotter, C. J. (2016b). Compatible finite element spaces for geophysical fluid dy-namics. arXiv preprint arXiv:1605.00551 .Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., McRae, A. T., Bercea, G.-T.,Markall, G. R., and Kelly, P. H. (2015). Firedrake: automating the finite element method bycomposing abstractions.
Submitted to ACM TOMS .Shutts, G. (2005). A kinetic energy backscatter algorithm for use in ensemble prediction systems.
Quarterly Journal of the Royal Meteorological Society , Computer Methods in AppliedMechanics and Engineering ,
3, 331–339.Tezduyar, T., Glowinski, R., and Liou, J. (1988). Petrov-Galerkin methods on multiply connecteddomains for the vorticity-stream function formulation of the incompressible Navier-Stokes equa-tions.
International journal for numerical methods in fluids , (10), 1269–1290.Thomas, S. J. and Loft, R. D. (2005). The NCAR spectral element climate dynamical core: Semi-implicit Eulerian formulation. Journal of Scientific Computing , Quarterly Journal of the Royal Meteorological Society , Journal of the Meteorological Society of Japan. Ser. II ,8