Spectral variational multi-scale method for parabolic problems: application to 1D transient advection-diffusion equations

In this work, we introduce a variational multi-scale (VMS) method for the numerical approximation of parabolic problems, where sub-grid scales are approximated from the eigenpairs of associated elliptic operator. The abstract method is particularized to the one-dimensional advection-diffusion equations, for which the sub-grid components are exactly calculated in terms of a spectral expansion when the advection velocity is approximated by piecewise constant velocities on the grid elements.We prove error estimates that in particular imply that when Lagrange finite element discretisations in space are used, the spectral VMS method coincides with the exact solution of the implicit Euler semi-discretisation of the advection-diffusion problem at the Lagrange interpolation nodes. We also build a feasible method to solve the evolutive advection-diffusion problems by means of an offline/online strategy with reduced computational complexity.We perform some numerical tests in good agreement with the theoretical expectations, that show an improved accuracy with respect to several stabilised methods.


Introduction
The Variational Multi-Scale is a general methodology to deal with the instabilities arising in the Galerkin discretisation of PDEs (Partial Differential Equations) with terms of different derivation orders (see Hughes (cf.[18,19,20])).
The VMS formulation is based upon the formulation of the Galerkin method as two variational problems, one satisfied by the resolved and another satisfied by the sub-grid scales of the solution.To build a feasible VMS method, the sub-grid scales problem is approximately solved by some analytic or computational procedure.In particular, an element-wise diagonalisation of the PDE operator leads to the Adjoint Stabilised Method, as well as to the Orthogonal Sub-Scales (OSS) method, introduced by Codina in [4].Within these methods, the effects of the sub-grid scales is modelled by means of a dissipative interaction of operator terms acting on and •, • is the duality pairing between X and X.We identify H with its topological dual H so that X ⊂ H ≡ H ⊂ X .Denote by L(X) the space of bilinear bounded forms on X and consider b ∈ L 1 (0, T ; L(X)) uniformly bounded and X-elliptic with respect to t ∈ (0, T ).
Given the data f ∈ L 2 (0, T ; X ) and u 0 ∈ H, we consider the following variational parabolic problem: Find u ∈ L 2 ((0, T ); X) ∩ C 0 ([0, T ]; H) such that, It is well known that this problem is well posed and, in particular, admits a unique solution [11].To discretize this problem, we proceed through the so-called Horizontal Method of Lines [1,2,13].First, we discretise in time by the Backward Euler scheme and then we apply a steady spectral VMS method to the elliptic equations appearing at each time step.Consider a uniform partition of the interval [0, T ], {0 = t 0 < t 1 < ... < t N = T }, with timestep size ∆t = T /N .The time discretization of problem (1) by the Backward Euler scheme gives the following family of stationary problems: given the initialization u 0 = u 0 , where b n+1 and f n+1 are some approximations of b(t; •, •) and f (t), respectively, at t = t n+1 .
To discretise in space problem (2), we assume that Ω is polygonal (when d = 2) or polyhedric (when d = 3), and consider a family of conforming and regular triangulations of Ω, {T h } h>0 , formed by simplycial elements, where the parameter h denotes the largest diameter of the elements of the triangulation T h .The VMS method is based on the decomposition, X = X h ⊕ X, where X h is a continuous finite element sub-space of X constructed on the grid T h , and X is a complementary, infinite-dimensional, sub-space of X.Notice that this is a multi-scale decomposition of the space X, being X h the large or resolved scale space and X the small or sub-grid scale space.This decomposition defines two projection operators P h : X → X h and P : X → X, by where v h and ṽ are the unique elements belonging to X h and X, respectively, such that v = v h + ṽ.Hence, one can decompose the solution of problem (2) as where u n+1 h = P h (u n+1 ) and ũn+1 = P (u n+1 ) satisfy the coupled problem, for all n = 0, 1, . . ., N − 1.The small scales component ũn+1 thus satisfies, where R n+1 (u n+1 h ), ṽ is the residual of the large scales component, defined as, R n+1 (u n+1 h ), ṽ := (u n h + ũn , ṽ) + ∆t f n+1 , ṽ − (u n+1 h , ṽ) − ∆t b n+1 (u n+1 h , ṽ), ∀ ṽ ∈ X.In condensed notation, this may be written as, where is the static condensation operator on X defined as, ( G, ṽ) + ∆t b n+1 ( G, ṽ) = g, ṽ ∀ ṽ ∈ X, for any g ∈ X .
Inserting expression (5) in the large scales equation (4.1), leads to the condensed VMS formulation of problem (2): with u 0 h = P h (u 0 ).This problem is an augmented Galerkin formulation, where the additional terms represents the effect of the small scales component of the solution (ũ n+1 ) on the large scales component (u n+1 h ).To build an approximation of the sub-grid scales, we use a spectral decomposition of the operator associated to the variational formulation on each grid element, at each discrete time.To apply this approximation to problem (6), the small scales space X is approximated by the "bubble" sub-spaces,
Hence, we approximate Then, problem (4) is approximated by the following family of decoupled problems, Let L n+1 : X → X be the operator defined by and let L n+1 K be the restriction of this operator to XK .Let us also consider the weighted L 2 space, where p is some measurable real function defined on K, which is positive a.e. on K.This is a Hilbert space endowed with the inner product We denote by • p the norm on L 2 p (K) induced by this inner product.Now, we can state the following result, which allows to compute the small scales on each grid element by means a spectral expansion.
Theorem 2.1.Let us assume that there exists a complete sub-set {z n,K j } j∈N on XK formed by eigenfunctions of the operator L n K , which is an orthonormal system in L 2 p n,K (K) for some weight function p n,K ∈ C 1 ( K).Then, where This is a rather straightforward application of Theorem 1 in [7], that we do not detail for brevity.

