Automatic Variationally Stable Analysis for Finite Element Computations: Transient Convection-Diffusion Problems

We establish stable finite element (FE) approximations of convection-diffusion initial boundary value problems using the automatic variationally stable finite element (AVS-FE) method. The transient convection-diffusion problem leads to issues in classical FE methods as the differential operator can be considered singular perturbation in both space and time. The unconditional stability of the AVS-FE method, regardless of the underlying differential operator, allows us significant flexibility in the construction of FE approximations. We take two distinct approaches to the FE discretization of the convection-diffusion problem: i) considering a space-time approach in which the temporal discretization is established using finite elements, and ii) a method of lines approach in which we employ the AVS-FE method in space whereas the temporal domain is discretized using the generalized-alpha method. In the generalized-alpha method, we discretize the temporal domain into finite sized time-steps and adopt the generalized-alpha method as time integrator. Then, we derive a corresponding norm for the obtained operator to guarantee the temporal stability of the method. We present numerical verifications for both approaches, including numerical asymptotic convergence studies highlighting optimal convergence properties. Furthermore, in the spirit of the discontinuous Petrov-Galerkin method by Demkowicz and Gopalakrishnan, 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 the resulting local restrictions of these estimates serve as error indicators in both space and time for which we present multiple numerical verifications adaptive strategies.


Introduction
Transient BVPs are commonplace in engineering applications and to date still pose significant challenges in numerical analysis. Time dependency in many BVPs, such as the heat equation, involve partial derivatives of the trial variable with respect to time and leads to numerical instabilities unless careful considerations are taken. The reason being 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 discretizations which lead to the development of the AVS-FE method in [1]. To overcome the stability issues in both space and time we propose two distinct approaches employing the AVS-FE method. First, we take a space-time approach in 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][8][9].
The use of space-time FE methods remains attractive as the approximations are standard FE approximations and therefore inherit attractive features of FE methods such as a priori and a posteriori error estimation and mesh adaptive strategies. Examples of space-time FE methods can be found in, e.g., [10][11][12]. The AVS-FE method [1] being unconditionally 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 consist of functions that are computed on the fly from Riesz representation problems [2][3][4][5][6]. In [13], the AVS-FE method is successfully employed in space and time for the Cahn-Hilliard BVP. The goal here was the extension of the AVS-FE method to a nonlinear BVP as well as an initial verification of AVS-FE space-time solutions. Similarly, in [14], the AVS-FE method is employed for space time solutions of a nonlinear transient wave propagation problem, the Korteweg de-Vries equation. Furthermore, its built-in a posteriori error estimate and their corresponding error indicators can be 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][16][17]. These space-time formulations are available in the 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 the spatial 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 introduced the generalized-α method in [9] to solve hyperbolic problems and extended it to parabolic differential equations such as Navier-Stokes equations in [19]. The method provides second-order accuracy in the temporal domain as well as unconditional stability. Although the method allows us to control the numerical dissipation in the high-frequency regions, it delivers adequately accurate results in low-frequency domains. Introduction of a user-defined parameter provides this control and includes the HHT-α method of Hilber, Hughes, Taylor [20] and the WBZ-α method of Wood, Bossak, and Zienkiewicz [21].
In the following, we introduce the AVS-FE method for transient BVPs by taking the two distinct approaches introduced above. In Section 2 we introduce our model problem and notations in addition to a review of the AVS-FE methodology and derivation of the AVS-FE weak formulation to be used. In this section we also present the discretization 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; and space-time AVS-FE method in Section 3.2. Results from numerical verifications for numerous PDEs and applications are presented in Section 4. Finally, we discuss conclusions and future work in Section 5.

The AVS-FE Methodology
The AVS-FE method [1] allows us to compute unconditionally stable FE approximations to BVPs for any differential operator, provided its kernel is trivial. In this section we introduce our model problem and briefly review the AVS-FE method, a thorough introduction can be found in [1].

