Automatic Variationally Stable Analysis for Finite Element Computations: Transient Convection-Diffusion Problems
Eirik Valseth, Pouria Behnoudfar, Clint Dawson, Albert Romkes
AAutomatic Variationally Stable Analysis for Finite Element Computations:Transient Convection-Diffusion Problems
Eirik Valseth a, ∗ , Pouria Behnoudfar b , Clint Dawson a , Albert Romkes c a Oden Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, TX 78712, USA b Curtin University Kent Street, Bentley, Western Australia 6101, Australia c Department of Mechanical Engineering, South Dakota School of Mines & Technology, Rapid City, SD 57701, USA
Abstract
We establish stable finite element (FE) approximations of convection-diffusion initial boundary value problems (IB-VPs) using an unconditionally stable FE method, the automatic variationally stable finite element (AVS-FE) method [1].The transient convection-diffusion problem leads to issues in classical FE methods as the differential operator can beconsidered singular perturbation in both space and time. The unconditional stability of the AVS-FE method, regardlessof the underlying differential operator, allows us significant flexibility in the construction of FE approximations. Thus,in this paper, we take two distinct approaches to the FE discretization of the convection-diffusion problem: i ) consider-ing a space-time approach in which the temporal discretization is established using finite elements, and ii ) a method oflines approach in which we employ the AVS-FE method in space whereas the temporal domain is discretized using thegeneralized- α method. We also introduce another space-time technique in which the temporal direction is partitioned,thereby leading to finite space-time ”slices” to reduce the computational cost of the space-time discretizations. Inthe generalized- α method, we discretize the temporal domain into finite sized time-steps and adopt the generalized- α method as time integrator. Then, we derive a corresponding norm for the obtained operator to guarantee the temporalstability of the method.We present numerical verifications for both approaches, including numerical asymptotic convergence studies high-lighting optimal convergence properties. Furthermore, in the spirit of the discontinuous Petrov-Galerkin (DPG)method by Demkowicz and Gopalakrishnan [2–6], the AVS-FE method also leads to readily available a posteriori error estimates through a Riesz representer of the residual of the AVS-FE approximations. Hence, the norm of theresulting local restrictions of these estimates serve as error indicators in both space and time for which we presentmultiple numerical verifications in mesh adaptive strategies. Keywords: stability, discontinuous Petrov-Galerkin, method of lines, space-time finite element method, andadaptivity ∗ Corresponding author
Email addresses:
[email protected] (Eirik Valseth),
[email protected] (Pouria Behnoudfar),
[email protected] (Clint Dawson),
[email protected] (Albert Romkes)
Preprint submitted to Elsevier October 2, 2020 a r X i v : . [ m a t h . NA ] S e p . Introduction Transient BVPs are commonplace in engineering applications and to date still pose significant challenges in nu-merical analysis. Time dependency in many BVPs, such as the heat equation, involve partial derivatives of the trialvariable with respect to time and leads to numerical instabilities unless careful considerations are taken. The reasonbeing that the time derivative is a convective transport term, i.e., transient problems lead to unstable discretizations.Additionally, the target problem of convection-diffusion also result in numerical instabilities in its spatial discretiza-tions which lead to the development of the AVS-FE method in [1]. To overcome the stability issues in both spaceand time we propose two distinct approaches employing the AVS-FE method. First, we take a space-time approachin which space and time are discretized directly considering time an additional dimension using the AVS-FE method.Second, we consider a method of lines to decouple the computations in space and time and employ a generalized α method for the temporal discretization [7–9].The use of space-time FE methods remains attractive as the approximations are standard FE approximations andtherefore inherit attractive features of FE methods such as a priori and a posteriori error estimation and mesh adap-tive strategies. Examples of space-time FE methods can be found in, e.g., [10–12]. The AVS-FE method [1] beingunconditionally stable for any differential operator is therefore a prime candidate for space-time FE discretizations.The unconditional stability is a consequence of the philosophy of the DPG method in which the test space consistof functions that are computed on the fly from Riesz representation problems [2–6]. In [13], the AVS-FE method issuccessfully employed in space and time for the Cahn-Hilliard BVP. The goal here was the extension of the AVS-FEmethod to a nonlinear BVP as well as an initial verification of AVS-FE space-time solutions. Similarly, in [14], theAVS-FE method is employed for space time solutions of a nonlinear transient wave propagation problem, the Kortewegde-Vries equation. Furthermore, its built-in a posteriori error estimate and their corresponding error indicators canbe directly applied to drive adaptivity. The DPG method has been successfully applied to several transient problems,e.g., convection-diffusion and the Navier-Stokes equations [15–17]. These space-time formulations are available inthe DPG FE code Camellia of Nathan Roberts [18].Alternatively, the method of lines can be employed to decouple the discretization of space and time where thespatial dimension is discretized first to obtain a semi-discretized system. Then, using an appropriate method, i.e.,time integrator, the discretization of the temporal domain subsequently results in a fully discrete system of equations.Here, we employ the AVS-FE method in space and the generalized- α method in time. Chung and Hulbert introducedthe generalized- α method in [9] to solve hyperbolic problems and extended it to parabolic differential equations suchas Navier-Stokes equations in [19]. The method provides second-order accuracy in the temporal domain as well asunconditional stability. Although the method allows us to control the numerical dissipation in the high-frequencyregions, it delivers adequately accurate results in low-frequency domains. Introduction of a user-defined parameterprovides this control and includes the HHT- α method of Hilber, Hughes, Taylor [20] and the WBZ- α method ofWood, Bossak, and Zienkiewicz [21].In the following, we introduce the AVS-FE method for transient BVPs by taking the two distinct approachesintroduced above. In Section 2 we introduce our model problem and notations in addition to a review of the AVS-2E methodology and derivation of the AVS-FE weak formulation to be used. In this section we also present thediscretization of the weak form, a priori error estimates, the alternative saddle point structure of the AVS-FE method,and its built-in a posteriori error estimate. In Section 3 we present the time discretization techniques which we employ.The method of lines using AVS-FE method in space and generalized- α method in time is presented in Section 3.1; andspace-time AVS-FE method in Section 3.2. Results from numerical verifications for numerous PDEs and applicationsare presented in Section 4. Finally, we discuss conclusions and future work in Section 5.
2. The AVS-FE Methodology
The AVS-FE method [1] allows us to compute unconditionally stable FE approximations to BVPs for any differ-ential operator, provided its kernel is trivial. In this section we introduce our model problem and briefly review theAVS-FE method, a thorough introduction can be found in [1].
Let Ω ⊂ R N , N ≤ ∂ Ω and outward unit normal vector n ,and let T be the final time. Then, define Ω T = Ω × ( , T ) to be the space time domain which is open and bounded witha Lipschitz boundary ∂ Ω T = Γ in ∪ Γ out ∪ Γ ∪ Γ T . Γ in and Γ out are the in and outflow boundaries, respectively, and Γ and Γ T are the initial and final time boundaries, respectively. The transient model problem is therefore the followinglinear convection-diffusion IBVP:Find u such that: ∂ u ∂ t − ∇∇∇ · ( D ∇∇∇ u ) + b · ∇∇∇ u = f , in Ω T , u = u in , on Γ in , D ∇∇∇ u · n = g , on Γ out , u = u , on Γ , (1)where D denotes the second order diffusion tensor, with symmetric, bounded, and elliptic coefficients D i j ∈ L ∞ ( Ω ) ; b ∈ [ L ( Ω )] the convection coefficient; f ∈ L ( Ω ) the source function; and g ∈ H − / ( Γ out ) the Neumann boundarydata. Note that the gradient operator ∇∇∇ refers to the spatial gradient operator, e.g., ∇∇∇ ( · ) = { ∂ ( · ) ∂ x , ∂ ( · ) ∂ y } T . We omit the full derivation of the weak formulation here and mention key points only. The derivation of a weakformulation for the AVS-FE method is shown in, e.g. [1]. To establish a weak formulation of (1), we need a regularpartition P h of Ω T into elements K m , such that: Ω T = int ( (cid:91) K m ∈ P h K m ) .
3e introduce an auxiliary variable q = D ∇∇∇ u , and recast (1) as a system of (distributional) first-order PDEs:Find ( u , q ) ∈ H ( Ω T ) × H ( div , Ω ) such that: D ∇∇∇ u − q = , in Ω T , ∂ u ∂ t − ∇∇∇ · q + b · ∇∇∇ u = f , in Ω T , u = u in , on Γ in , q · n = g , on Γ out , u = u , on Γ . (2)Note that the flux variable q does not explicitly depend on time and therefore, in the weak enforcement of the PDE, itbelongs to H ( div , Ω ) and not H ( div , Ω T ) .To derive the AVS-FE weak formulation, we enforce the PDEs (2) weakly on each element K m ∈ P h , applyGreen’s identity to shift all derivatives to the test functions except the time derivative, and subsequent summation ofthe local contributions we arrive at the global variational formulation:Find ( u , q ) ∈ U ( Ω T ) such that: B (( u , q ) ; ( v , w )) = F (( v , w )) , ∀ ( v , w ) ∈ V ( P h ) , (3)In (3), the bilinear form, B : U ( Ω T ) × V ( P h ) −→ R , and linear functional, F : V ( P h ) −→ R , are defined: B (( u , q ) ; ( v . w )) def = ∑ K m ∈ P h (cid:26) (cid:90) K m (cid:20) − u D ∇∇∇ · w m − q · w m + ∂ u ∂ t v m + q · ∇∇∇ v m − ( b · ∇∇∇ v m ) u (cid:21) d x + (cid:73) ∂ K m \ ∂ Ω (cid:20) ( b · n ) γ m ( u ) γ m ( v m ) + γ m n ( w m ) γ m ( u ) − γ m n ( q ) γ m ( v m ) (cid:21) d s + (cid:73) ∂ K m \ Γ in (cid:20) ( b · n ) γ m ( u ) γ m ( v m ) + γ m n ( w m ) (cid:21) d s − (cid:73) ∂ K m \ Γ out (cid:20) γ m n ( q ) γ m ( v m ) (cid:21) d s (cid:27) , F (( v , w )) def = ∑ K m ∈ P h (cid:26) (cid:90) K m f v m d x + (cid:73) ∂ K m ∩ Γ out g γ m ( v m ) d s − (cid:73) ∂ K m ∩ Γ in (cid:20) u in γ m n ( w m ) + ( b · n ) u in γ m ( v m ) (cid:21) d s (cid:27) , (4)where the continuous trial and broken test function spaces, U ( Ω T ) and V ( P h ) , are defined as follows: U ( Ω T ) def = (cid:26) ( u , q ) ∈ H ( Ω T ) × H ( div , Ω ) : u Γ = u (cid:27) , V ( P h ) def = (cid:26) ( v , w ) ∈ H ( P h ) × H ( div , P h ) (cid:27) . (5)The broken Hilbert spaces are defined: H ( P h ) def = (cid:26) v ∈ L ( Ω T ) : v m ∈ H ( K m ) , ∀ K m ∈ P h (cid:27) , H ( div , P h ) def = (cid:26) v ∈ [ L ( Ω )] : w m ∈ H ( div , K m ) , ∀ K m ∈ P h (cid:27) , (6)and the norms on these spaces (cid:107)·(cid:107) U ( Ω T ) : U ( Ω T ) −→ [ , ∞ ) and (cid:107)·(cid:107) V ( P h ) : V ( P h ) −→ [ , ∞ ) are defined as follows: (cid:107) ( u , q ) (cid:107) U ( Ω T ) def = (cid:115) (cid:90) Ω (cid:20) ∇∇∇ u · ∇∇∇ u + u + ( ∇∇∇ · q ) + q · q (cid:21) d x . (cid:107) ( v , w ) (cid:107) V ( P h ) def = (cid:115) ∑ K m ∈ P h (cid:90) K m (cid:20) h m ∇∇∇ v m · ∇∇∇ v m + v m + h m ( ∇∇∇ · w m ) + w m · w m (cid:21) d x . (7)4he operators γ m : H ( K m ) : −→ H / ( ∂ K m ) and γ m n : H ( div , K m ) −→ H − / ( ∂ K m ) denote the trace and normal traceoperators (e.g., see [22]) on K m .The bilinear form and linear functional in (4) differs from the ones presented in [1] due to the term ∂ u ∂ t and theapplication of Green’s Identity to all terms involving spatial derivatives. The goal of this is to apply the boundaryconditions on both u and q weakly as indicated by the boundary integrals intersecting in and outflow portions of theglobal boundary in the linear functional. This choice has been made based on extensive numerical experimentationand the analysis of various AVS-FE weak formulations is left for future works.The weak formulation (3) represents a DPG formulation as the test and trial spaces are of different regularity.However, the fact that the trial space consist of global Hilbert spaces ensures the existence of the trace operatorsneeded on each element as well as the continuity of the trial spaces. Thus, we can employ classical FE approximationfunctions consisting of C ( Ω ) polynomials for u and, e.g, Raviart-Thomas polynomials for the flux q . In the followingwe review important key points of the AVS-FE method and present a general well-posedness result. Remark 2.1
The kernel of the differential operator of the underlying transient convection-diffusion problem is trivial.Hence, the solution of the AVS-FE weak formulation (3) is unique.
A key point in the well posedness of the AVS-FE weak formulation is the existence of an equivalent norm on the trialspace U ( Ω T ) . Since the kernel of B ( · , · ) is trivial, we introduce the following energy norm (cid:107)·(cid:107) B : U ( Ω T ) −→ [ , ∞ ) : (cid:107) ( u , q ) (cid:107) B def = sup ( v , w ) ∈ V ( P h ) \{ ( , ) } | B (( u , q ) ; ( v , w )) |(cid:107) ( v , w ) (cid:107) V ( P h ) . (8)As in the DPG method, the energy norm of ( u , q ) ∈ U ( Ω T ) can be identified by functions ( p , r ) ∈ V ( P h ) that aresolutions of the following Riesz representation problem: ( ( p , r ) ; ( v , w ) ) V ( P h ) = B ( ( u , q ) ; ( v , w ) ) , ∀ ( v , w ) ∈ V ( P h ) , (9)where ( · ; · ) V ( P h ) : V ( P h ) × V ( P h ) −→ R , is an inner product on V ( P h ) : ( ( r , z ) ; ( v , w ) ) V ( P h ) def = ∑ K m ∈ P h (cid:90) K m (cid:20) h m ∇∇∇ r m · ∇∇∇ v m + r m v m + h m ( ∇∇∇ · z m )( ∇∇∇ · w m ) + z m · w m (cid:21) d x . (10)Due to the Riesz representation problem (9) we can establish the equivalence between the energy norm of trial func-tions ( u , q ) ∈ U ( Ω T ) , and the norm of the Riesz representers ( p , r ) ∈ V ( P h ) : (cid:107) ( u , q ) (cid:107) B = (cid:107) ( p , r ) (cid:107) V ( P h ) . (11)Finally, the well posedness of the AVS-FE weak formulation in terms of the energy norm (8) is established by thefollowing lemma: Lemma 2.1
Let f ∈ ( H ( P h )) (cid:48) and g ∈ H − / ( Γ out ) . Then, the weak formulation (3) has a unique solution and iswell posed.Proof : Due to the energy norm definition and the Riesz problem (9), the application of Babuˇska Lax-Milgram Theo-rem [23] leads to inf-sup and continuity constants equal to unity. (cid:3) .3. AVS-FE Discretization In this section, we present a review of the AVS-FE discretization and for the sake of brevity we consider only spatialdiscretizations here. The discretization of the time domain is presented separately in Section 3. Hence, we suppress no-tation related to time dependency in this and the following section. To establish FE approximations ( u h , q h ) of ( u , q ) theAVS-FE method follows the classical FE method and represent the FE approximations u h and q h as linear combinationsof basis functions and their corresponding degree of freedom. In the case of Ω ⊂ R , ( e i ( x ) , ( E jx ( x ) , E ky ( x ))) ∈ U h ( Ω ) and { u hi ∈ R , i = , , . . . , N } , { q h , jx ∈ R , j = , , . . . , N } , and { q h , ky ∈ R , k = , , . . . , N } ; i.e., u h ( x ) = N ∑ i = u hi e i ( x ) , q hx ( x ) = N ∑ j = q h , jx E jx ( x ) , q hy ( x ) = N ∑ k = q h , ky E ky ( x ) . (12)Proper choices of bases are, e.g., continuous polynomials for the base variable u h ∈ P p ( Ω ) and Raviart-Thomas poly-nomials for the flux q h ∈ RT p ( Ω ) .As the test space V ( P h ) is broken, the test functions are to be piecewise discontinuous and are constructed byemploying the DPG philosophy [2–6, 24]. Hence, each basis function in the trial space U h ( Ω ) is paired with a (vectorvalued) test function. In the same spirit as ( p , r ) are the Riesz representers of ( u , q ) in (9), ( ˜ e i , ˜E i ) are the Rieszrepresenters of the basis functions ( e i , ( E jx ( x ) , E ky ( x ))) through (9), e.g, for a basis function e i for the scalar valued trialvariable: (cid:16) ( r , z ) ; ( ˜ e i , ˜E i ) (cid:17) V ( P h ) = B ( ( e i , ) ; ( r , z ) ) , ∀ ( r , z ) ∈ V ( P h ) , i = , . . . , N , (13)where the LHS is the inner product on V ( P h ) . Hence, the the Riesz representation problem, dubbed the test functionproblem, (13) are well posed. Clearly, (13) is of infinite dimension and must be approximated. Due to the brokennature of V ( P h ) we can solve local counterparts of (13) in a decoupled fashion element-by-element by computingpiecewise polynomial approximations of the same degree as the bases in the discrete trial space [1].Finally, the FE discretization of (3) governing the FE approximation ( u h , q h ) ∈ U h ( Ω T ) is:Find ( u h , q h ) ∈ U h ( Ω T ) such that: B (( u h , q h ) ; ( v ∗ , w ∗ )) = F (( v ∗ , w ∗ )) , ∀ ( v ∗ , w ∗ ) ∈ V ∗ ( P h ) , (14)where the finite dimensional subspace of test functions V ∗ ( P h ) ⊂ V ( P h ) is spanned by the numerical approximationsof the test functions. Thanks to the DPG methodology we employ to construct the optimal test space, the discreteproblem (14) is unconditionally stable for any mesh parameters h m and p m . The AVS-FE discretization (14) can be implemented in existing FE software by redefining routines that computethe element stiffness matrices. However, in several commonly used FE solvers, such as FEniCS [25] or Firedrake [26],manipulations of the element assembly routines may not as easily be performed. Thus, to enable straightforwardimplementation into these FE solvers, we will introduce an equivalent interpretation of the AVS-FE method as aglobal saddle point problem. This alternative interpretation is commonly employed for the DPG method [27] and is6lso the approach taken by Calo et. al. in [28]. We omit several details here and highlight only key features of thisinterpretation, interested readers are referred to [27] for a complete presentation.The AVS-FE method is a minimum residual method, in the sense that their solution realizes the minimum of afunctional according to the following principle: u h = arg min v h ∈ U h ( Ω T ) (cid:107) B v h − F (cid:107) V ( P h ) (cid:48) , (15)where B and F are operators induced by the bilinear and linear forms, respectively. Due to the Riesz representationproblem (9) and energy norm, we can relate the norm on the dual space V ( P h ) (cid:48) (cid:107)·(cid:107) V ( P h ) (cid:48) to the energy norm (cid:107)·(cid:107) B .Thus, we can consider a Riesz representer of the approximation error ( u − u h , q − q h ) , which we refer to as an errorrepresentation function [3]. This error representation function ( ˆ e , ˆ E ) is then defined as the solution of the followingweak problem: Find ( ˆ e , ˆ E ) ∈ V ( P h ) such that: (( ˆ e , ˆ E ) , ( v , w )) V ( P h ) = F ( v , w ) − B ( ( u h , q h ) ; ( v , w ) ) (cid:124) (cid:123)(cid:122) (cid:125) Residual ∀ ( v , w ) ∈ V ( P h ) . (16)Essentially, this error representation function is a solution of an AVS-FE weak form for the specific case in which theload/RHS is the residual in ( u h , q h ) .The energy norm of ( u − u h , q − q h ) can be identified by the V ( P h ) norm of the error representation function: Proposition 2.1
Let ( u , q ) ∈ U ( Ω T ) be the solution of the AVS-FE weak form (3) and ( u h , q h ) ∈ U h ( Ω T ) its corre-sponding AVS-FE approximation through (14) . Then, the the energy norm of ( u − u h , q − q h ) is identical to the V ( P h ) norm of ( ˆ e , ˆ E ) : (cid:107) ( u − u h , q − q h ) (cid:107) B = (cid:107) ( ˆ e , ˆ E ) (cid:107) V ( P h ) . (17) Proof : The identity is a consequence of the norm equivalence in (11), the definition of the energy norm (8) and theweak problem governing the error representation function (16). (cid:3)
The norm of approximate error representation function ( ˆ e h , ˆ E h ) is therefore an a posteriori error estimate, i.e, (cid:107) ( u − u h , q − q h ) (cid:107) B ≈ (cid:107) ( ˆ e h , ˆ E h ) (cid:107) V ( P h ) . (18)Furthermore, its local restriction can be computed element-wise as the space V ( P h ) is broken to yield the errorindicator: η = (cid:107) ( ˆ e h , ˆ E h ) (cid:107) V ( K m ) . (19)This type of error indicator has been applied with great success to multiple problems (see, e.g., [3, 6, 28, 29]), and weshow several numerical experiments using this indicator for the AVS-FE method in Section 4.The minimum residual interpretation allows us to consider the following AVS-FE saddle point formulation towhich we seek the solution ( u , q ) under the constraint of the error representation function minimizes the residual of7he AVS-FE method, see (16):Find ( u h , q h ) ∈ U h ( Ω T ) , ( ˆ e h , ˆ E h ) ∈ V h ( P h ) such that: (cid:0) ( ˆ e h , ˆ E h ) ; ( v , w ) (cid:1) V ( P h ) − B (( u h , q h ) ; ( v h , w h )) = − F ( v h , w h ) , ∀ ( v h , w h ) ∈ V h ( P h ) , B (( p h , r h ) ; ( ˆ e h , ˆ E h )) = , ∀ ( p h , r h ) ∈ U h ( Ω T ) . (20)Solving (20) gives both the AVS-FE solution for ( u h , q h ) and its error representation functions ( ˆ e h , ˆ E h ) in a singleglobal solution step. This is very convenient as we now have a built-in a posteriori error estimate and error indicatorsimmediately upon solving (20). However, the computational cost of doing so has been shifted from local computationsfor optimal test functions to the global cost of a larger system of equations. However, the global nature of (20)allows for very simple implementation of the AVS-FE method in readily available FE solvers like FEniCS [25] andFiredrake [26].
3. Time Discretization
In the weak formulation (3) we have made no assumptions on the type of discretization of the time domain. Dueto the flexibility of the AVS-FE method, we are lead to consider two distinct approaches. In both cases the spatialdiscretizations are performed with finite elements and the AVS-FE methodology. First, we consider a discretizationof the time domain by employing the method of lines to decouple the spatial and time discretization and subsequentlyemploying the generalized- α method. Second, the unconditional discrete stability of the AVS-FE method, allows usto discretize the time domain with finite elements without the use of a CFL condition, i.e., a space-time approach.Hence, giving us freedom in the choice of mesh parameters in our numerical experiments in the FE discretization ofthe time portion of the domain Ω T . In this section, we first discuss the method in an abstract setting before proceeding to the particular case of theAVS-FE method and generalized- α methods. To this end, we define two Hilbert spaces U and V , and introduce awell-posed weak formulation for a transient BVP, e.g., the convection-diffusion problem of Section 2.1:Find u ∈ U such that: b ( u , v ) = l ( v ) , ∀ v ∈ V , (21)where the bilinear form b contains all spatial and temporal terms. We then denote by L the time derivative operator,and modify the bilinear form to contain only spatial terms, i.e., b h . To seek approximations of (21) we consider FEpolynomial subspaces of U and V , i.e., U h and V h and introduce a semi-discrete formulation:Find u h ∈ U h such that: ( L ( u h ) , v h ) L ( Ω ) + b h ( u h , v h ) = l ( v h ) , ∀ v h ∈ V h , (22)8here ( · , · ) L ( Ω ) denotes the L ( Ω ) inner product. The semi-discrete formulation is well-posed, and the modifiedbilinear form b h satisfies the inf-sup condition with respect to the norm (cid:107)·(cid:107) V h .To advance the solution in time, we consider a uniform partition of the time domain from t = t N = T , where each time step t i is of width τ , and we compute approximations to u h at each step. In particular, weemploy the second-order accurate generalized- α methods of [9, 19]. For parabolic or first-order hyperbolic problems,the generalized- α method for the transient term L ( u h ) in (22) is to find u n + h ∈ U h , such that: ( ϑ n + α m h , v h ) L ( Ω ) + b h ( u n + α f h , v h ) = l n + α f ( , v h ) , ∀ v h ∈ V h , (23)where u nh , ϑ nh are the approximations to u ( ., t n ) and ∂ u ( ., t n ) ∂ t , respectively, and: l n + α f = l ( t n + α f ) , η n + α g = η n + α g δ ( η n ) , where η = u , ϑ , g = m , f , δ ( η n ) = η n + − η n . (24)By a Taylor expansion, we obtain u n + = u n + τϑ n + τγδ ( ϑ n ) as a linear combination of u n , ϑ n . Substitution of theexpressions in (24) into (23) gives for the LHS: ( ϑ n + h , v h ) L ( Ω ) + b h ( ζ t ϑ n + h , v h ) = ( α m l n + , v h ) , ∀ v h ∈ V h , (25)where ζ = τγα f α m , and the RHS is: l n + = l n + α f + ( α m − ) ( ϑ nh , v h ) L ( Ω ) + τα f ( γ − ) b h ( ϑ nh , v h ) − b h ( u nh , v h ) . (26)It can be shown that this scheme is formally second order accurate (see [19]) if: γ = + α m − α f . (27)Finally, to control high frequency damping, the two parameters α m and α f are defined in terms of the spectral radius ρ ∞ corresponding to an infinite time step: α m = (cid:16) − ρ ∞ + ρ ∞ (cid:17) , α f = + ρ ∞ . (28) Remark 3.1
The generalized- α method requires additional initial data for ϑ h . This value is obtained by setting α f = α m = n = and solving (23) .3.1.1. Generalized- α and The AVS-FE Method Having introduced the generalized- α method for a well defined weak formulation, we now extend it to the AVS-FEmethod for our model IBVP of convection-diffusion. Hence, let us consider the AVS-FE weak formulation (3), andthe trial and test spaces U ( Ω T ) and V ( P h ) (5). The generalized- α method for the AVS-FE method is:Find ( ϑ n + h , q n + h ) ∈ U h ( Ω ) such that: ( ϑ n + h , v h ) L ( Ω ) + b h (( ζ ϑ n + h , q n + h ) ; ( v h , w h )) = ( α m (cid:96) n + (( v h , w h )) , ∀ ( v h , w h ) ∈ V ∗ ( P h ) , (29)9here the operators are defined: b h (( u , q ) ; ( v . w )) def = ∑ K m ∈ P h (cid:26) (cid:90) K m (cid:20) − u D ∇∇∇ · w m − q · w m + q · ∇∇∇ v m − ( b · ∇∇∇ v m ) u (cid:21) d x + (cid:73) ∂ K m \ ∂ Ω (cid:20) ( b · n ) γ m ( u ) γ m ( v m ) + γ m n ( w m ) γ m ( u ) − γ m n ( q ) γ m ( v m ) (cid:21) d s + (cid:73) ∂ K m \ Γ in (cid:20) ( b · n ) γ m ( u ) γ m ( v m ) + γ m n ( w m ) (cid:21) d s − (cid:73) ∂ K m \ Γ out (cid:20) γ m n ( q ) γ m ( v m ) (cid:21) d s (cid:27) ,(cid:96) n + (( v h , w h )) def = ∑ K m ∈ P h (cid:26) ( f n + α f v h ) d x + (cid:73) ∂ K m ∩ Γ out g γ m ( v m ) d s − (cid:73) ∂ K m ∩ Γ in (cid:20) u in γ m n ( w m ) + ( b · n ) u in γ m ( v m ) (cid:21) d s (cid:27) + ( α m − ) ( ϑ nh , v h ) L ( Ω ) + τα f ( γ − ) · b h (( ϑ nh , ) ; ( v h , w h )) − b h (( u nh , q nh ) ; ( v h , w h )) . (30)To establish the solutions to (29) we shall take the same residual minimization approach introduced in Section 2.4and define a saddle point system similar to (20). The major difference between the ”original” weak form (3) and theone corresponding to the generalized- α method, i.e., (29) other than the adjusted bilinear and linear forms, is the term ( ϑ n + h , v h ) L ( Ω ) . The approximations of (29) are governed by the following residual minimization problem:Find ( ϑ n + h , q n + h ) ∈ U h ( Ω ) ⊂ U ( Ω ) , such that: ( ϑ n + h , q n + h ) = arg min z h ∈ U h ( Ω ) (cid:107) l n + − ( M + ζ B h ) z h (cid:107) V (cid:48) h , (31)where the operators B h and l n + correspond to the actions of the adjusted forms B h and (cid:96) n + , respectively, and M to the new term ( ϑ n + h , v h ) L ( Ω ) . Thankfully, the Riesz map (induced by the equivalent of the Riesz representationproblem (9) for (29)) allows us to relate the norm on the dual space (cid:107)·(cid:107) V (cid:48) h to the energy norm on U ( Ω ) exactly asin (11) and establish a saddle point system like (20). Hence, we define the following error representation function:Find ( ε n + h , ψ n + h ) ∈ V ( P h ) such that: (( ε n + h , ψ n + h ) , ( v , w ))) V ( P h ) = (cid:96) n + ( v , w ) − ( ϑ n + h , v ) L ( Ω ) + b h (( ζ ϑ n + h , q n + h ) ; ( v , w )) (cid:124) (cid:123)(cid:122) (cid:125) Residual , ∀ ( v , w ) ∈ V ( P h ) . (32)which now measures how far we are from the best approximation of ( ϑ n + h , q n + h ) at the current time step. In the samefashion as in Section 2.4, the norm of this function is an a posteriori error estimate and its restriction to each K m ∈ P h an error indicator. Thus, we can introduce the saddle point problem for each time step:Find ( ϑ n + h , q n + h ) ∈ U h ( Ω ) , ( ε n + h , ψ n + h ) ∈ V h ( P h ) such that: (( ε n + h , ψ n + h ) , ( v h , w h )) V h + (( ϑ n + h , ) , ( v h , w h )) L ( Ω ) + ζ · b h (( ϑ n + h , q n + h ) , ( v h , w h )) = α m (cid:96) n + (( v h , w h )) , ∀ ( v h , w h ) ∈ V h ( P h ) , (( z h , r h ) , ( ε n + h , )) L ( Ω ) + ζ · b h (( z h , r h ) , ( ε n + h , ψ n + h )) = , ∀ ( z h , r h ) ∈ U h ( Ω ) , (33)10here the inner product ( · , · ) V h is defined: (( ε n + h , ψ n + h ) , ( v h , w h )) V h = (( ζ · ε n + h , ψ n + h ) , ( v h , w h )) V ( P h ) + (( ε n + h , ) , ( v h , w h )) L ( Ω ) . (34)Computing ϑ n + h from (33), we readily obtain u n + h from a Taylor expansion and we solve this at each time step. Theoverall procedure only requires the inversion of a single matrix at each time step and two explicit updates. The defini-tion of the inner product (34) allows us to ensure a stabilized solution by a bound on the operator (( ϑ n + h , ) , ( v h , w h ))+ ζ · b h (( ϑ n + h , q n + h ) , ( v h , w h )) at each time step. Additionally, we maintain the consistency of problem which can bereadily checked by setting the time step: τ → As pointed out in Remark 3.1, we need to retrieve the additional initial data ϑ h to solve (33). Hence, we set α f = α m = ( ϑ h , q h ) ∈ U h ( Ω ) , ( ε h , ψ h ) ∈ V h ( P h ) such that: (( ε h , ψ h ) , ( v h , w h )) V h + (( ϑ h , ) , ( v h , w h )) L ( Ω ) = F (( v h , w h )) − ζ · b h (( u , q h ) , ( v h , w h )) , ∀ ( v h , w h ) ∈ V h ( P h ) , (( z h , r h ) , ( ε h , )) L ( Ω ) = , ∀ ( z h , r h ) ∈ U h ( Ω ) , (35)where u , q , and F (( v h , w h )) correspond to the initial data.To ascertain that (35) is well posed we propose a lemma: Lemma 3.1
Let ( v h , w h ) ∈ V h be arbitrary test functions. Then, ϑ h ∈ U h ( Ω ) exists and is unique. We omit the proof here as it is trivial to show that (( ϑ h , ) , ( v h , w h )) L ( Ω ) , i,e, the L ( Ω ) inner product, satisfies thefollowing three properties: • Stability: There exist a constant C sta > (cid:54) = z h ∈ U h ( Ω ) sup (cid:54) = v h ∈ V h | (( z h , ) , ( v h , w h )) L ( Ω ) |(cid:107) z h (cid:107) L ( Ω ) (cid:107) v h (cid:107) L ( Ω ) ≥ C sta . (36) • Consistency: Employing a similar argument as [28] to study the consistency of the saddle-point problem, wecan state the consistency as: (( ϑ h , ) , ( v h , w h )) L ( Ω ) = ( f , ( v h , w h )) − ζ · b h (( u h , q h ) , ( v h , w h )) , ∀ ( v h , w h ) ∈ V h (37) • Boundedness: There exists a constant C bnd < ∞ , uniformly with respect to the mesh size, such that: (( z , ) , ( v h , w h )) L ( Ω ) ≤ C bnd (cid:107) z (cid:107) L ( Ω ) (cid:107) v h (cid:107) L ( Ω ) , ∀ ( z , v h ) ∈ U × V h . (38)See [30] for details on these conditions. (cid:3) Thus, using (35), we have a stable and adaptive method to find the initial data which is critical for the generalized- α method to ensure second-order accuracy in time. 11 .2. Space-Time FE Approach The use of FE discretizations for transient problems is commonly avoided due to the inherently unstable natureof transient problems. The discretizations must be very carefully constructed to achieve discrete stability using theclassical FE method. Alternatively, conditionally stable FE methods can be applied. The latter leads to conditionallystable discretizations requiring arduous a priori analyses to ascertain stabilization parameters. However, the AVS-FEmethod is unconditionally stable for any choice of mesh parameters and allows us to discretize the entire space-timedomain with finite elements without the need for arduous stability analyses. Thus, a posteriori error estimates anderror indicators are immediately available to us as error indicators are obtained directly in the saddle point approachof the AVS-FE method (20).To establish AVS-FE space-time approximations of weak formulation (3) or (20), we pick appropriate discretiza-tions of the space H ( Ω T ) × H ( div , Ω ) . For H ( Ω T ) , the choice is classical (conforming) FE basis functions that are C continuous functions such as Lagrange or Legendre polynomials. For H ( div , Ω ) , which does not explicitly requirea three dimensional basis (due to the definition of the flux variable), a conforming choice of basis is, e.g., a Raviart-Thomas basis. However, as in [1], we employ approximations for q h by vector valued C ( Ω ) polynomials as this hasshown to yield superior results for convex domains. As an alternative to the space-time discretization of the full space-time domain Ω T , in this section we introducea time slice approach for the AVS-FE method. While the space-time approach introduced in the preceding sectionallows straightforward implementation of the AVS-FE method and its ”built-in” error indicator can drive mesh adaptiverefinements, the large number of degrees of freedom quickly makes the method intractable. In an effort to reduce thecomputational cost of the space-time approach, we propose to partition the space-time domain into ”space-time slices”.The slices can be constructed in a number of ways, from uniformly to a graded mesh structure as considered in [15, 31]for the DPG method.As solution information does not need to travel backwards in time, a solution can be obtained on a slice which canbe transferred to the neighboring slice as an initial condition. Hence, we can perform mesh adaptive refinements oneach slice to ensure the complete resolution of any interior or boundary layer (i.e., physical features) before proceedingto the next. This is of particular interest in applications in which physical parameters are time dependent leading towidely different solution features as time progresses. Hence, it gives an even greater flexibility in the choice of meshparameters than the ”global” space-time approach. In Figure 1 an arbitrary domain Ω T is shown and is partitioneduniformly into two space-time slices Ω ∩ ( , T slice ) and Ω ∩ ( , T ) . In Section 4 we compare two possible approachesemploying time slices and mesh adaptive refinements. A Priori
Error Estimates
In this section, we present a priori error estimates for the space-time version of the AVS-FE method in terms ofnorms of the approximation errors u − u h and q − q h . We present estimates in terms of the energy norm and classical12 x t Ω yT slice Figure 1: Partition of space-time domain into slices.
Sobolev norms on the trial space U ( Ω T ) . The a priori error estimates of the AVS-FE method for stationary convection-diffusion problems are established in [32] and form the basis of our results here. To keep this presentation brief, wedo not present extensive proofs but we present outlines and references to the appropriate literature. We assume herethat both q h and u h are discretized by C ( Ω ) polynomials and note that approximations using other bases for q h canbe established with minor efforts employing the work of Brezzi and Fortin [33].First we present the a priori estimate in terms of the energy norm. Proposition 3.1
Let ( u , q ) ∈ U ( Ω T ) be the exact solution of the AVS-FE weak formulation (3) and ( u h , q h ) ∈ U h ( Ω T ) its corresponding AVS-FE approximation. Then: ∃ C > (cid:107) ( u − u h , q − q h ) (cid:107) B ≤ C h µ − , (39) where h is the maximum element diameter, µ = min ( p u + , r ) , p u the minimum polynomial degree of approximationof u h in the mesh, and r the minimum regularity of the solution components ( p , r ) of the distributional PDE underlyingthe Riesz representation problem (9) .Proof : Because the AVS-FE approximation satisfies a best approximation property in terms of the energy norm (8),the norm equaivalence in (11), and the convergence properties of piecewise polynomial interpolants [34–36] used inthe FE approximation of the optimal test functions (9), we establish the bound (39). (cid:3) The energy norm cannot be computed exactly due to the supremum in its definition (8) and must be computed approx-imately by computing the approximate error representation function and Proposition 2.1. To establish bounds on theapproximation errors in terms of classical Sobolev norms on the trial space, we first note that since the energy norm isan equivalent norm on U ( Ω T ) , we have that: ∃ C , C ≥ C (cid:107) ( u , q ) (cid:107) U ( Ω T ) ≤ (cid:107) ( u , q ) (cid:107) B ≤ C (cid:107) ( u , q ) (cid:107) U ( Ω T ) ∀ ( u , q ) ∈ U ( Ω T ) . (40)Due to this norm equivalence, we can readily establish the following a priori estimates:13 emma 3.2 Let ( u , q ) ∈ U ( Ω T ) be the exact solution of the AVS-FE weak formulation (3) and ( u h , q h ) ∈ U h ( Ω T ) itscorresponding AVS-FE approximation (14) . Then: (cid:107) u − u h (cid:107) L ( Ω T ) ≤ C h µ , (cid:107) u − u h (cid:107) H ( Ω T ) ≤ C h µ − , (cid:107) q − q h (cid:107) H ( div , Ω ) ≤ C h µ − , (cid:107) ( u − u h , q − q h ) (cid:107) U ( Ω T ) ≤ C ( h µ − + h µ − ) . (41) where h is the maximum element diameter, µ = min ( p u + , r u ) , p u the minimum polynomial degree of approximationof u h in the mesh, µ = min ( p q + , r q ) , p q the minimum polynomial degree of approximation of q h , and r u , r q theregularities of the solutions u and q of the governing first order system of distributional PDEs (2) , respectively.Proof : Because the AVS-FE approximation satisfies a quasi-best approximation property in terms of (cid:107) ( · , · ) (cid:107) U ( Ω T ) dueto the norm equivalence on U ( Ω T ) , the best approximation property of the AVS-FE method (in terms of the energynorm), and the convergence properties of piecewise polynomial interpolants [34, 35], we establish the last boundin (41). The bounds on individual solution components in (41) are direct consequences of this last bound. The boundon the L ( Ω ) norm of u − u h in (41) is a direct consequence of an application of the the Aubin-Nitche lift [37, 38]. (cid:3) Remark 3.2
Establishing the optimal bound on (cid:107) u − u h (cid:107) L ( Ω T ) can be done using well established techniques found inFE literature. Unfortunately, a similar, optimal, bound cannot be established for (cid:107) q − q h (cid:107) L ( Ω ) as the required dualityresult for the flux variable can only be established for specially designed mesh geometries (see [39]). However, we donote that extensive numerical experimentation indicate optimal convergence behavior of the flux in the case of convexcomputational domains. For other approximation spaces, such as Raviart-Thomas, optimal bounds can be establishedusing the techniques found in the text of Brezzi and Fortin [33].
4. Numerical Verifications
To conduct numerical verifications, we consider the following simplified form of our model scalar-valued convec-tion diffusion problem (1) with Dirichlet boundary conditions: ∂ u ∂ t − ε ∆ u + b · ∇∇∇ u = f , in Ω , u = u , on ∂ Ω , u = u intial , on ∂ Ω ∩ { t = } , (42)where the coefficient ε ∈ L ∞ ( Ω ) is a scalar-valued isotropic diffusion coefficient. First, we verify the convergenceproperties of the AVS-FE method for both time discretization schemes in Section 4.1. The use of the time sliceapproach is also considered here, and we present verifications of two distinct refinement strategies for this approach.14n Section 4.2, we consider a hyperbolic problem of purely convective transport. Last, in Section 4.3 we presentverifications of a problem with both a hyperbolic and a parabolic part, i.e., a transient convection-diffusion problem.The particular case we investigate corresponds to a challenging physical application, a shock wave problem.In all the presented numerical experiments we use the saddle point description in (20) implemented in FEniCS [25].The verifications in which we employ adaptive refinements all use the same criterion as in [28], i.e., the built-in errorindicator (19) as well as the marking strategy of [40] using the approximate energy error computed using (17). Alsonote that in all cases where we report the number of degrees of freedom, we include the degrees of freedom for theerror representation function in the saddle point systems (20) and (33). The polynomial degree of approximation usedfor this error representation function is always identical to the degree of the trial space. To numerically verify the predicted rates of convergence from the a priori bounds above we perform numericalconvergence studies. This is accomplished by considering a well-known example of transient convection-diffusion,the Eriksson-Johnson problem [41]: ∂ u ∂ t − ε ∆ u + ∂ u ∂ x = f , in Ω T . (43)The exact solution of this problem is: u ex ( x ) = el t (cid:16) e λ x − e λ x (cid:17) + cos ( π y ) es x − er xe − s − e − r , (44)where l =
2, and: λ , = − ± √ − ε l − ε , r = + (cid:112) + π ε ε , s = − (cid:112) + π ε ε . (45)The problem domain Ω T = ( − , ) × ( − . , . ) × ( , ) . For these studies we consider the moderately convectiondominated case of (43) with ε = . N , which increases at O ( h − ) ,i.e., the h − convergence rates of the FE approximations can be extracted from these by a simple adjustment. For thelinear case of (cid:107) u − u h (cid:107) L ( Ω ) , we get O ( N − ) = O ( h ) = O ( h p + ) . Thus, the corresponding rates of convergence areas predicted by the a priori error estimates of Section 3.2.2 for (cid:107) u − u h (cid:107) L ( Ω ) and the energy norm in Figure 2. Theobserved rates for (cid:107) q − q h (cid:107) L ( Ω ) are as expected for the continuous polynomial approximations we use for linear poly-nomials approximations and is one half order lower than expected for the quadratic case. We expect the optimal rateto be recovered upon further mesh refinements but the computational cost is too large for our available computationalresources and we note that for other problems we have indeed observed the expected rates.Analogously, in Figure 3, the convergence plots for generalized- α are presented for to study the convergence ofthe method at the final time T = s with time step of τ = − . The observed rates of convergence in Figure 3 are15 -3 -2 -1 E rr o r n o r m DOFs (a) Linear polynomial approximations. -5 -4 -3 -2 -1 DOFs E rr o r n o r m (b) Quadratic polynomial approximations.Figure 2: Convergence histories for the space-time convergence study. as expected from the polynomial approximations employed. Comparison of the results in Figures 2 and 3 for thetwo methods reveal that the number of degrees of freedom is significantly larger for the space-time approach. This isbecause the discrete problem for the space-time method only requires a single global solve, whereas the generalized- α requires 1000 solves. The accuracy of the two methods is not straightforward to compare, as the errors reported inFigure 2 are global and the errors in Figure 3 are at the final time step. However, since the generalized- α we useis second order accurate, we expect for p ≥ In this section we present results for the time slice approach introduced in 3.2. In particular, we consider twodifferent approaches to mesh adaptive refinements for the time slice approach. We again consider the Eriksson-Johnson problem of the preceding section, now with ε = .
05. The time domain is partitioned into two slices from 0to 0 . s and from 0 . s to 1 . s , the polynomial degree of approximation is chosen to be one for both variables in bothapproaches.First, the space-time AVS-FE method is implemented on the bottom slice and we perform adaptive mesh refine-ments based on the built in error indicator (19) and we employ the marking strategy of [40] using the energy errorfrom the slice to ascertain which elements to refine. We perform a total of 13 steps of adaptation before the solution isused as in an initial condition for the final time slice to which we perform the same number of mesh refinements. Thesecond approach does not perform mesh refinements on each slice individually but rather computes the solution on thefirst slice and then moves directly to the second slice. Then, we again employ the marking strategy of [40] based on16 DOFs 10 − − − E rr o r no r m (a) Linear polynomial approximations, ρ ∞ = DOFs 10 − − − E rr o r no r m (b) Linear polynomial approximations, ρ ∞ = . DOFs 10 − − − − E rr o r no r m (c) Quadratic polynomial approximations, ρ ∞ = DOFs 10 − − − − E rr o r no r m (d) Quadratic polynomial approximations, ρ ∞ = . α method for time discretization at final time T = . s . the global energy error and we perform a total of 6 refinement steps.In Figure 4 we present the energy error and the L ( Ω ) error of both solution variables. The reported degrees offreedom in these figures are the maximum degrees of freedom for each slice at each refinement level. Comparisonbetween figures 4(a) and 4(b) reveal that the global error for both methods are of similar magnitude with the approachof adapting between the slices performing slightly better. However, we do not advocate one approach over the otherbut have demonstrated that both are capable of delivering accurate solutions and prospective users should consider17 OFs E rr o r n o r m -3 -2 -1 (a) Adaptivity between slices. DOFs E rr o r n o r m -3 -2 -1 (b) Adaptivity after both slices.Figure 4: Convergence histories for the time slice approach for the Eriksson-Johnson problem (43). both approaches. In this verification, we consider the problem of purely convective transport, i.e., ε = α method. Weconsider the spatial domain to be the unit square Ω = ( , ) ⊂ R , the convection vector is chosen to be variable, b = ( − y , x ) T . The source, boundary conditions are such that the exact solution u ex is: u ex ( x , y , t ) = t (cid:16) + tanh (cid:16) M (cid:16) . − (cid:12)(cid:12)(cid:12) . − (cid:112) x + y (cid:12)(cid:12)(cid:12)(cid:17)(cid:17)(cid:17) , (46)and the initial solution is u =
0. We choose the time-step τ = × − and the final time T final =
1. We set theparameter M =
500 in (46) and ρ ∞ = . α method and employ linear polynomial elements.To show the effectiveness of the built-in error estimate of the AVS-FE method, we again employ a mesh adaptiverefinement strategy. The initial mesh is coarse and is such that no elements in the mesh align with the expected interiorlayers of (46) and is shown in Figure 5(a). In Figure 6(b), the converged solution at the final time T = . s is presentedand in Figure 6(a) the corresponding convergence history is shown. In Figure 5(b), we present the final adapted mesh.Clearly, the error representation function provides accurate error indicators that are highly effective leading to meshrefinements focused near the circular interior layer. As a final numerical verification, we present a consideration of (42) in which the solution behaves as two shockstraveling through the space-time domain while rotating about the origin. Furthermore, the choices we make for the18 a) Initial mesh (b) Final meshFigure 5: Adaptive mesh guided using the obtained error representation. DOFs 10 − − E rr o r no r m (a) L and energy errors. (b) Solution u h at final time.Figure 6: Estimated solution and convergence histories at the final time T = . s for the pure convection problem. problem parameters are such that the interface of the shock is skewed and rotates in the space-time domain as t → T final .19hus, we have the following choices: b = {− x + y , } T , u = , ε = − , u = , f = − x ε + x ( − y ) . T final = . s Ω = ( − , ) × ( − , ) (47)For this particular problem, we present the time slice approach in which we perform mesh adaptations between eachslice. Experience has shown that the slice containing the initial condition is critical to the proper resolution of thespace-time process. Thus, we consider the case of two space-time slices, the first from 0 s to 0 . s and the final from0 . s to 1 . s . In Figures 7 and 8 we present the AVS-FE solution for the base variable at different time steps. Asexpected, two shock-waves originate at the boundaries of x = ±
1, and as time progress, the two waves approachthe center of the domain while rotating. The adaptively refined meshes shown in Figures 7(a) and 8(b) at the finaltimes of each slice show that the mesh refinements are focused at the interfaces of the shocks, further indicating theapplicability of the built-in error indicators. (a) Solution u h at t = . s with final adapted mesh. (b) Solution u h at t = . s .Figure 7: AVS-FE approximations of the shock problem, i.e., (42) with parameters from 47.
5. Conclusions
The AVS-FE method is a Petrov-Galerkin method which uses classical continuous FE trial basis functions, whilethe test space consist of functions that are discontinuous across element edges. This broken topology in the test spaceallows us to employ the DPG philosophy and introduce an equivalent saddle point problem which we implement usinghigh level FE solvers. We have introduced two distinct approaches to transient problems using the AVS-FE method.20 a) Solution u h at t = . s . (b) Solution u h at t = . s with final adapted mesh.Figure 8: AVS-FE approximations of the shock problem, i.e., (42) with parameters from 47. First, we take a space-time approach in which the entire space-time domain in discretized using finite elements, andsecond, using the method of lines to discretize the spatial domain independently to deliver a semi-discretized system.Then, using a time-marching method, we obtain a fully discrete system.The space-time method allows us to exploit the unconditional stability of the AVS-FE method and perform asingle global solve governing the FE approximation. As the AVS-FE approximations computed from the saddle pointsystem (31) come with built-in error indicators, we are capable of utilizing mesh adaptive strategies in space and time.For the convection-diffusion IBVP we consider here, we establish a priori error estimates to the space-time AVS-FEmethod in Section 3.2.2 as for the classical Galerkin FE method due to the best approximation property of the AVS-FEmethod. In an effort to reduce the computational cost of solving the global system we have consider a time sliceapproach in which the space-time domain is partitioned into finite sized space-time slices on which we employ theAVS-FE method. The advantage here is that the size of the global system is reduced and we are able to employ meshadaptive strategies on each slice.The method of lines, in which we use the AVS-FE method for the spatial discretization and a generalized- α method to derive a fully-discretized system. In this case, the discrete stability in the temporal domain is ensured bythe generalized- α method leading to highly efficient stable FE computations. We show that the AVS-FE method usesa corresponding norm as a function of the time-step. Another distinguishing feature of this method is that due tothe influence of the initial data on the accuracy of the solution, we find a stable approximation for ∂ u ∂ t at the initialtime. Accordingly, at each time step, one requires to solve a system with a smaller number of degrees of freedom incomparison with the space-time approach.Numerical verifications for several cases of the transient convection-diffusion IBVP show that both methods exhibitoptimal asymptotic convergence behavior as well as similar norms of the numerical approximation error. For degreesof approximation above 2, the space-time approach becomes more accurate as it is not limited to the second-order21ccuracy of the generalized- α method. However, we do not advocate one method over the other but we point out thesedifferences for potential users as their available computational resources will likely dictate which approach to use. Forboth cases, we present additional numerical verifiactions highlighting the adaptive mesh refinement capabilities. Infuture efforts, we expect to pursue alternative error estimators and indicators as well as the AVS-FE approximation ofchallenging transient physical phenomena. Acknowledgements
Authors Valseth and Dawson have been supported by the United States National Science Foundation - NSFPREEVENTS Track 2 Program,under NSF Grant Number 1855047. Authors Valseth and Romkes have been sup-ported by the United States National Science Foundation - NSF CBET Program, under NSF Grant Number 1805550.The authors also gratefully acknowledge the assistance of Austin Kaul of the Department of Mechanical Engineeringat South Dakota School of Mines and Technology in performing the numerical verifications of Section 4.1 and 4.1.1.
References [1] V. M. Calo, A. Romkes, E. Valseth, Automatic variationally stable analysis for FE computations: an introduction,in: Barrenechea G., Mackenzie J. (eds) Boundary and Interior Layers, Computational and Asymptotic MethodsBAIL 2018, Springer, 2020, pp. 19–43.[2] L. Demkowicz, J. Gopalakrishnan, A class of discontinuous Petrov-Galerkin methods. Part I: The transportequation, Computer Methods in Applied Mechanics and Engineering 199 (23) (2010) 1558–1572.[3] C. Carstensen, L. Demkowicz, J. Gopalakrishnan, A posteriori error control for DPG methods, SIAM Journal onNumerical Analysis 52 (3) (2014) 1335–1353.[4] L. Demkowicz, J. Gopalakrishnan, Analysis of the DPG method for the Poisson equation, SIAM Journal onNumerical Analysis 49 (5) (2011) 1788–1809.[5] L. Demkowicz, J. Gopalakrishnan, A class of discontinuous Petrov-Galerkin methods. II. Optimal test functions,Numerical Methods for Partial Differential Equations 27 (1) (2011) 70–105.[6] L. Demkowicz, J. Gopalakrishnan, A class of discontinuous Petrov-Galerkin methods. Part III: Adaptivity, Ap-plied numerical mathematics 62 (4) (2012) 396–427.[7] Q. Deng, P. Behnoudfar, V. M. Calo, High-order generalized- α methods, arXiv preprint arXiv:1902.05253(2019).[8] P. Behnoudfar, Q. Deng, V. M. Calo, Higher-order generalized- α methods for hyperbolic problems, arXivpreprint arXiv:1906.06081 (2019). 229] J. Chung, G. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation:the generalized- α method, Journal of Applied Mechanics 60 (2) (1993) 371–375.[10] T. J. R. Hughes, J. R. Stewart, A space-time formulation for multiscale phenomena, Journal of Computationaland Applied Mathematics 74 (1996) 217–229.[11] T. J. Hughes, G. M. Hulbert, Space-time finite element methods for elastodynamics: formulations and errorestimates, Computer methods in applied mechanics and engineering 66 (3) (1988) 339–363.[12] A. K. Aziz, P. Monk, Continuous finite elements in space and time for the heat equation, Mathematics of Com-putation 52 (186) (1989) 255–274.[13] E. Valseth, A. Romkes, A. R. Kaul, A stable FE method for the space-time solution of the cahn-hilliard equation,Preprint Submitted to Journal of Computational Physics, arXiv preprint arXiv:2006.02283 (2020).[14] E. Valseth, C. Dawson, An unconditionally stable spacetime FE method for the Kortewegde Vries equation,Computer Methods in Applied Mechanics and Engineering 371 (2020) 113297. doi:https://doi.org/10.1016/j.cma.2020.113297 .[15] T. E. Ellis, L. Demkowicz, J. Chan, R. D. Moser, Space-time DPG: Designing a method for massively parallelCFD, ICES report, The Institute for Computational Engineering and Sciences, The University of Texas at Austin(2014) 14–32.[16] T. Ellis, J. Chan, L. Demkowicz, Robust DPG methods for transient convection-diffusion, in: Building bridges:connections and challenges in modern approaches to numerical partial differential equations, Springer, 2016, pp.179–203.[17] N. V. Roberts, L. Demkowicz, R. Moser, A discontinuous Petrov–Galerkin methodology for adaptive solutionsto the incompressible Navier–Stokes equations, Journal of Computational Physics 301 (2015) 456–483.[18] N. V. Roberts, Camellia: A software framework for discontinuous Petrov–Galerkin methods, Computers & Math-ematics with Applications 68 (11) (2014) 1581–1604.[19] K. E. Jansen, C. H. Whiting, G. M. Hulbert, A generalized- α method for integrating the filtered Navier–Stokesequations with a stabilized finite element method, Computer Methods in Applied Mechanics and Engineering190 (3-4) (2000) 305–319.[20] H. M. Hilber, T. J. Hughes, R. L. Taylor, Improved numerical dissipation for time integration algorithms instructural dynamics, Earthquake Engineering & Structural Dynamics 5 (3) (1977) 283–292.[21] W. Wood, M. Bossak, O. Zienkiewicz, An alpha modification of Newmark’s method, International Journal forNumerical Methods in Engineering 15 (10) (1980) 1562–1566.2322] V. Girault, P.-A. Raviart, Finite element methods for Navier-Stokes equations; theory and algorithms, in: SpringerSeries in Computational Mathematics, Vol. 5, Springer-Verlag, 1986.[23] I. Babuˇska, Error-bounds for finite element method, Numerische Mathematik 16 (1971) 322–333.[24] A. H. Niemi, N. O. Collier, V. M. Calo, Automatically stable discontinuous Petrov–Galerkin methods for sta-tionary transport problems: Quasi-optimal test space norm, Computers & Mathematics with Applications 66 (10)(2013) 2096–2113.[25] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N.Wells, The fenics project version 1.5, Archive of Numerical Software 3 (100) (2015) 9–23.[26] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. McRae, G.-T. Bercea, G. R. Markall,P. H. Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Transactions onMathematical Software (TOMS) 43 (3) (2017) 24.[27] L. F. Demkowicz, J. Gopalakrishnan, An overview of the discontinuous Petrov-Galerkin method, in: Recentdevelopments in discontinuous Galerkin finite element methods for partial differential equations, Springer, 2014,pp. 149–180.[28] V. M. Calo, A. Ern, I. Muga, S. Rojas, An adaptive stabilized conforming finite element method via residualminimization on dual discontinuous Galerkin norms, Computer Methods in Applied Mechanics and Engineering363 (2020) 112891.[29] F. Fuentes, B. Keith, L. Demkowicz, P. Le Tallec, Coupled variational formulations of linear elasticity and theDPG methodology, Journal of Computational Physics 348 (2017) 715–731.[30] D. A. Di Pietro, A. Ern, Mathematical aspects of discontinuous Galerkin methods, Vol. 69, Springer Science &Business Media, 2011.[31] T. E. Ellis, Space-time discontinuous Petrov-Galerkin finite elements for transient fluid mechanics, Ph.D. thesis(2016).[32] E. Valseth, Automatic Variationally Stable Analysis for Finite Element Computations, Ph.D. thesis (2019).[33] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Vol. 15, Springer-Verlag, 1991.[34] I. Babuˇska, M. Suri, The hphp