Application to transient advection-diffusion problems
In this section, we apply the abstract spectral VMS method introduced in the previous section to transient advection-diffusion equations, that we state with homogeneous boundary conditions, where a ∈ L ∞ (0, T ; W 1,∞ (Ω)) d is the advection velocity field, µ > 0 is the diffusion coefficient, f ∈ L 2 ((0, T ); L 2 (Ω)) is the source term and u 0 ∈ L 2 (Ω) is the initial data.Different boundary conditions may be treated as well, as these also fit into the general spectral VMS method introduced in the previous section.The weak formulation of problem ( 14) reads, Problem (15) admits the abstract formulation (1) with In practice, we replace the velocity field a by a h , the piecewise constant function defined a. e. on Ω such that a h = a K on the interior of each element K ∈ T h .Then, we apply the spectral VMS method to the approximated problem, with u 0 = u 0 .In this case, L n w = a n h • ∇w − µ∆w is the advection-difusion operator.Proposition 1 in [7] proved that the eigenpairs ( wn,K j , λ n,K j ) of operator L n K can be obtained from the eigenpairs ( W K j , σ K j ) of the Laplace operator in H 1 0 (K), in the following way: Moreover, for the weight function is a complete and orthonormal system in L 2 p n,K (K) (see Theorem 2 in [7]).Then, Theorem 2.1 holds and it is possible to apply the method (12) to problem (16).

One dimension problems
The eigenpairs of the Laplace operator can be exactly computed for grid elements with simple geometrical forms, as it is the case of parallelepipeds.In the 1D case, the elements Solutions of this problem are As the function p n,K defined in ( 18) is unique up to a constant factor, to express the eigenpairs in terms of non-dimensional parameters, we replace p n,K by (we still denote it in the same way), where is the element Péclet number.Then, from expressions ( 17) and ( 19), It follows where is a non-dimensional parameter that represents the relative strength of the time derivative and diffusion terms in the discrete equations, at element K.