Model Problem and Notation
Let Ω ⊂ R N , N ≤ 2 be an open bounded domain with Lipschitz boundary ∂ Ω and outward unit normal vector n, and let T be the final time. Then, define Ω T = Ω × (0, T ) to be the space time domain which is open and bounded with a Lipschitz boundary ∂ Ω T = Γ in ∪ Γ out ∪ Γ 0 ∪ Γ T . Γ in and Γ out are the in and outflow boundaries, respectively, and Γ 0 and Γ T are the initial and final time boundaries, respectively. The transient model problem is therefore the following linear convection-diffusion IBVP: Find u such that: where D denotes the second order diffusion tensor, with symmetric, bounded, and elliptic coefficients D i j ∈ L ∞ (Ω); b ∈ [L 2 (Ω)] 2 the convection coefficient; f ∈ L 2 (Ω) the source function; and g ∈ H −1/2 (Γ out ) the Neumann boundary data. Note that the gradient operator ∇ ∇ ∇ refers to the spatial gradient operator, e.g.,

Weak Formulation
We omit the full derivation of the weak formulation here and mention key points only. The derivation of a weak formulation for the AVS-FE method is shown in, e.g. [1]. To establish a weak formulation of (1), we need a regular partition P h of Ω T into elements K m , such that: We introduce an auxiliary variable q = D∇ ∇ ∇u, and recast (1) as a system of (distributional) first-order PDEs: Find (u, q) ∈ H 1 (Ω T ) × H(div, Ω) such that: Note that the flux variable q does not explicitly depend on time and therefore, in the weak enforcement of the PDE, it belongs 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 , apply Green's identity to shift all derivatives to the test functions except the time derivative, and subsequent summation of the local contributions we arrive at the global variational formulation: Find (u, q) ∈ U(Ω T ) such that: In (3), the bilinear form, B : U(Ω T ) ×V (P h ) −→ R, and linear functional, F : V (P h ) −→ R, are defined: where the continuous trial and broken test function spaces, U(Ω T ) and V (P h ), are defined as follows: The broken Hilbert spaces are defined: and the norms on these spaces · U(Ω T ) : U(Ω T )−→[0, ∞) and · V (P h ) : V (P h )−→[0, ∞) are defined as follows: The operators γ m 0 : H 1 (K m ) :−→ H 1/2 (∂ K m ) and γ m n : H(div, K m ) −→ H −1/2 (∂ K m ) denote the trace and normal trace operators (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 the application of Green's Identity to all terms involving spatial derivatives. The goal of this is to apply the boundary conditions on both u and q weakly as indicated by the boundary integrals intersecting in and outflow portions of the global boundary in the linear functional. This choice has been made based on extensive numerical experimentation and 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 operators needed on each element as well as the continuity of the trial spaces. Thus, we can employ classical FE approximation functions consisting of C 0 (Ω) polynomials for u and, e.g, Raviart-Thomas polynomials for the flux q. In the following we 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 trial space U(Ω T ). Since the kernel of B(·, ·) is trivial, we introduce the following energy norm · B : U(Ω T ) −→ [0, ∞): As in the DPG method, the energy norm of (u, q) ∈ U(Ω T ) can be identified by functions (p, r) ∈ V (P h ) that are solutions of the following Riesz representation problem: where ( ·; · ) V (P h ) : V (P h ) ×V (P h ) −→ R, is an inner product on V (P h ): Due to the Riesz representation problem (9) we can establish the equivalence between the energy norm of trial functions (u, q) ∈ U(Ω T ), and the norm of the Riesz representers (p, r) ∈ V (P h ): Finally, the well posedness of the AVS-FE weak formulation in terms of the energy norm (8) is established by the following lemma: Lemma 2.1 Let f ∈ (H 1 (P h )) and g ∈ H −1/2 (Γ out ). Then, the weak formulation (3) has a unique solution and is well posed.
Proof : Due to the energy norm definition and the Riesz problem (9), the application of Babuška Lax-Milgram Theorem [23] leads to inf-sup and continuity constants equal to unity.

