Robust error estimates in weak norms for advection dominated transport problems with rough data
NNovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100
Mathematical Models and Methods in Applied Sciencesc (cid:13)
World Scientific Publishing Company
Robust error estimates in weak norms for advection dominatedtransport problems with rough data
Erik Burman
Department of Mathematics,University College London, Gower Street, London,UK–WC1E 6BT,United [email protected]
Received (Day Month Year)Revised (Day Month Year)Communicated by (xxxxxxxxxx)We consider transient convection–diffusion equations with a velocity vector field withmultiscale character and rough data. We assume that the velocity field has two scales, acoarse scale with slow spatial variation, which is responsible for advective transport anda fine scale with small amplitude that contributes to the mixing. For this problem weconsider the estimation of filtered error quantities for solutions computed using a finiteelement method with symmetric stabilization. A posteriori error estimates and a priorierror estimates are derived using the multiscale decomposition of the advective velocityto improve stability. All estimates are independent both of the P´eclet number and of theregularity of the exact solution.
Keywords : passive transport; convection–diffusion ; stabilized finite element method ;error estimates.AMS Subject Classification: 65M12, 65M60, 65M20, 65M15
1. Introduction
In spite of much progress in recent years the problem of deriving a posteriori errorestimates and a priori error estimates for transient convection–diffusion equations isnot fully understood. The difficulty is related to the wide range of different problemscovered by the equation class. Indeed depending on the characteristics of the velocityfield and the molecular diffusion the solutions may feature very different behaviour.From a computational point of view the complexity will depend mainly on thesmoothness of data and on the P´eclet numberPe L := U Lµ , where U is a characteristic velocity, L is a lengthscale and µ the molecular diffu-sion. If the variations of the transport velocity are small and the data are smooth,the solution will be smooth, with moderate Sobolev norm (at least the H -norm) a r X i v : . [ m a t h . NA ] M a y ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman independent of the viscosity. Then a standard Galerkin method for low P´eclet num-ber flows and a stabilized finite element method for high P´eclet number flows willyield accurate results. Most difficult is the case of a high P´eclet number and astrongly varying, or even turbulent, velocity field, transporting a concentration thatis strongly fragmented and may dilute or concentrate. In this case a computationcan experience strong amplification of errors due to repeated bifurcation of stream-lines and inexact representation of internal layers caused by spurious oscillationsor numerical diffusion. Sometimes this regime is referred to as scalar turbulenceand LES-models have been derived for the modelization of the passive scalar us-ing filtering . Such models encounter a similar Reynolds stress conundrum as thefiltering of the full Navier-Stokes’ equations, and hence the modeling error is dif-ficult to quantify. Another approach that has been attempted for this problem isheterogeneous multiscale methods however in that case the underlying theory isbased on homogenization and depending on a periodicity hypothesis that in mostapplications will not hold.In this paper our approach is to apply a stabilized finite element method tothe computation of the solution of the standard physical model, instead of a coarsegrained model. The accuracy of the large scales is measured by estimating theregularized, or filtered error, related to estimating the error in local averages of thesolution.The combination of these two ingredients allows us to derive error estimateswith an order in h , for a norm that is in a certain sense in between H − and L ,but which contains the L -norm of a filtered error. These estimates are robust in thesense that they do not depend on any high order Sobolev norm of the exact solutionand they only have exponential growth depending on the maximum gradient of thecoarse scale velocity field, under a certain scale separation assumption given below.This means that the filtered quantities considered are robust under diverging finescale characteristic trajectories. In some sense we extract the coarse scales for whichwe have some (provable) accuracy from the computation.
2. The transient advection–diffusion equation
The problem that we will consider takes the following form. Let Ω be an open, convexpolygonal/polyhedral subset of R d , d = 2 , ∂ Ω and associatedoutward pointing normal n ∂ Ω . We will denote the computational time interval by I := [0 , T ] and the space-time domain by Q := I × Ω. Assuming homogeneousDirichlet boundary conditions u | ∂ Ω = 0 we may formally write our problem, for t > u ∈ H (Ω) such that u ( x,
0) = u ( x ) in Ω and ∂ t u + β · ∇ u − µ ∆ u = f, in Ω . (2.1)For the problem data we consider u ∈ L (Ω), f ∈ L ( Q ) and let the velocity field β ∈ [ C (0 , T ; W , ∞ (Ω))] d , such that ∇ · β = 0, β · n ∂ Ω | ∂ Ω = 0, and the moleculardiffusivity µ ∈ R with µ >
0. The L -scalar product over X , where X can beovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport either a space domain of R d or a space-time domain, will be denoted by ( · , · ) X ,the L -scalar product over subsets X of R d − will be denoted (cid:104)· , ·(cid:105) X . In both casesthe corresponding L -norm is denoted by (cid:107) · (cid:107) X . We will use the notation a (cid:46) b to denote a ≤ Cb with C a moderate constant independent of the mesh-parameterand the physical parameters of the problem, (except those that are assumed to beunity). We will also use the notation a ∼ b for a (cid:46) b and b (cid:46) a .In this paper the analysis will be restricted to velocity fields with a particularmultiscale character. We will assume that the problem is normalized so that U := (cid:107) β (cid:107) L ∞ ( Q ) = 1 and we assume that the characteristic lengthscale is given by L := 1,similarly we assume that T ∼ L/ (cid:107) β (cid:107) L ∞ ( Q ) = 1. Instead of making the standardassumption that (cid:107) β (cid:107) W , ∞ (Ω) is small, we assume that there is a decomposition ofthe velocity field, β = β + β (cid:48) , where, for all t , (cid:107) β (cid:107) W , ∞ (Ω) ∼ (cid:107) β (cid:48) (cid:107) L ∞ (Ω) ∼ µ . This allows us to define atimescale for the flow relating to the coarse scale spatial variation and the fine scaleamplitude, ( τ F ) − := sup t ∈ I ( (cid:107) β ( t ) (cid:107) W , ∞ (Ω) + (cid:107) β (cid:48) ( t ) (cid:107) L ∞ (Ω) /µ ) ∼ . (2.2)Essentially we assume that the velocity vectorfield can be decomposed in a coarsescale, responsible for transport, that is slowly varying in space and a fine scale,responsible for mixing, that has small amplitude but may have strong spatial vari-ation. Expressed in P´eclet numbers this means that the coarse scale P´eclet numbermay be arbitrarily high, whereas the fine scale P´eclet number, based on T , β (cid:48) and µ , must be O (1). A sharper value of τ F given a molecular diffusion µ and a velocityfield β may be obtained by solving a certain minimzation problem in the L ∞ -normthat will be detailed in the a priori analysis.We will assume that the coarse scale velocity satisfies a pointwise non-penetration condition on the domain boundary, β · n ∂ Ω = 0.The rough initial data or source term together with the high P´eclet numberand the multiscale character of the velocity field may lead to complex, low regular-ity solutions. More precisely solutions are smooth, due to parabolic regularity, butwith large Sobolev norms, rendering standard a priori error estimates based on ap-proximation theory worthless. Indeed classical global estimates for stabilized finiteelement methods for time dependent convection–diffusion equations , , yield thehigh P´eclet number error estimate:sup t ∈ I (cid:107) ( u − u h )( t ) (cid:107) Ω + (cid:107) µ ∇ ( u − u h ) (cid:107) Q ≤ Ch (1 + Pe − / h ) | u | L ( I ; H (Ω)) (2.3)where u h denotes the approximate solution using piecewise affine approximation.Even though | u | H (Ω) is huge in the presence of layers, stabilized methods are rele-vant for this case since one may derive localized error estimates, showing that per-turbations can not spread too far upwind of crosswind in the stationary case , , ,ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman or spread too far across characteristics in the transient case . The reason this worksis that | u | H ( ω ) can be assumed to be small in a part of the domain ω ⊂ Ω pro-vided u has no layers in a neighbourhood of ω . In the estimates the bad part canbe cut away using a suitably chosen weight function. This technique is not appli-cable in the present case, since the strong oscillations of the velocity field and thenonsmoothness of u and f , makes it unrealistic to assume that | u | H ( ω ) is smallin any part of the domain. The aim of this paper is to show that also in this case,stabilized finite element methods produce an improved solution compared to thatof standard Galerkin and that we can indeed derive error estimates for some largescale quantities defined by differential filtering.Drawing from earlier ideas on a posteriori error estimation , we propose toestimate a regularized error, or in other terms, work in a norm in between H − and L . Indeed let the regularized error ˜ e be defined by the partial differential equation − h ∆˜ e + ˜ e = ( u − u h )( · , T ) , (2.4)with ˜ e | ∂ Ω = 0, h ∈ R + . We then prove that (cid:107) ˜ e (cid:107) h := (cid:16) (cid:107) h ∇ ˜ e (cid:107) + (cid:107) ˜ e (cid:107) (cid:17) (cid:46) C T,f,u (cid:18) h h (cid:19) . The constant C T,f,u of the estimates is independent both of the P´eclet number andof the regularity of the exact solution. It depends on the data of the problem andthe main time dependence is a factorexp (cid:18) c Λ Tτ F (cid:19) where T denotes the final time of the simulation and c Λ is a moderate constant.This factor limits the applicability of the analysis to time intervals of the size τ F and hence limits the validity of the arguments to smooth vectorfields β or thosesatisfying the scale separation case discussed above.The parameter h can be related to a filter width δ by h := δ , or an artificialdiffusivity for elliptic smoothing. In both cases we see that the hidden constantincludes dimensional quantities of order unity. For h constant the results can beinterpreted as H − -norm error estimates and in this norm the convergence rate h is most likely optimal.Below we will first recall the standard L -norm error analysis for high P´ecletnumber flow and see what is required of data to obtain an estimate with an orderthat is fully independent of µ (that is, we control suitable Sobolev norms of u apriori). Then, we consider error estimation of filtered quantities and we obtain aposteriori error estimates. A priori estimates for rough solutions in the general casefollow directly using standard stability estimates of the discrete solution.The weak formulation of equation (2.1) reads, for t >
0, find u ∈ H (Ω) suchthat u ( x,
0) = u ( x ) and( ∂ t u, v ) Ω + a ( u, v ) = ( f, v ) Ω , ∀ v ∈ H (Ω) , (2.5)ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport where a ( · , · ) is defined by: a ( u, v ) := ( β · ∇ u, v ) Ω + ( µ ∇ u, ∇ v ) Ω .
3. Stability and regularity of the solution
Since the constant in the estimate (2.3) depends on the H -norm of the solution itis important to understand what quantities it is possible to control, if robustnessin µ is required. In this section we will first show the standard global regularityestimate for parabolic problems detailing the dependence of the constant on µ . Wewill then show an estimate where the constant is independent of µ , but not of τ F .In both cases we will need to assume that data have some smoothness. Lemma 3.1. (Standard energy estimate) Let u be the solution of (2.5) , then thereholds sup t ∈ I (cid:107) u ( t ) (cid:107) Ω + (cid:107) µ ∇ u (cid:107) Q (cid:46) (cid:90) I (cid:107) f ( t ) (cid:107) Ω d t + (cid:107) u (cid:107) Ω . (3.1) Proof.
We take v = u in (2.5) and then use that (cid:90) I ( f, u ) Ω d t ≤ sup t ∈ I (cid:107) u ( · , t ) (cid:107) Ω (cid:90) I (cid:107) f (cid:107) Ω d t from which the bound on sup t ∈ I (cid:107) u ( t ) (cid:107) Ω follows. This bound is then used to get the H -bound, which is a consequence of the equation (2.5) and the equality (cid:107) µ ∇ u (cid:107) Q = (cid:90) I a ( u, u ) d t. Theorem 3.1.
Let u be the weak solution of (2.1) , with f ∈ L ( Q ) and u ∈ H (Ω) . Assume that Ω is a convex domain. Then there holds sup t ∈ I (cid:107) µ ∇ u ( t ) (cid:107) Ω + | µu | L ( I ; H (Ω)) (cid:46) µ − / (cid:18)(cid:90) I (cid:107) f ( t ) (cid:107) Ω d t + (cid:107) u (cid:107) Ω (cid:19) + (cid:107) f (cid:107) Q + (cid:107) µ ∇ u (cid:107) Ω . Proof.
Multiplying (2.1) with − µ ∆ u , integrating over the space time domain Q ∗ := (0 , t ∗ ) × Ω, with t ∗ < T and applying the Cauchy-Schwarz inequality and thearithmetic-geometric inequality gives (cid:107) µ ∇ u ( t ∗ ) (cid:107) + (cid:107) µ ∆ u (cid:107) Q ∗ ≤ (cid:107) β (cid:107) L ∞ ( Q ∗ ) µ − (cid:107) µ ∇ u (cid:107) Q ∗ + 14 (cid:107) µ ∆ u (cid:107) Q ∗ + (cid:107) f (cid:107) Q ∗ + (cid:107) µ ∇ u (cid:107) . It follows from (3.1) thatsup t ∈ I (cid:107) µ ∇ u ( t ) (cid:107) Ω (cid:46) (cid:107) β (cid:107) L ∞ ( Q ) µ − / ( (cid:90) I (cid:107) f ( t ) (cid:107) Ω d t + (cid:107) u (cid:107) Ω ) + (cid:107) f (cid:107) Q + (cid:107) µ ∇ u (cid:107) Ω ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman and (cid:107) µ ∆ u (cid:107) Q (cid:46) (cid:107) β (cid:107) L ∞ ( Q ) µ − / ( (cid:90) I (cid:107) f ( t ) (cid:107) Ω d t + (cid:107) u (cid:107) Ω ) + (cid:107) f (cid:107) Q + (cid:107) µ ∇ u (cid:107) Ω . We conclude using elliptic regularity, recalling that the advection velocity is nor-malised.
Remark 3.1.
Observe that it follows that (cid:107)∇ u (cid:107) Q (cid:46) µ − / , sup t ∈ I (cid:107)∇ u (cid:107) Ω (cid:46) µ − and (cid:107) u (cid:107) L ( I ; H (Ω)) (cid:46) µ − / . Hence we conclude that the regularity estimates ofLemma 3.1 and Theorem 3.1 both are sensitive to the variation of the diffusivityand can only be used when Pe h ≤
1. Also some regularity of the inital data isneeded, u ∈ H (Ω). Revisiting the estimate (2.3) in this context we get (cid:107) µ ∇ ( u − u h ) (cid:107) Q (cid:46) (cid:104)(cid:18) hµ (cid:19) / + (cid:18) hµ (cid:19)(cid:105) ( (cid:90) I (cid:107) f ( t ) (cid:107) Ω d t + (cid:107) u (cid:107) Ω + µ (cid:107) f (cid:107) Q + (cid:107) µ ∇ u (cid:107) Ω ) . (3.2)This estimate is clearly useless when Pe h > O ( µ ) width layers in the solution u . The corresponding growth in the H -norm,makes global estimates independent of µ impossible. Theorem 3.2.
Let u be the weak solution of (2.1) , with Ω convex, f ∈ L ( I ; H (Ω)) and u ∈ H (Ω) . Then for < h , µ < , there holds sup t ∈ I (cid:107) u ( · , t ) (cid:107) h + | ( h µ ) u | L ( I ; H (Ω)) + T − (cid:107) h ∇ u (cid:107) Q + T − (cid:107) h ∂ t u (cid:107) Q (cid:46) C T ( h (cid:107) f (cid:107) L ( I ; H (Ω)) + (cid:90) I (cid:107) f (cid:107) d t + (cid:107) u (cid:107) h ) , (3.3) with C T = e (cid:16) c Λ TτF (cid:17) , where τ F is given by (2.2) and c Λ is a moderate constant . Proof.
First observe that we know from Theorem 3.1 that u ∈ L ( I ; H (Ω)), sothe important addition in Theorem 3.2 is the control of (cid:107) u ( · , t ) (cid:107) h that is robust µ .Multiply equation (2.1) by − h ∆ u and integrate over Q ∗ := (0 , t ∗ ) × Ω with t ∗ < T ,to obtain12 (cid:107) h ∇ u ( · , t ∗ ) (cid:107) − ( β · ∇ u, h ∆ u ) Q ∗ + (cid:107) ( µ h ) ∆ u (cid:107) Q ∗ (cid:46) ( ∇ f, h ∇ u ) Q ∗ + 12 (cid:107) h ∇ u (cid:107) . ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport The second term in the left hand side does not have a sign and requires some furtherconsideration. First split the velocity field in the large and the fine scale component, − ( β · ∇ u, h ∆ u ) Q ∗ = − ( β · ∇ u, h ∆ u ) Q ∗ − ( β (cid:48) · ∇ u, h ∆ u ) Q ∗ , then integrate by parts in the term representing the large scale transport, notingthat if t and t denotes the two orthonormal tangential vectors to ∂ Ω, β · ∇ u | ∂ Ω = β · n ∂ Ω (cid:124) (cid:123)(cid:122) (cid:125) =0 ∇ u · n ∂ Ω | ∂ Ω + (cid:88) i =1 β · t i ( ∇ u · t i (cid:124) (cid:123)(cid:122) (cid:125) =0 ) | ∂ Ω = 0 . Recalling that u ∈ L ( I ; H (Ω)) and β ∈ C ( I ; W , ∞ (Ω)) we have β · ∇ u ∈ L ( I ; H (Ω)). Then the following integration by parts is justified − ( β · ∇ u, h ∆ u ) Q ∗ = ( ∇ ( β · ∇ u ) , h ∇ u ) Q ∗ . Note that by the product rule( ∇ ( β · ∇ u ) , h ∇ u ) Q ∗ = d (cid:88) i =1 (( ∂ x i β ) · ∇ u, h ∂ x i u ) Q ∗ + d (cid:88) i =1 ( β · ( ∂ x i ∇ u ) , h ∂ x i u ) Q ∗ . (3.4)For the first sum of the right hand side we have d (cid:88) i =1 (( ∂ x i β ) · ∇ u, h ∂ x i u ) Q ∗ = (( ∇ S β ) · ∇ u, h ∇ u ) Q ∗ where ∇ S denotes the symmetric part of the gradient tensor. Similarly we obtainfor the second part d (cid:88) i =1 ( β · ( ∂ x i ∇ u ) , h ∂ x i u ) Q ∗ = d (cid:88) i =1 d (cid:88) j =1 ( β j ( ∂ x i ∂ x j u ) , h ∂ x i u ) Q ∗ = d (cid:88) i =1 d (cid:88) j =1 ( β j ( ∂ x j ∂ x i u ) , h ∂ x i u ) Q ∗ = d (cid:88) i =1 ( β · ∇ ∂ x i u, h ∂ x i u ) Q ∗ . (3.5)By the divergence theorem, recalling that β · n ∂ Ω = 0, we have d (cid:88) i =1 ( β · ∇ ∂ x i u, h ∂ x i u ) Q ∗ = − d (cid:88) i =1 ( ∇ · β ∂ x i u, h ∂ x i u ) Q ∗ . We conclude that, with I denoting the identity matrix, − ( β · ∇ u, h ∆ u ) Q ∗ = (( ∇ S β − ∇ · β I ) ∇ u, h ∇ u ) Q ∗ (3.6)Observing that( β (cid:48) · ∇ u, h ∆ u ) Q ∗ ≤ (cid:107) h | β (cid:48) | µ − / ∇ u (cid:107) Q ∗ (cid:107) ( h µ ) ∆ u (cid:107) Q ∗ ≤ (cid:107) h | β (cid:48) | µ − / ∇ u (cid:107) Q ∗ + 12 (cid:107) ( h µ ) ∆ u (cid:107) Q ∗ (3.7)ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman we have,12 (cid:107) h ∇ u ( · , t ∗ ) (cid:107) + 12 (cid:107) ( µ h ) ∆ u (cid:107) Q ∗ ≤ (cid:90) t ∗ ( ∇ f, h ∇ u ) Ω + 12 (cid:107) h ∇ u (cid:107) + (Λ( β , β (cid:48) , µ ) ∇ u, h ∇ u ) Ω (3.8)where Λ( β , β (cid:48) , µ ) := ∇ S β − ∇ · β I + 12 | β (cid:48) | µ − . Defining ˜ τ − F := sup t ∈ I inf β ∈ W , ∞ (Ω) β · n =0 on ∂ Ω σ p (Λ( β , β (cid:48) , µ )) (3.9)where σ p ( A ) denotes the largest positive eigenvalue of the matrix A , we may write (cid:107) h ∇ u ( · , t ∗ ) (cid:107) + (cid:107) ( µ h ) ∆ u (cid:107) Q ∗ ≤ t ∗ (cid:107) h ∇ f (cid:107) Q ∗ + (( t ∗ ) − + 2˜ τ − F ) (cid:107) h ∇ u (cid:107) Q ∗ + (cid:107) h ∇ u (cid:107) . (3.10)The claim regarding the control of the space derivatives now follows after an appli-cation of Gronwall’s lemma, yielding (cid:107) h ∇ u ( · , t ∗ ) (cid:107) + (cid:107) ( h µ ) ∆ u (cid:107) Q ∗ (cid:46) e t ∗ / ˜ τ F (cid:16) (cid:107) h ∇ f (cid:107) Q ∗ + (cid:107) h ∇ u (cid:107) (cid:17) . (3.11)Combining this estimate with estimate (3.1) yields the claim for the first two termsin the left hand side of equation (3.3). Note that we also have T − (cid:107) h ∇ u (cid:107) Q (cid:46) sup t ∈ I (cid:107) h ∇ u ( · , t ) (cid:107) (cid:46) e T/ ˜ τ F (cid:16) (cid:107) h ∇ f (cid:107) Q + (cid:107) h ∇ u (cid:107) (cid:17) . (3.12)For the control of the time derivative, simply multiply the equation with h ∂ t u andintegrate to obtain12 (cid:107) h ∂ t u (cid:107) Q ≤ − ( β · ∇ u, h ∂ t u ) Q + 12 (cid:107) h f (cid:107) Q + 12 (cid:107) ( h µ ) ∇ u (cid:107) ≤ (cid:107) β (cid:107) L ∞ ( Q ) (cid:107) h ∇ u (cid:107) Q (cid:107) h ∂ t u (cid:107) Q + 12 (cid:107) h f (cid:107) Q + 12 (cid:107) ( h µ ) ∇ u (cid:107) and we conclude using the estimate (3.12) and the observation that comparing thedefinitions (2.2) and (3.9) we have τ F ∼ ˜ τ F , with a moderate constant. Since Ω isconvex, elliptic regularity holds | ( µ h ) u | L ( I ; H (Ω)) (cid:46) (cid:107) ( µ h ) ∆ u (cid:107) Q . The claim follows.ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100
Robust error estimates in weak norms for advection dominated transport
4. Finite element discretizations
Let {T h } be a family of nonoverlapping conforming, quasi uniform triangulations, T h := { K } h where the simplices K have diameter h K and that is indexed by h := max h K . We let the set of interior faces { F } h of a triangulation T h be denotedby F .We will consider a standard finite element space of piecewise polynomial, con-tinuous functions V h := { v h ∈ H (Ω) : v h | K ∈ P k ( K ) , ∀ K ∈ T h } , where P k ( K ) denotes the polynomials of degree less than or equal to k on K . Thefollowing inverse inequalities are known to hold on V h , (cid:107)∇ v h (cid:107) K ≤ c i h − K (cid:107) v h (cid:107) K (4.1)and (cid:107) v h (cid:107) ∂K ≤ c t h − / K (cid:107) v h (cid:107) K . (4.2)We let π h : L (Ω) → V h denote the L -projection defined by π h v ∈ V h such that( π h v, w h ) Ω = ( v, w h ) Ω for all w h ∈ V h .The standard finite element method is then obtained by restricting the weakformulation (2.5) to the discrete space V h . For t > u h ∈ V h such that u h ( x,
0) = π h u ( x ) and( ∂ t u h , v h ) Ω + a ( u h , v h ) = ( f, v h ) Ω , ∀ v h ∈ V h . (4.3)Taking v h = u h we immediately get the stability estimatesup t ∈ (0 ,T ] (cid:107) u h ( t ) (cid:107) Ω + (cid:32)(cid:90) T (cid:107) µ ∇ u h (cid:107) d t (cid:33) (cid:46) (cid:90) T (cid:107) f (cid:107) Ω d t + (cid:107) u (cid:107) Ω . When the P´eclet number is low this approximations of the parabolic equationmay be analysed using well known finite element techniques, see Thom´ee .However when the local P´eclet number is high, the problem is a singularlyperturbed parabolic problem and the stability properties of the standard Galerkinmethod are in general insufficient for optimal convergence. In particular in thepresence of layers the whole computational domain may be polluted by spuriousoscillations, but as can be seen in Figure 1, convergence order is lost also for smoothsolutions. In this case we study a Gaussian function convected one turn in a disc.Time discretization is performed with the Crank-Nicolson method and we comparethe result of the standard Galerkin method with those obtained using symmetricstabilization methods.In this paper we will consider symmetric stabilization methods only, however atleast for time constant β one can obtain similar results for the SUPG-method.ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman
Fig. 1. Convergence in the L -norm plotted against the timestep τ = h , of the stabilized (dashedline for explicit and dotted line for implicit treatment of the stabilization operator) and unstabilized(dash-dot) Crank-Nicolson method. The full line shows optimal second order convergence. Symmetric stabilization methods
In the last ten years there has been an important development in the field of highorder symmetric stabilization methods. Such methods are typically obtained by theaddition of a weakly consistent, dissipative operator to the formulation. Some ofthe more important symmetric stabilization methods present in the literature arethe subgrid viscosity method suggested by Guermond , the orthogonal subscalemethod proposed by Codina , the local project method introduced by Becker andBraack , the discontinuous Galerkin method and the continuous interior penalty(CIP) method suggested by Douglas and Dupont and analysed by Burman andHansbo . The analysis that we propose herein will rely on an orthogonality argu-ment and hence is valid only for the orthogonal subscales, the discontinuous Galerkinmethod and for the CIP-method. For simplicity we will focus on the latter below.Herein we will assume that the following holds:(1) s h ( · , · ) : V h × V h (cid:55)→ R is a symmetric, bilinear and positive semi definiteoperator.(2) For the L -projection π h : L (Ω) (cid:55)→ V h there holds– Approximability, ∀ u ∈ H (Ω) ∩ H k +1 (Ω) (cid:107) h − / ( u − π h u ) (cid:107) Ω + (cid:107) µ ∇ ( u − π h u ) (cid:107) Ω (cid:46) ( h + µ ) h k | u | H k +1 (Ω) , (4.4)– Stability and weak consistency of the stabilization operator s h ( u h , π h v ) (cid:46) s h ( u h , u h ) h (cid:107)∇ v (cid:107) Ω , ∀ u h ∈ V h , v ∈ H (Ω) . (4.5) s h ( π h u, π h u ) (cid:46) h k +1 / | u | H k +1 (Ω) , ∀ u ∈ H k +1 (Ω) . (4.6)ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport – Enhanced continuity of the convection term: for all v h ∈ V h and β ∈ W , ∞ (Ω), with β · n ∂ Ω | ∂ Ω = 0, there holds | ( u − π h u, β · ∇ v h ) Ω | (cid:46) (cid:107) β (cid:107) W , ∞ (Ω) (cid:107) u − π h u (cid:107) Ω (cid:107) v h (cid:107) Ω + (cid:107) β (cid:107) L ∞ (Ω) (cid:107) h − / ( u − π h u ) (cid:107) Ω s h ( v h , v h ) . (4.7)As an example consider the CIP-method. Here s h ( · , · ) consists in a penalty on thejump of the gradient over element faces, and takes the form s h ( u h , v h ) := γ (cid:88) F ∈F (cid:10) h F (cid:107) β · n F (cid:107) L ∞ ( F ) [[ ∇ u h · n F ]] , [[ ∇ v h · n F ]] (cid:11) F where F denotes the faces in the mesh, [[ x ]] the jump of x over F , the orientationis not important, n F a fixed but arbitrary normal associated to each face. In theanalysis below we will for simplicity use this stabilization operator.The stabilized finite element method then takes the form, for t > u h ∈ V h such that u h ( x,
0) = π h u ( x ) and( ∂ t u h , v h ) Ω + a ( u h , v h ) + s h ( u h , v h ) = ( f, v h ) Ω , ∀ v h ∈ V h . (4.8)For the satisfaction of the assumptions (4.4)–(4.7) we refer to Ref. 4. Althoughin that reference Nitsche-type boundary conditions are used, the same argumentmay be shown to work whenever β · ∇ v h | ∂ Ω = 0, which is the case herein. Forcompleteness we show how to obtain (4.7) in our case. We also show that theflow field β in the stabilization operator may be replaced by β at the cost of anonessential perturbation. The key result is the following Lemma. Lemma 4.1.
Assume that for k ≥ no element in T h has more than one faceintersecting ∂ Ω . Then there holds for all u h ∈ V h , inf v h ∈ V h (cid:107) h ( π β · ∇ u h − v h ) (cid:107) Ω (cid:46) s h ( u h , u h ) + h (cid:107) β (cid:107) W , ∞ (Ω) (cid:107) u h (cid:107) Ω , where π β is a piecewise constant approximation of β that will be defined below. Proof.
Let π β be the projection onto element wise constants such that for every K that has no face on the boundary there holds (cid:90) K π β d x = (cid:90) K β d x. On elements adjacent to the boundary define π β by (cid:90) ∂K ∩ ∂ Ω π β · n ∂ Ω d x = (cid:90) ∂K ∩ ∂ Ω β · n ∂ Ω d x and (cid:90) ∂K ∩ ∂ Ω π β · t i d x = (cid:90) ∂K ∩ ∂ Ω β · t i d x, i = 1 , ..., d − . It follows from standard approximation that for all K ∈ T h , (cid:107) β − π β (cid:107) L ∞ ( K ) (cid:46) (cid:107) β (cid:107) W , ∞ ( K ) h K . ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman
Note that for any element K with one face on the boundary there holds π β · ∇ u h | ∂ Ω ∩ ∂K = π β · n ∂ Ω (cid:124) (cid:123)(cid:122) (cid:125) =0 ∇ u h · n ∂ Ω | ∂ Ω ∩ ∂K + d (cid:88) i =1 π β · t i ( ∇ u h · t i (cid:124) (cid:123)(cid:122) (cid:125) =0 ) | ∂ Ω ∩ ∂K = 0 . For k = 1, if an element K has two faces on the boundary, then u h ≡ π β · ∇ u h | ∂ Ω ∩ ∂K = 0. It then follows using the same argumentsas in Ref. 4 thatinf v h ∈ V h (cid:107) h ( π β · ∇ u h − v h ) (cid:107) Ω ≤ (cid:32) (cid:88) F ∈F (cid:107) h [[ π β · ∇ u h ]] (cid:107) F (cid:33) . By adding and subtracting β inside the jump and by applying a triangle inequalityfollowed by a trace inequality we have (cid:32) (cid:88) F ∈F (cid:107) h [[ π β · ∇ u h ]] (cid:107) F (cid:33) (cid:46) s h ( u h , u h ) + h (cid:107) β (cid:107) W , ∞ (Ω) (cid:107) u h (cid:107) Ω and the proof is finished.The continuity (4.7) is now obtained by adding and subtracting π β in the rightslot of the left hand side of the equation | ( u − π h u, β ·∇ v h ) Ω | ≤ | ( u − π h u, ( β − π β ) ·∇ v h ) Ω | + inf w h ∈ V h | ( u − π h u, π β ·∇ v h − w h ) | Ω (cid:46) (cid:107) β (cid:107) W , ∞ (Ω) (cid:107) u − π h u (cid:107) Ω (cid:107) v h (cid:107) Ω + (cid:107) β (cid:107) L ∞ (Ω) (cid:107) h − / ( u − π h u ) (cid:107) Ω s h ( v h , v h ) . (4.9)where a Cauchy-Schwarz inequality, an inverse inequality, the approximation prop-erties of π β and Lemma 4.1 have been used in the last inequality. Lemma 4.2.
Let ¯ s h ( u h , v h ) := γ (cid:88) F ∈F (cid:10) h F (cid:107) β · n F (cid:107) L ∞ ( F ) [[ ∇ u h · n F ]] , [[ ∇ v h · n F ]] (cid:11) F . Then s h ( u h , u h ) (cid:46) ¯ s h ( u h , u h ) + h (cid:107) µ ∇ u h (cid:107) and ¯ s h ( u h , u h ) (cid:46) s h ( u h , u h ) + h (cid:107) µ ∇ u h (cid:107) . Proof.
The proof for the upper and the lower bounds are similar, so we consideronly the first inequality. By using the decomposition of β and a triangular inequalitywe have s h ( u h , u h ) (cid:46) γ (cid:88) F ∈F (cid:107) h F (cid:107) β · n F (cid:107) L ∞ ( F ) [[ ∇ u h ]] (cid:107) F + γ (cid:88) F ∈F (cid:107) h F (cid:107) β (cid:48) · n F (cid:107) L ∞ (Ω) [[ ∇ u h ]] (cid:107) F . ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport Using now the size constraint on β (cid:48) and a trace inequality we conclude γ (cid:88) F ∈F (cid:107) h F (cid:107) β (cid:48) · n F (cid:107) L ∞ (Ω) [[ ∇ u h ]] (cid:107) F (cid:46) γ (cid:88) K ∈T h (cid:107) h K µ ∇ u h (cid:107) K (cid:46) γh (cid:107) µ ∇ u h (cid:107) . Stability and convergence results
We will now recall the main results on stability and convergence of symmetricstabilization methods. These results are minor modifications of those in Ref. 3. Forconvenience we introduce a triple norm associated to the stabilized method. Let |(cid:107) u h (cid:107)| h := (cid:90) I (cid:16) (cid:107) µ ∇ u h (cid:107) + s h ( u h , u h ) (cid:17) d t. Following the proof of Lemma 3.1 it is straightforward to derive the following sta-bility and error estimates.
Lemma 4.3. (Stability)Let u h be the solution of (4.8) , with γ ≥ . Then sup t ∈ I (cid:107) u h ( t ) (cid:107) Ω + |(cid:107) u h (cid:107)| h (cid:46) (cid:90) T (cid:107) f (cid:107) Ω d t + (cid:107) u (cid:107) Ω . Applying stability to the error equation leads to the following error estimate:
Theorem 4.1. (Convergence CIP-method)Let u ∈ L ( I ; H r (Ω)) , r ≥ , be the solution of (2.5) and u h the solution of (4.8) with γ > . Then there holds sup t ∈ I (cid:107) ( u − u h )( t ) (cid:107) Ω + |(cid:107) u − u h (cid:107)| h (cid:46) ( τ − F T h + ( τ − F h + γ + γ − ) h + µ ) h s − | u | L ( I ; H s (Ω)) , with s = min( r, k + 1) and where we used that (cid:107) β (cid:107) L ∞ ( Q ) ∼ . Proof.
For simplicity we only give the proof in the form of a final time result. Let e h := u h − π h u , using the same arguments as for the stability Lemma 3.1 we have12 (cid:107) e h ( T ) (cid:107) + |(cid:107) e h (cid:107)| h = (cid:90) I [( ∂ t e h , e h ) Ω + a ( e h , e h ) + s h ( e h , e h )] d t. Using Galerkin orthogonality and the orthogonality of the L -projection we have12 (cid:107) e h ( T ) (cid:107) + |(cid:107) e h (cid:107)| h = (cid:90) I [ a ( u − π h u, e h ) − s h ( π h u, e h )] d t. Using the decomposition of the velocity we have a ( u − π h u, e h ) (cid:46) ( u − π h u, β ·∇ e h ) Ω +( u − π h u, β (cid:48) ·∇ e h ) Ω + (cid:107) µ ∇ ( u − π h u ) (cid:107) Ω (cid:107) µ ∇ e h (cid:107) Ω ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman and using the continuity (4.7) in the first term of the right hand side and the boundon β (cid:48) in the second we have a ( u − π h u, e h ) (cid:46) (cid:107)∇ β (cid:107) L ∞ (Ω) (cid:107) ( u − π h u ) (cid:107) Ω (cid:107) e h (cid:107) Ω + ( (cid:107) β (cid:107) L ∞ (Ω) (cid:107) h − / ( u − π h u ) (cid:107) Ω + τ − F (cid:107) u − π h u (cid:107) Ω + (cid:107) µ ∇ ( u − π h u ) (cid:107) Ω ) |(cid:107) e h (cid:107)| h . It follows after a Cauchy-Schwarz inequality and an arithmetic-geometric inequalityin the right hand side that (cid:107) e h ( T ) (cid:107) + |(cid:107) e h (cid:107)| h (cid:46) s h ( π h u, π h u ) + (cid:107)∇ β (cid:107) L ∞ ( Q ) T (cid:107) h − / ( u − π h u ) (cid:107) Q + ( τ − F + (cid:107) β (cid:107) L ∞ ( Q ) h − ) (cid:107) u − π h u (cid:107) Q + (cid:107) µ ∇ ( u − π h u ) (cid:107) Q + T − (cid:107) e h (cid:107) Q . We conclude by applying Gronwall’s lemma and the approximation results of (4.4)and (4.6).For high mesh P´eclet numbers this estimate is sub optimal optimal with O ( h )in the L ∞ ( I ; L (Ω))-norm and for low mesh P´eclet numbers it is optimal in the L ( I ; H (Ω))-norm. In the latter case the convergence in the L -norm can be im-proved under certain assumptions on the time variation of β . Note that thisestimate does not have exponential growth, however in the high P´eclet case thatfactor is hidden in the Sobolev norm of the exact solution. Combining this conver-gence result with the regularity result of Theorem 3.2 we may prove the followingestimate that is fully independent of µ , in the sense that we also control the Sobolevnorm in the constant. Note however that smoothness of the source term and theinitial data is required. Corollary 4.1.
Let f ∈ L ( I ; H (Ω)) and u ∈ H (Ω) , let u denote the solutionof (2.5) and u h the solution of (4.8) , with γ > and assume that Pe h >> . Then (cid:107) ( u − u h )( T ) (cid:107) Ω (cid:46) C T ( T (1 + τ − F ) + 1) h ( (cid:107) f (cid:107) L ( I ; H (Ω)) + (cid:107)∇ u (cid:107) Ω ) . Proof.
An immediate consequence of the Theorem 4.1 and Theorem 3.2.
5. Perturbation equation and the dual problem
The error analysis in weak norms uses a perturbation equation and an associateddual problem. Taking the difference of the two formulations (2.5) and (4.3), setting e = u − u h and integrating by parts we obtain( ∂ t e, ϕ ) Ω + a ( e, ϕ ) = ( e ( T ) , ϕ ( T )) Ω − ( e (0) , ϕ (0)) Ω − ( e, ∂ t ϕ + β ·∇ ϕ ) Ω +( µ ∇ e, ∇ ϕ ) Ω . This suggests the adjoint equation, find ϕ ∈ H (Ω) such that − ∂ t ϕ − β · ∇ ϕ − µ ∆ ϕ = 0 in Qϕ = 0 on ∂ Ω × Iϕ ( · , T ) = Ψ( · ) in Ω , (5.1)ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport where Ψ ∈ H (Ω). Then the following error representation holds( e ( T ) , Ψ) Ω = ( e (0) , ϕ (0)) Ω + ( ∂ t e, ϕ ) Ω + a ( e, ϕ ) . (5.2)We will now proceed and discuss the choice of Ψ and the associated stability esti-mates on ϕ . Regularization of the error and weak norms
Since it appears not to be possible to prove a posteriori error estimates in the L -norm independently of the P´eclet number, unless one ressorts to a saturationassumption, we will here consider a regularized error, where a parameter h (thatmay ultimately depend on h ) sets the scale of the regularization. We recall theproblem (2.4) for a given computational error e := u − u h , find ˜ e such that ˜ e | ∂ Ω = 0and − h ∆˜ e + ˜ e = e ( · , T ) . On weak form the problem writes: find ˜ e ∈ H (Ω) such that( h ∇ ˜ e, ∇ v ) Ω + (˜ e, v ) Ω = ( e, v ) Ω , ∀ v ∈ H (Ω) . (5.3)This is what is commonly called a Helmholtz filter or a differential filter, althoughit is not properly speaking a filter. The key observation here is that when h is small,the filtered solution is close to the solution where ever the solution is smooth. Closeto layers or other strongly localised features of e , ˜ u − u may be large locally, alsofor h small. Associated with the problem (5.3) we have the norm (cid:107) ˜ e (cid:107) h := (cid:107) h ∇ ˜ e (cid:107) + (cid:107) ˜ e (cid:107) and an associated relation, obtained by testing (5.3) with v = ˜ e , (cid:107) ˜ e (cid:107) h = ( e, ˜ e ) Ω . (5.4)We deduce from (5.4) that choosing Ψ = ˜ e in (5.1) above leads to the followingerror representation for (cid:107) ˜ e (cid:107) h : (cid:107) ˜ e (cid:107) h = ( e (0) , ϕ (0)) Ω + ( ∂ t e, ϕ ) Ω + a ( e, ϕ ) . (5.5) Stability of the dual solution
The advantage of using the dual technique is that instead of relying on regularityestimates for the exact solutions we can use regularity of the adjoint equation, whichmay be better behaved provided the data is well chosen. The following Theoremgives a precise characterization of the regularity of the dual problem in the multiscaleframework.
Theorem 5.1.
Let ϕ ( x, t ) be the weak solution of (5.2) , with Ψ = ˜ e , where ˜ e isdefined by (5.3) , with u and u h solutions of (2.5) and (4.8) respectively. Then thereholds sup t ∈ I (cid:107) ϕ ( · , t ) (cid:107) h + T − (cid:107) h ∇ ϕ (cid:107) Q + T − (cid:107) h ∂ t ϕ (cid:107) Q + | ( h µ ) ϕ | L ( I ; H (Ω)) (cid:46) C T (cid:107) ˜ e (cid:107) h , ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman with C T = e (cid:16) c Λ TτF (cid:17) where τ F is given by (2.2) and c Λ is a moderate constant. Proof.
The dual problem is equivalent to the forward problem after a change ofvariable ˜ t = T − t and ˜ x = − x , with the source term f = 0 and the initial data u = ˜ e ∈ H (Ω). The result then follows from Theorem 3.2.
6. Error estimates
In this section we will prove estimates for (cid:107) ˜ e (cid:107) h where the constant is robust in µ (under our assumptions on the data).We will only consider the case of semidiscretization in space and show how to prove an a posteriori error estimate, wherethe stability constant is essentially C T of Theorem 5.1. The a priori error estimatethen follows using the fact that the a posteriori residuals are a priori controlledby the discrete stability estimate of Lemma 4.3. Then we will consider the case ofinsufficient data, i.e. when only β is known, and show that under our assumptionson data we can obtain an upper bound of the error also in this case, where anonconsistent part limits the asymptotic convergence. In the regime that we areinterested in however this part is smaller than the discretization error.In the low mesh P´eclet number regime, we need to modify our estimate to obtainoptimality. In particular we need to use elliptic regularity to obtain optimality. Wehave not written the two estimates in a unified manner, since the natural form of theresidual quantities uses different norms. We outline the differences for estimationin the low mesh P´eclet number regime in a remark below. A posteriori and a priori error estimates
To prove a posteriori error estimates we use a duality technique together with the apriori control of the dual solution. In practice one may need to ressort to numericalsolution of the dual problem.
Theorem 6.1. (A posteriori error estimate) Let ˜ e be defined by (5.3) , with u and u h solutions of (2.5) and (4.8) respectively. Then there holds (cid:107) ˜ e (cid:107) h (cid:46) C T (cid:18) h h (cid:19) (cid:16)(cid:90) I inf v h ∈ V h (cid:107) h ( β · ∇ u h − v h ) (cid:107) Ω d t + (cid:90) I (cid:32) inf v h ∈ V h (cid:88) K ∈T h (cid:107) h ( f + µ ∆ u h − v h ) (cid:107) K (cid:33) d t + (cid:90) I (cid:32) (cid:88) F ∈F (cid:107) µ [[ ∇ u h · n F ]] (cid:107) F (cid:33) d t + (cid:90) I s h ( u h , u h ) d t + h (cid:107) u − π h u (cid:107) Ω (cid:17) . (6.1) The constant C T is defined in Theorem 3.2. Proof.
Starting from (5.5) and using Galerkin orthogonality and the orthogonalityovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100
Robust error estimates in weak norms for advection dominated transport of the L -projection we have (cid:107) ˜ e (cid:107) h = ( e (0) , ϕ (0) − π h ϕ ) Ω + ( ∂ t e, ϕ − π h ϕ ) Q + (cid:90) I ( a ( e, ϕ − π h ϕ ) + s h ( u h , π h ϕ )) d t. (6.2)Using the weak formulation (2.5) and the orthogonality of the L -projection wemay write (cid:107) ˜ e (cid:107) h = ( e (0) , ϕ (0) − π h ϕ (0)) Ω + ( β · ∇ u h − v h , π h ϕ − ϕ ) Q + ( µ ∇ u h , ∇ ( π h ϕ − ϕ )) Q + ( f, ϕ − π h ϕ ) Q + (cid:90) I s h ( u h , π h ϕ ) d t. (6.3)Considering the right hand side term by term we get using Cauchy-Schwarz inequal-ity, approximation and the stability of Theorem 5.1( e (0) , ϕ (0) − π h ϕ (0)) Ω ≤ h h − / (cid:107) e (0) (cid:107) Ω sup t ∈ I (cid:107) h ∇ ϕ (cid:107) Ω (cid:46) C T (cid:18) h h (cid:19) h (cid:107) e (0) (cid:107) Ω (cid:107) ˜ e (cid:107) h , (6.4)( β · ∇ u h − v h , π h ϕ − ϕ ) Q ≤ h h − / (cid:90) I (cid:107) h ( β · ∇ u h − v h ) (cid:107) Ω d t sup t ∈ I (cid:107) h ∇ ϕ ( · , t ) (cid:107) Ω (cid:46) C T (cid:18) h h (cid:19) (cid:90) I (cid:107) h ( β · ∇ u h − v h ) (cid:107) Ω d t (cid:107) ˜ e (cid:107) h . (6.5)In the third term we first integrate by parts on each element and then proceed withtrace inequalities, followed by approximation for the dual solution,( f, ϕ − π h ϕ ) Q − ( µ ∇ u h , ∇ ( ϕ − π h ϕ )) Q (cid:46) (cid:90) I (cid:32) inf v h ∈ V h (cid:88) K ∈T h (cid:107) h ( f + µ ∆ u h − v h ) (cid:107) K (cid:33) + (cid:32) (cid:88) F ∈F (cid:107) µ [[ ∇ u h · n F ]] (cid:107) F (cid:33) d t (cid:124) (cid:123)(cid:122) (cid:125) R ( u h ,µ ) × (cid:18) sup t ∈ I (cid:107) h − / ( π h ϕ − ϕ )( · , t ) (cid:107) Ω + sup t ∈ I (cid:107) ( π h ϕ − ϕ )( · , t ) (cid:107) F (cid:19) (cid:46) (cid:18) h h (cid:19) R ( u h , µ ) sup t ∈ I (cid:107) h ∇ ϕ ( · , t ) (cid:107) Ω (cid:46) C T (cid:18) h h (cid:19) R ( u h , µ ) (cid:107) ˜ e (cid:107) h . (6.6)Finally the stabilization term is handled using a Cauchy-Schwarz inequality, fol-ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman lowed by a trace inequality and the H -stability of the L -projection. (cid:90) I s h ( u h , π h ϕ ) d t (cid:46) (cid:107) β (cid:107) L ∞ ( Q ) (cid:90) I s ( u h , u h ) d t h sup t ∈ I (cid:107)∇ ϕ ( · , t ) (cid:107) Ω (cid:46) C T (cid:18) h h (cid:19) (cid:90) I s ( u h , u h ) d t (cid:107) ˜ e (cid:107) h . (6.7)The claim follows by collecting the upper bounds (6.4)-(6.7) and dividing by (cid:107) ˜ e (cid:107) h . Theorem 6.2. (A priori error estimate) Let ˜ e be defined by (5.3) , with u thesolution of (2.5) and and u h the solution (4.8) , with γ > . Assume that Pe h > then there holds (cid:107) ˜ e (cid:107) h (cid:46) C T (cid:18) h h (cid:19) (1 + h + T ) (cid:32)(cid:90) T (cid:107) f (cid:107) Ω d t + (cid:107) u (cid:107) Ω (cid:33) . The constant C T is defined in Theorem 3.2. Proof.
The result follows from estimate (6.1) by bounding all the residual termsusing Lemma 4.3. First observe that since the P´eclet number is high we may useinverse inequalities and trace inequalities to show that R ( u h , µ ) (cid:46) h (cid:90) T (cid:107) f (cid:107) Ω d t + T (cid:107) µ ∇ u h (cid:107) Q . For the contributions on the faces we have used (cid:90) I (cid:107) µ [[ ∇ u h · n F ]] (cid:107) F d t (cid:46) (cid:107) β (cid:107) L ∞ (Ω) h T (cid:32)(cid:90) I (cid:88) F ∈F (cid:107) µ [[ ∇ u h · n F ]] (cid:107) F d t (cid:33) (cid:46) (cid:107) µ ∇ u h (cid:107) Q (cid:46) |(cid:107) u h (cid:107)| h . (6.8)Similarly using a Cauchy-Schwarz inequality in time and the stability (4.3) we have (cid:90) I s ( u h , u h ) d t ≤ T (cid:90) I s ( u h , u h ) d t. We finally consider the first term on the right hand side of (6.1). First note that (cid:90) I inf v h ∈ V h (cid:107) h ( β ·∇ u h − v h ) (cid:107) Ω d t ≤ T inf v h ∈ V h (cid:107) h ( β ·∇ u h − v h ) (cid:107) Q + T (cid:107) h β (cid:48) ·∇ u h (cid:107) Q . By Lemma 4.1 and Lemma 4.2 we haveinf v h ∈ V h (cid:107) h ( β · ∇ u h − v h ) (cid:107) Q (cid:46) h (cid:107) β (cid:107) W , ∞ (Ω) (cid:107) u h (cid:107) Q + (cid:18)(cid:90) I ¯ s h ( u h , u h ) d t (cid:19) (cid:46) h (cid:107) β (cid:107) W , ∞ (Ω) (cid:107) u h (cid:107) Q + h (cid:107) µ ∇ u h (cid:107) Q + (cid:18)(cid:90) I s h ( u h , u h ) d t (cid:19) . ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport It follows thatinf v h ∈ V h (cid:107) h ( β · ∇ u h − v h ) (cid:107) Q (cid:46) max( h T (cid:107) β (cid:107) W , ∞ (Ω) , (cid:107) β (cid:107) L ∞ (Ω) )(sup t ∈ I (cid:107) u h ( · , t ) (cid:107) Ω + |(cid:107) u h (cid:107)| h ) . (6.9)Using the assumption on the small scale fluctuations (cid:107) β (cid:48) (cid:107) L ∞ ( Q ) (cid:46) µ we have (cid:107) h β (cid:48) · ∇ u h (cid:107) Q (cid:46) h (cid:107) µ ∇ u h (cid:107) Q ≤ h |(cid:107) u h (cid:107)| h . We conclude by collecting terms and applying Lemma 4.3.
Remark 6.1. (The necessity of stabilization for robustness) Note that the stabilityof the dual problem holds regardless of the numerical method used. The stabilizationin the numerical method allows us to control the first residual in the a posteriorierror estimate, by using the discrete stability estimate of Lemma 4.3. If no stabiliza-tion is present there is no control of the streamline derivative making it impossibleto obtain uniformity in µ . Another observation that is worthwhile is that the abovea priori error estimate is valid only for high mesh P´eclet number. This is becauseTheorem 6.1 is optimal only in this regime. For low mesh P´eclet number we mayinstead use the stability | ( h µ ) ϕ | L ( I ; H (Ω)) ≤ (cid:107) ˜ e (cid:107) h in the various bounds above.We only detail how equation (6.5) and (6.6) are modified( f − β · ∇ u h , ϕ − π h ϕ ) Q − ( µ ∇ u h , ∇ ( ϕ − π h ϕ )) Q (cid:46) (cid:32)(cid:90) I inf v h ∈ V h (cid:88) K ∈T h (cid:107) h ( f − β · ∇ u h + µ ∆ u h − v h ) (cid:107) K + (cid:88) F ∈F (cid:107) µh [[ ∇ u h · n F ]] (cid:107) F d t (cid:33) (cid:124) (cid:123)(cid:122) (cid:125) R ( u h ,µ ) × (cid:18)(cid:90) I (cid:107) h − ( π h ϕ − ϕ )( · , t ) (cid:107) + (cid:107) h − ( π h ϕ − ϕ ) (cid:107) F d t (cid:19) (cid:46) (cid:18) h h µ (cid:19) R ( u h , µ ) | ( µ h ) ϕ | L ( I ; H (Ω)) (cid:46) C T (cid:18) h h µ (cid:19) R ( u h , µ ) (cid:107) ˜ e (cid:107) h . (6.10)Finally using this modified version of Theorem 6.1, an a priori result with conclusionsimilar to that of Theorem 6.2 holds for u h solution of (4.8) with γ ≥
0. Thestabilization may be omitted when Pe h < The degenerate case of unknown β (cid:48) In many relevant cases β (cid:48) may be unknown or only partially known. If the statisticsof β (cid:48) are known some stochastic method may be used to recover expectancy valuesfor the solution. In this section we will consider the situation, that β (cid:48) is simplyovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman excluded from the computation and we will show that under our assumptions onthe small scale velocity fluctuations the error estimates still hold for high meshP´eclet numbers. Indeed in the high P´eclet number regime the consistency errormade by dropping the fine scale fluctuations of β is smaller than the discretizationerror.Here we use an advective field β that we assume is divergence free and let β bereplaced by β in (4.8). Theorem 6.3.
Let u be the solution of (2.5) and u h the solution of (4.8) , with β instead of β for the advective field and γ > , assume that Pe h >> , then (cid:107) ˜ e (cid:107) h (cid:46) C T (cid:18) h h (cid:19) (1 + h + T + T ) (cid:32)(cid:90) T (cid:107) f (cid:107) Ω d t + (cid:107) u (cid:107) Ω (cid:33) . Proof.
We proceed as in the proofs of Theorems 6.1 and 6.2 (cid:107) ˜ e (cid:107) h = ( e (0) , ( ϕ − π h ϕ )(0)) Ω + ( ∂ t e, ϕ − π h ϕ ) Q + (cid:90) I (cid:16) a ( e, ϕ − π h ϕ ) + ¯ s h ( u h , π h ϕ ) + ( β · ∇ u h − β · ∇ u h (cid:124) (cid:123)(cid:122) (cid:125) = − β (cid:48) ·∇ u h , π h ϕ ) Ω dt (cid:17) . (6.11)The only thing that differs from the previous analysis is the last term in the righthand side. Except for this term the proof proceeds as before. We will therefore hereonly show how this term may be bounded. After an integration by parts we have( β (cid:48) · ∇ u h , π h ϕ ) Q = ( u h , β (cid:48) · ∇ π h ϕ ) Q ≤ sup t ∈ I (cid:107) u h ( · , t ) (cid:107) Ω T (cid:107) β (cid:48) (cid:107) L ∞ ( Q ) h − / sup t ∈ I (cid:107) h ∇ π h ϕ ( · , t ) (cid:107) Ω (cid:46) sup t ∈ I (cid:107) u h ( · , t ) (cid:107) Ω T (cid:107) β (cid:48) (cid:107) L ∞ ( Q ) h − / sup t ∈ I (cid:107) h ∇ ϕ ( · , t ) (cid:107) Ω . Note that under our assumptions (cid:107) β (cid:48) (cid:107) L ∞ ( Q ) ∼ µ ≤ (cid:107) β (cid:107) L ∞ ( Q ) h leading to( β (cid:48) · ∇ u h , π h ϕ ) Q (cid:46) (cid:18) h h (cid:19) T C T (cid:107) ˜ e (cid:107) h . It follows that the consistency error is of the same order as the discretizationerror. If the P´eclet number is large, the contribution from the discretization errorcan be assumed to be dominating and the same order of convergence as in theunperturbed case should be observed, until the P´eclet number becomes so smallthat the inconsistency dominates. This means that in the high P´eclet regime, ifdata are known to be rough, noise in the velocities satisying the constraint on β (cid:48) may be neglected.ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 Robust error estimates in weak norms for advection dominated transport Conclusion
We have derived robust a posteriori and a priori error estimates for transientconvection–diffusion equations. The upshot is that the estimates are completelyrobust with respect to the P´eclet number, in the sense that we also control theSobolev norms of the exact solution in the error constant. The estimates allow forlow regularity data and multiscale advection, that may have strong spatial variationon the fine scale under a special scale separation condition. The aim of this workwas to take a first step towards an understanding of what transport problems arecomputable in the high P´eclet regime, beyond the standard assumption of smoothdata. These results also set a baseline for what should be achieved theoreticallyin the analysis of more involved methods, such as multiscale methods, in order toclaim that they produce an accuracy beyond what is obtained using a standardstabilized finite element method.
Acknowledgment
The author wishes to thank Professor Vivette Girault and Professor Alexandre Ernfor helpful advice.
References
1. R. Becker and M. Braack. A two-level stabilization scheme for the Navier-Stokes equa-tions. In
Numerical mathematics and advanced applications pages 123–130 (Springer,Berlin, 2004).2. E. Burman.
Adaptive finite element methods for compressible two-phase flows . PhDthesis, Chalmers University of Technology, 1998.3. E. Burman and M. A. Fern´andez. Finite element methods with symmetric stabilizationfor the transient convection-diffusion-reaction equation.
Comput. Methods Appl. Mech.Engrg. (2009) 2508–2519.4. E. Burman, M. A. Fern´andez, and P. Hansbo. Continuous interior penalty finite ele-ment method for Oseen’s equations.
SIAM J. Numer. Anal. (2006) 1248–1274.5. E. Burman, J. Guzm´an, and D. Leykekhman. Weighted error estimates of the con-tinuous interior penalty method for singularly perturbed problems. IMA J. Numer.Anal. (2009) 284–314.6. E. Burman and P. Hansbo. Edge stabilization for Galerkin approximations ofconvection-diffusion-reaction problems. Comput. Methods Appl. Mech. Engrg. (2004) 1437–1453.7. E. Burman and G. Smith. Analysis of the space semi-discretized SUPG method fortransient convection-diffusion equations.
Math. Models Methods Appl. Sci. (2011)2049–2068.8. R. Codina. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Comput. Methods Appl. Mech. Engrg. (2000)1579–1599.9. J. Douglas and T. Dupont. Interior penalty procedures for elliptic and parabolicGalerkin methods. In
Computing methods in applied sciences (Second Internat. Sym-pos., Versailles, 1975) , pages 207–216. (Lecture Notes in Phys., Vol. 58. Springer,Berlin, 1976). ovember 12, 2018 11:21 WSPC/INSTRUCTION FILEBurman˙final˙Paper01-13-100 E. Burman
10. J.-L. Guermond. Stabilization of Galerkin approximations of transport equations bysubgrid modeling.
M2AN Math. Model. Numer. Anal. , (1999) 1293–1316.11. J.-L. Guermond. Subgrid stabilization of Galerkin approximations of linear contrac-tion semi-groups of class C in Hilbert spaces. Numer. Methods Partial DifferentialEquations (2001) 1–25.12. J. Guzm´an. Local analysis of discontinuous Galerkin methods applied to singularlyperturbed problems. J. Numer. Math. (2006) 41–56.13. P. Henning and M. Ohlberger. The heterogeneous multiscale finite element method foradvection-diffusion problems with rapidly oscillating coefficients and large expecteddrift. Netw. Heterog. Media (2010) 711–7440.14. P. Houston, J. A. Mackenzie, E. S¨uli, and G. Warnecke. A posteriori error analysis fornumerical approximations of Friedrichs systems. Numer. Math. (1999) 433–470.15. C. Johnson, U. N¨avert, and J. Pitk¨aranta. Finite element methods for linear hyperbolicproblems. Comput. Methods Appl. Mech. Engrg. (1984) 285–312.16. C. Johnson and J. Pitk¨aranta. An analysis of the discontinuous Galerkin method fora scalar hyperbolic equation. Math. Comp. (1986) 1–26.17. M. Martins Afonso, A. Celani, R. Festa, and A. Mazzino. Large-eddy-simulation clo-sures of passive scalar turbulence: a systematic approach. J. Fluid Mech. (2003)355–364.18. G. Smith.
Global and local estimates for SUPG discretizations of the transientconvection–diffusion equation . PhD thesis, University of Sussex, 2013. in preparation.19. V. Thom´ee.
Galerkin finite element methods for parabolic problems , volume 25 of