Error analysis
We afford in this section the error analysis for the solution of the 1D evolutive convectiondiffusion problem by the spectral VMS method (12).Let {α i } I i=0 ∈ Ω be the Lagrange interpolation nodes of space X h .Let Xi .
Observe that this decomposition generalises (7) with Xh = X.Moreover, when operator in (10) can be exactly decoupled into the family of problems (9).In particular, if the projection operator P h in (3) is the Lagrange interpolate on X h , then (12).Notice that thanks to the spectral expansion, the sub-grid scales contribution in method (12), when the advection velocity is element-wise constant, is exactly computed, and then, the discretisation error only is due to the time discretisation and the approximation of the advection velocity a, but not to the space discretisation.
Therefore, to analize the discretisation error we compare the solution of problem ( 16) to the solution of the implicit Euler time semi-discretisation of problem (15), We assume that a h restricted to each K is extended by continuity to ∂K.Given a sequence We shall use the following discrete Gronwall's lemma, whose proof is standard, and so we omit it.Lemma 3.1.Let α n , β n , γ n , n = 1, 2, ... be non-negative real numbers such that for some σ ≥ 0, τ ≥ 0. Assume that σ ∆t ≤ 1 − δ for some δ > 0. Then it holds and , where we recall that U n is the solution of the discrete problem (16), and denote δ t e n+1 = e n+1 − e n ∆t .It holds the following result.
Corollary 3.3.Under the hypotheses of Proposition 3.2, it holds for some constant C > 0. Moreofer, if a is constant, then the solution U n h of the spectral VMS method (12) coincides with the solution u n of the implicit Euler time semi-discretisation (23) at the Lagrange interpolation nodes of space X h .Proof.In one space dimension H 1 (Ω) is continuously injected in L ∞ (Ω).Then estimate (28) follows from estimate (24).
If a is constant obviously ) at the Lagrange interpolation nodes α i , i = 1, . . ., I, then U n h coincides with u n at these nodes.
4 Feasible method: offline/online strategy Building the spectral VMS method using the formulation (12) requires quite large computing times, due to the summation of the spectral expansions that yield the coefficients of the matrices that appear in the algebraic expression of the method.In order to reduce this time, we shall neglect the dependency of method (12) w.r.t.ũn−1 .Then, our current discretization of problem (1) is the following, where ũn+1 h is given by ( 13), but ũn h is defined from an approximated residual: with Rn h (u n h ), ṽ = (u n−1 h , ṽ) + ∆t f n , ṽ − (u n h , ṽ) − ∆t b n (u n h , ṽ), ∀ṽ ∈ X. Neglecting the dependency of method (12) w.r.t.ũn−1 allows to eliminate the recurrence in time of the sub-grid scales.Thanks to this fact, problem (29) is equivalent to a linear system (that we describe in detail in Appendix), whose coefficients only depend on non-dimensional parameters.

Application to 1D transient advection-diffusion problems
In this case the coefficients of the linear system equivalent to problem (29) only depend on two non-dimensional parameters, as we confirm below.
As we can see in Appendix, if {ϕ m } L+1 m=1 is a basis of the space X h associated to a partition {x 1 < x 2 < . . .< x L+1 } of Ω, the solution u n+1 h of (29) can be written as Then, the unknown vector is the solution of the linear system where the matrix and second term are defined in (43) from matrices A n+1 i and B n+1 i given by ( 35)-( 38) and ( 39)-(42).
We focus, for instance, on the coefficients of matrix A n 1 : From expressions ( 20) and ( 21), p n,K and zn,K j depend on the element non-dimensional parameters P n,K and S K and the non-dimensional variable readily proves that these expressions (up to a factor depending on h) can be written as functions of S K and P K .Further, by (22) the coefficients β n,K j also depend on P n,K and S K .Then, for each K ∈ T h the spectral expansion that determines the element contribution to coefficient (A n 1 ) lm , that is, is a function of P n,K and S K , up to a factor depending on h.This also holds for the coefficients of all other matrices that defines the linear system (31), A n i and B n i , as these are built from the basic values (ϕ m , p n,K zn,K j ), (z n,K j , ϕ l ), b n (ϕ m , p n,K zn,K j ) and b n (z n,K j , ϕ l ).We take advantage of this fact to compute these matrices in a fast way, by means of an offline/online computation strategy.

Offline stage
In the offline stage we compute the element contribution to the coefficients of all matrices appearing in system (31) as a function of the two parameters P and S, that take values at the nodes of a uniform grid, between minimum and maximum feasible values of these parameters.That is, In order to set these values, we consider the piecewise affine finite element functions associated to a uniform partition of Ω with step h.In practical applications the advection dominates and P takes values larger than 1.Also, taking usual values of diffusion coefficient and h ∆t, S takes low positive values.Moreover, when we compute the spectral series that determines the coefficients of the system matrices as functions of P and S, we observe that these values are nearly constant as P and S approaches 20.For instance, we can see in figures 1, 2 and 3 how the spectral series for the diagonal coefficient of A 3 matrix tend to a constant value as P or S increase to 20.Therefore, in numerical tests, we will consider a step ∆ = 0.02 and M = 1000 in (32).To do the computations in this stage, in order to avoid computational roundoff problems due to large velocities, we express the eigenfunctions of the advection-diffusion operator given in (19) in terms of the midpoint of the grid elements xl,l+1 That is, we consider h K , for any j ∈ N.
We further truncate the spectral series neglecting all the terms following to the first term that reaches an absolute value less than a prescribed threshold ε.Actually, we have taken ε = 10 −10 .In Figure 4 we represent the number of these summands needed to reach a first term with absolute value smaller than this ε for the series defining the diagonal coefficient of A 3 matrix.As we can see, more terms are needed as P increases and as S decreases to 0.