AVS-FE Discretization
In this section, we present a review of the AVS-FE discretization and for the sake of brevity we consider only spatial discretizations here. The discretization of the time domain is presented separately in Section 3. Hence, we suppress notation related to time dependency in this and the following section. To establish FE approximations (u h , q h ) of (u, q) the AVS-FE method follows the classical FE method and represent the FE approximations u h and q h as linear combinations of basis functions and their corresponding degree of freedom. In the case of Ω ⊂ R 2 , (e i (x), (E j x (x), E k y (x))) ∈ U h (Ω) and {u h i ∈ R, i = 1, 2, . . . , N}, {q h, j x ∈ R, j = 1, 2, . . . , N}, and {q h,k y ∈ R, k = 1, 2, . . . , N}; i.e., Proper choices of bases are, e.g., continuous polynomials for the base variable u h ∈ P p (Ω) and Raviart-Thomas polynomials 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 by employing the DPG philosophy [2][3][4][5][6]24]. Hence, each basis function in the trial space U h (Ω) is paired with a (vector valued) test function. In the same spirit as (p, r) are the Riesz representers of (u, q) in (9), (9), e.g, for a basis function e i for the scalar valued trial variable: where the LHS is the inner product on V (P h ). Hence, the the Riesz representation problem, dubbed the test function problem, (13) are well posed. Clearly, (13) is of infinite dimension and must be approximated. Due to the broken nature of V (P h ) we can solve local counterparts of (13) in a decoupled fashion element-by-element by computing piecewise 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: where the finite dimensional subspace of test functions V * (P h ) ⊂ V (P h ) is spanned by the numerical approximations of the test functions. Thanks to the DPG methodology we employ to construct the optimal test space, the discrete problem (14) is unconditionally stable for any mesh parameters h m and p m .

Saddle Point Problem
The AVS-FE discretization (14) can be implemented in existing FE software by redefining routines that compute the 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 straightforward implementation into these FE solvers, we will introduce an equivalent interpretation of the AVS-FE method as a global saddle point problem. This alternative interpretation is commonly employed for the DPG method [27] and is also the approach taken by Calo et. al. in [28]. We omit several details here and highlight only key features of this interpretation, 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 a functional according to the following principle: where B and F are operators induced by the bilinear and linear forms, respectively. Due to the Riesz representation problem (9) and energy norm, we can relate the norm on the dual space Thus, we can consider a Riesz representer of the approximation error (u − u h , q − q h ), which we refer to as an error representation function [3]. This error representation function (ê,Ê) is then defined as the solution of the following weak problem: Essentially, this error representation function is a solution of an AVS-FE weak form for the specific case in which the The energy norm of (u − u h , q − q h ) can be identified by the V (P h ) norm of the error representation function: sponding AVS-FE approximation through (14). Then, the the energy norm of (u norm of (ê,Ê): Proof : The identity is a consequence of the norm equivalence in (11), the definition of the energy norm (8) and the weak problem governing the error representation function (16).
The norm of approximate error representation function (ê h ,Ê h ) is therefore an a posteriori error estimate, i.e, Furthermore, its local restriction can be computed element-wise as the space V (P h ) is broken to yield the error indicator: This type of error indicator has been applied with great success to multiple problems (see, e.g., [3,6,28,29]), and we show 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 to which we seek the solution (u, q) under the constraint of the error representation function minimizes the residual of the AVS-FE method, see (16): Solving (20) gives both the AVS-FE solution for (u h , q h ) and its error representation functions (ê h ,Ê h ) in a single global solution step. This is very convenient as we now have a built-in a posteriori error estimate and error indicators immediately upon solving (20). However, the computational cost of doing so has been shifted from local computations for 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] and Firedrake [26].

Time Discretization
In the weak formulation (3)

Method of Lines
In this section, we first discuss the method in an abstract setting before proceeding to the particular case of the AVS-FE method and generalized-α methods. To this end, we define two Hilbert spaces U and V, and introduce a well-posed weak formulation for a transient BVP, e.g., the convection-diffusion problem of Section 2.1: 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 FE polynomial 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: where (·, ·) L 2 (Ω) denotes the L 2 (Ω) inner product. The semi-discrete formulation is well-posed, and the modified bilinear form b h satisfies the inf-sup condition with respect to the norm · V h .
To advance the solution in time, we consider a uniform partition of the time domain from t 0 = 0 to the final time t N = T , where each time step t i is of width τ, and we compute approximations to u h at each step. In particular, we employ 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+1 h ∈ U h , such that: where u n h , ϑ n h are the approximations to u(.,t n ) and ∂ u(.,t n ) ∂t , respectively, and: By a Taylor expansion, we obtain u n+1 = u n + τϑ n + τγδ (ϑ n ) as a linear combination of u n , ϑ n . Substitution of the expressions in (24) into (23) gives for the LHS: where ζ = τγα f α m , and the RHS is: It can be shown that this scheme is formally second order accurate (see [19]) if: 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: Remark 3.1 The generalized-α method requires additional initial data for ϑ 0 h . This value is obtained by setting α f = α m = n = 0 and solving (23).