Online stage
In the online stage, for each grid element K we compute the contribution of this element to the coefficients of all matrices appearing in system (31).Then, we sum up over grid elements, to calculate these coefficients.
For that, we determine P K and S K and find the indices i, j ∈ 1, . . ., M such that (P K , S K ) belongs to [P i , P i+1 ] × [S j , S j+1 ].In other case, if P k < ∆ we set i = 1 and if P K > ∆M we set i = M − 1, and similarly for j in terms of S K .
As we see above, each element contribution is a function of P K and S K that we denote C(P K , S K ) in a generic way.For instance, for matrix A n 1 , Then, we compute C(P k , S K ) by the following second-order interpolation formula: where the α k are the four corners of the cell [P i , P i+1 ] × [S j , S j+1 ], Q = ∆ 2 is its area and the Q k are the areas of the four rectangles in which the cell is split by (P k , S K ) (see Figure 5).

Numerical Tests
In this section, we present the numerical results obtained with the spectral method to solve 1D advection-diffusion problems.Our purpose, on the one hand, is to confirm the theoretical results stated in Corollary 3.3 for the spectral VMS method and, on the other hand, test the accuracy of the spectral VMS and feasible spectral VMS methods for problems with strong advection-dominance, in particular by comparison with several stabilised methods.

Test 1: Accuracy of spectral VMS method for constant advection velocity
To test the property stated in Corollary 3.3, we consider the following advection-diffusion problem: whose exact solution is given by exp(x + (µ − a)t).
We set T = 0.1, a = 1 and µ = 20.We apply the spectral VMS method (12) to solve this problem with time step ∆t = 0.01 and piecewise affine finite element space on a uniform partition of interval (0, 1) with steps h = 0.05/(2 i ) for i = 2, 3, ..., 7. We have truncated the spectral expansions that yield the small scales ũn h to 10 eigenfunctions.The errors in l ∞ (L 2 ) and l 2 (H 1 ) norms computed at grid nodes are represented in Figure 6.We observe that, indeed, the errors quite closely do not depend on the space step h.
Moreover, we have computed the convergence orders in time, obtaining very closely order 1 in l 2 (H 1 ) norm and order 2 in l ∞ (L 2 ) norm, as could be expected.
In the following numerical experiments we consider the 1D problem (14) setting Ω = (0, 1), with constant velocity field a, source term f = 0 and the hat-shaped initial condition We also set X h to be the piecewise affine finite element space constructed on a uniform partition of interval (0, 1) with step size h.

Test 2: Accuracy of spectral VMS method Very large Péclet numbers
In this test we examine the accuracy of spectral VMS method for very high Péclet numbers.To do that, we set a = 1000 and µ = 1, and solve this problem by the spectral VMS method ( 12), truncating to 150 spectral basis functions the series (13) that yield the sub-grid components.
The solution interacts with the boundary condition at x = 1 in times of order 1/a, that is, 10 −3 .
We then set a time-step ∆t = 10 −3 .Moreover, we set h = 0.02 that corresponds to P = 10 and S = 2.5.We present the results obtained in Figure 7, where we represent the Galerkin solution (in red) on the left panels and the spectral solution (in cyan) on the right panels, both with the exact solution (in blue): in (a) the first 4 time-steps, in (b) time-steps from 5 to 7 and in c times-steps 8 and 9.By Corollary 3.3 the discrete solution coincides at the grid nodes with the exact solution of the implicit Euler semi-discretisation, the expected errors at grid nodes are of order ∆t = 10 −3 .We can see that the spectral solution indeed is very close to the exact solution at grid nodes.
As the discrete solution coincides at the grid nodes with the exact solution of the implicit Euler semi-discretisation and u 0 is exact, then u 1 h should coincide with the exact solution at grid nodes.This can already be observed in Figure 7 (a).We also test this result with different discretisation parameters.We actually set ∆t = 10 −5 and h = 0.02 that corresponds to P = 10 and S = 0.025.The solution in the first time-step is represented in Figure 8 (a) and a zoom around x = 0.7 in depicted in (b).Indeed the discrete solutions coincides with the exact one at grid nodes.