Generalized-α and The AVS-FE Method
Having introduced the generalized-α method for a well defined weak formulation, we now extend it to the AVS-FE method for our model IBVP of convection-diffusion. Hence, let us consider the AVS-FE weak formulation (3), and the trial and test spaces U(Ω T ) and V (P h ) (5). The generalized-α method for the AVS-FE method is: (Ω) such that: where the operators are defined: To establish the solutions to (29) we shall take the same residual minimization approach introduced in Section 2.4 and define a saddle point system similar to (20). The major difference between the "original" weak form (3) and the one corresponding to the generalized-α method, i.e., (29) other than the adjusted bilinear and linear forms, is the term . The approximations of (29) are governed by the following residual minimization problem: where the operators B h and l n+1 correspond to the actions of the adjusted forms B h and n+1 , respectively, and M to the new term (ϑ n+1 h , v h ) L 2 (Ω) . Thankfully, the Riesz map (induced by the equivalent of the Riesz representation problem (9) for (29)) allows us to relate the norm on the dual space · V h to the energy norm on U(Ω) exactly as in (11) and establish a saddle point system like (20). Hence, we define the following error representation function: which now measures how far we are from the best approximation of (ϑ n+1 h , q n+1 h ) at the current time step. In the same fashion 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: where the inner product (·, ·) V h is defined: Computing ϑ n+1 h from (33), we readily obtain u n+1 h from a Taylor expansion and we solve this at each time step. The overall procedure only requires the inversion of a single matrix at each time step and two explicit updates. The definition of the inner product (34) allows us to ensure a stabilized solution by a bound on the operator ((ϑ n+1 at each time step. Additionally, we maintain the consistency of problem which can be readily checked by setting the time step: τ → 0.

Retrieving initial data
As pointed out in Remark 3.1, we need to retrieve the additional initial data ϑ 0 h to solve (33). Hence, we set α f = α m = 0 and get: where u 0 , q 0 , and F 0 ((v h , w h )) correspond to the initial data.
To ascertain that (35) is well posed we propose a lemma: be arbitrary test functions. Then, ϑ 0 h ∈ U h (Ω) exists and is unique.
We omit the proof here as it is trivial to show that ((ϑ 0 h , 0), (v h , w h )) L 2 (Ω) , i,e, the L 2 (Ω) inner product, satisfies the following three properties: • Stability: There exist a constant C sta > 0 independent of the mesh size, such that: • Consistency: Employing a similar argument as [28] to study the consistency of the saddle-point problem, we can state the consistency as: • Boundedness: There exists a constant C bnd < ∞, uniformly with respect to the mesh size, such that: See [30] for details on these conditions.
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.

Space-Time FE Approach
The use of FE discretizations for transient problems is commonly avoided due to the inherently unstable nature of transient problems. The discretizations must be very carefully constructed to achieve discrete stability using the classical FE method. Alternatively, conditionally stable FE methods can be applied. The latter leads to conditionally stable discretizations requiring arduous a priori analyses to ascertain stabilization parameters. However, the AVS- shown to yield superior results for convex domains.

Time Slice Approach
As an alternative to the space-time discretization of the full space-time domain Ω T , in this section we introduce a time slice approach for the AVS-FE method. While the space-time approach introduced in the preceding section allows straightforward implementation of the AVS-FE method and its "built-in" error indicator can drive mesh adaptive refinements, the large number of degrees of freedom quickly makes the method intractable. In an effort to reduce the computational 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 can be transferred to the neighboring slice as an initial condition. Hence, we can perform mesh adaptive refinements on each slice to ensure the complete resolution of any interior or boundary layer (i.e., physical features) before proceeding to the next. This is of particular interest in applications in which physical parameters are time dependent leading to widely different solution features as time progresses. Hence, it gives an even greater flexibility in the choice of mesh parameters than the "global" space-time approach. In Figure 1 an arbitrary domain Ω T is shown and is partitioned uniformly into two space-time slices Ω ∩ (0, T slice ) and Ω ∩ (0, T ). In Section 4 we compare two possible approaches employing 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 of norms of the approximation errors u − u h and q − q h . We present estimates in terms of the energy norm and classical Sobolev norms on the trial space U(Ω T ). The a priori error estimates of the AVS-FE method for stationary convectiondiffusion problems are established in [32] and form the basis of our results here. To keep this presentation brief, we do not present extensive proofs but we present outlines and references to the appropriate literature. We assume here that both q h and u h are discretized by C 0 (Ω) polynomials and note that approximations using other bases for q h can be 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. its corresponding AVS-FE approximation. Then: where h is the maximum element diameter, µ = min (p u + 1, r), p u the minimum polynomial degree of approximation of u h in the mesh, and r the minimum regularity of the solution components (p, r) of the distributional PDE underlying the 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][35][36] used in the FE approximation of the optimal test functions (9), we establish the bound (39).
The energy norm cannot be computed exactly due to the supremum in its definition (8) and must be computed approximately by computing the approximate error representation function and Proposition 2.1. To establish bounds on the approximation errors in terms of classical Sobolev norms on the trial space, we first note that since the energy norm is an equivalent norm on U(Ω T ), we have that: Due to this norm equivalence, we can readily establish the following a priori estimates: Lemma 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 ) its corresponding AVS-FE approximation (14). Then: where h is the maximum element diameter, µ 1 = min (p u +1, r u ), p u the minimum polynomial degree of approximation of u h in the mesh, µ 2 = min (p q + 1, r q ), p q the minimum polynomial degree of approximation of q h , and r u , r q the regularities 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 (·, ·) U(Ω T ) due to the norm equivalence on U(Ω T ), the best approximation property of the AVS-FE method (in terms of the energy norm), and the convergence properties of piecewise polynomial interpolants [34,35], we establish the last bound in (41). The bounds on individual solution components in (41) are direct consequences of this last bound. The bound on the L 2 (Ω) norm of u − u h in (41) is a direct consequence of an application of the the Aubin-Nitche lift [37,38].
Remark 3.2 Establishing the optimal bound on u − u h L 2 (Ω T ) can be done using well established techniques found in FE literature. Unfortunately, a similar, optimal, bound cannot be established for q − q h L 2 (Ω) as the required duality result for the flux variable can only be established for specially designed mesh geometries (see [39]). However, we do note that extensive numerical experimentation indicate optimal convergence behavior of the flux in the case of convex computational domains. For other approximation spaces, such as Raviart-Thomas, optimal bounds can be established using the techniques found in the text of Brezzi and Fortin [33].