Very small time steps
We test here the arising of spurious oscillations due to extra small time-steps.These spurious oscillations occur in the solutions provided by the Galerkin discretisation when CF L < CF L bound = P/(3(1 − P )) (see [13]).For that, we consider the same problem as in this section but with a = 20, h = 0.01 and the time-step ∆t is chosen such that CF L/CF L bound = 1/2.We obtain the results shown in Figure 9, where we have represented the first five time-steps.As one can see the spectral solution does not present any oscillation.

Test 3: Accuracy of the feasible spectral VMS method. Comparison with other stabilised methods
We next proceed to compare the results obtained with the feasible spectral VMS method (29) with those obtained by several stabilised methods.Stabilised methods add specific stabilising terms to the Galerkin discretisation, generating the following matrix scheme, where M and R n are, respectively, mass and stiffness matrices, while M s is a tridiagonal matrix defined by (M s ) i,i = 2 h , (M s ) i+1,i = (M s ) i,i+1 = − 1 h .Each stabilised method is determined by the stabilised coefficient τ .In particular, we consider: 1.The optimal stabilisation coefficient for 1D steady advection-diffusion equation [10,23],  2. The stabilisation coefficient based on orthogonal sub-scales proposed by Codina in [4], 3. The stabilisation coefficient based on L 2 proposed by Hauke et.al. in [17], , ∆t .
4. The stabilisation coefficient separating the diffusion-dominated from the convection-dominated regimes proposed by Franca in [8], where P > 0 is a threshold separating the diffusion dominated (P ≤ P ) to the advection dominated (P e > P ) regimes.
In figures 10, 11 and 12, we show the solutions of each method for different values of P and S, always for advection-dominated regime P > 1.We also display the errors in l ∞ (L 2 ) and l 2 (H 1 ) norms for the solutions of these problems in tables 1, 2 and 3.As it can be observed in the three tables, spectral method reduces the error between 10 and 100 times compared to the stabilised methods, without presenting oscillations.Table 1: l ∞ (L 2 ) and l 2 (H 1 ) errors for the solutions represented in Figure 10.
Next, we consider the same tests performed in Section 5.2, but applying the feasible spectral VMS method.
Firstly, we check the behaviour of the feasible spectral VMS method (29) for very large Péclet numbers.In Figure 13 we represent the solution of same problem as in Figure 7 obtained with this method.We show solutions in time-steps 1 to 4 in (a), times-steps 5 to 7 in (b) and time-steps 8 and 9 in (c).As we can observe, the spectral method is the closest to the reference solution without presenting any spurious oscillations.
Secondly, Figure 14 is the analogous to Figure  the spectral method, we can see on the right figure that this approximation does not satisfy the Maximun Principle.Finally, we illustrate the fact that the feasible spectral VMS method is the only method among those studied that does not have oscillations for small time steps, when CF L < CF L bound .In Figure 15, we can see the first five time-steps solutions obtained with each method using a time-step that verifies CF L/CF L bound = 1/2.Regarding computing times, by means of the offline/online strategy the feasible spectral VMS method requires somewhat larger computing times than the remaining stabilised methods, due to the interpolation step to build the system matrices.

Figure 1 :
Figure 1: Values of the spectral series to compute the diagonal coefficient of matrix A 3 for each pair (P, S).

Figure 4 :
Figure 4: Number of summands needed to reach a first term with absolute value lower than ε = 10 −10 for the series defining the diagonal coefficient of matrix A 3 , in terms of (P, S).

Figure 5 :
Figure 5: Splitting of interpolation cell for online computation of matrices coefficients.

Figure 8 :Figure 9 :
Figure 8: Solution of problem (14) for a = 1000, µ = 1, f = 0 and u 0 given by (34) with ∆t = 10 −5 and h = 0.02 (P = 10, S = 0.025).The spectral VMS solution is compared to the exact solution and the Galerkin solution at first time step.Figures (a) and (b) respectively show these solutions in the whole domain and a zoom around x = 0.7.

Table 2 :
8, but comparing the feasible spectral VMS with different stabilised methods.Although Hauke's solution is closer to the exact solution than Comparison of different stabilised methods to solve problem (14) when P = 1 and S = 5 with ∆t = 10 −3 and h = 10 −2 .Solutions in the three first time-steps.Right: zoom around x = 0.7.l ∞ (L 2 ) and l 2 (H 1 ) errors for the solutions represented in Figure11.