Numerical Verifications
To conduct numerical verifications, we consider the following simplified form of our model scalar-valued convection diffusion problem (1) with Dirichlet boundary conditions: where the coefficient ε ∈ L ∞ (Ω) is a scalar-valued isotropic diffusion coefficient. First, we verify the convergence properties of the AVS-FE method for both time discretization schemes in Section 4.1. The use of the time slice approach is also considered here, and we present verifications of two distinct refinement strategies for this approach.
In Section 4.2, we consider a hyperbolic problem of purely convective transport. Last, in Section 4.3 we present verifications 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 error indicator (19) as well as the marking strategy of [40] using the approximate energy error computed using (17). Also note that in all cases where we report the number of degrees of freedom, we include the degrees of freedom for the error representation function in the saddle point systems (20) and (33). The polynomial degree of approximation used for this error representation function is always identical to the degree of the trial space.

Convergence Studies
To numerically verify the predicted rates of convergence from the a priori bounds above we perform numerical convergence studies. This is accomplished by considering a well-known example of transient convection-diffusion, the Eriksson-Johnson problem [41]: The exact solution of this problem is: where l = 2, and: The problem domain Ω T = (−1, 0) × (−0.5, 0.5) × (0, 1). For these studies we consider the moderately convection dominated case of (43) with ε = 0.1.
In Figure 2 the convergence plots for linear and quadratic polynomial degrees for the space-time approach are shown. In Figure 2, we plot error norms versus the number of degrees of freedom N, which increases at O(h −2 ), i.e., the h−convergence rates of the FE approximations can be extracted from these by a simple adjustment. For the Thus, the corresponding rates of convergence are as predicted by the a priori error estimates of Section 3.2.2 for u − u h L 2 (Ω) and the energy norm in Figure 2. The observed rates for q − q h L 2 (Ω) are as expected for the continuous polynomial approximations we use for linear polynomials approximations and is one half order lower than expected for the quadratic case. We expect the optimal rate to be recovered upon further mesh refinements but the computational cost is too large for our available computational resources 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 of the method at the final time T = 1s with time step of τ = 10 −3 . The observed rates of convergence in Figure 3 are

Time Slices
In this section we present results for the time slice approach introduced in 3.2. In particular, we consider two different approaches to mesh adaptive refinements for the time slice approach. We again consider the Eriksson-Johnson problem of the preceding section, now with ε = 0.05. The time domain is partitioned into two slices from 0 to 0.5s and from 0.5s to 1.0s, the polynomial degree of approximation is chosen to be one for both variables in both approaches.
First, the space-time AVS-FE method is implemented on the bottom slice and we perform adaptive mesh refinements based on the built in error indicator (19) and we employ the marking strategy of [40] using the energy error from the slice to ascertain which elements to refine. We perform a total of 13 steps of adaptation before the solution is used as in an initial condition for the final time slice to which we perform the same number of mesh refinements. The second approach does not perform mesh refinements on each slice individually but rather computes the solution on the first slice and then moves directly to the second slice. Then, we again employ the marking strategy of [40] based on      both approaches.

Pure Convection
In this verification, we consider the problem of purely convective transport, i.e., ε = 0 with a time-dependent source term. To present a numerical verification, we use the method described in (33), the generalized-α method. We consider the spatial domain to be the unit square Ω = (0, 1) 2 ⊂ R 2 , 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: and the initial solution is u 0 = 0. We choose the time-step τ = 5 × 10 −3 and the final time T f inal = 1. We set the parameter M = 500 in (46) and ρ ∞ = 0.9 for the generalized-α 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 adaptive refinement strategy. The initial mesh is coarse and is such that no elements in the mesh align with the expected interior layers of (46) and is shown in Figure 5(a). In Figure 6(b), the converged solution at the final time T = 1.0s is presented and 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 mesh refinements focused near the circular interior layer.

Shock Problem
As a final numerical verification, we present a consideration of (42) in which the solution behaves as two shocks traveling through the space-time domain while rotating about the origin. Furthermore, the choices we make for the   problem parameters are such that the interface of the shock is skewed and rotates in the space-time domain as t → T f inal .
Thus, we have the following choices: . 35s For this particular problem, we present the time slice approach in which we perform mesh adaptations between each slice. Experience has shown that the slice containing the initial condition is critical to the proper resolution of the space-time process. Thus, we consider the case of two space-time slices, the first from 0s to 0.2s and the final from 0.2s to 1.35s. In Figures 7 and 8

Conclusions
The AVS-FE method is a Petrov-Galerkin method which uses classical continuous FE trial basis functions, while the test space consist of functions that are discontinuous across element edges. This broken topology in the test space allows us to employ the DPG philosophy and introduce an equivalent saddle point problem which we implement using high level FE solvers. We have introduced two distinct approaches to transient problems using the AVS-FE method.  First, we take a space-time approach in which the entire space-time domain in discretized using finite elements, and second, 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 a single global solve governing the FE approximation. As the AVS-FE approximations computed from the saddle point system (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-FE method in Section 3.2.2 as for the classical Galerkin FE method due to the best approximation property of the AVS-FE method. In an effort to reduce the computational cost of solving the global system we have consider a time slice approach in which the space-time domain is partitioned into finite sized space-time slices on which we employ the AVS-FE method. The advantage here is that the size of the global system is reduced and we are able to employ mesh adaptive 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 by the generalized-α method leading to highly efficient stable FE computations. We show that the AVS-FE method uses a corresponding norm as a function of the time-step. Another distinguishing feature of this method is that due to the influence of the initial data on the accuracy of the solution, we find a stable approximation for ∂ u ∂t at the initial time. Accordingly, at each time step, one requires to solve a system with a smaller number of degrees of freedom in comparison with the space-time approach.
Numerical verifications for several cases of the transient convection-diffusion IBVP show that both methods exhibit optimal asymptotic convergence behavior as well as similar norms of the numerical approximation error. For degrees of approximation above 2, the space-time approach becomes more accurate as it is not limited to the second-order accuracy of the generalized-α method. However, we do not advocate one method over the other but we point out these differences for potential users as their available computational resources will likely dictate which approach to use. For both cases, we present additional numerical verifiactions highlighting the adaptive mesh refinement capabilities. In future efforts, we expect to pursue alternative error estimators and indicators as well as the AVS-FE approximation of challenging transient physical phenomena.