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

Computational and Applied Mathematics (2023) 42:43 https://doi.org/10.1007/s40314-022-02174-w Spectral variational multi-scale method for parabolic problems: application to 1D transient advection-diffusion equations Tomás Chacón Rebollo1·Soledad Fernández-García1·David Moreno-Lopez1· Isabel Sánchez Muñoz2 Received: 13 July 2022 / Revised: 14 December 2022 / Accepted: 21 December 2022 / Published online: 13 January 2023 © The Author(s) 2023 Abstract 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 advectiondiffusion 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. Keywords Variational multi-scale ·Parabolic problems ·Transient advection-diffusion · Stabilized method ·Spectral approximation Communicated by Forrest Carpenter. BSoledad Fernández-García [email protected] Tomás Chacón Rebollo [email protected] David Moreno-Lopez [email protected] Isabel Sánchez Muñoz [email protected] 1Dpto. EDAN & IMUS, University of Seville, Campus de Reina Mercedes, 41012 Sevilla, Spain 2Dpto. Matemática Aplicada I, University of Seville, Ctra de Utrera s/n, 41013 Sevilla, Spain 123 43 Page 2 of 27 T.C. Rebollo et al. Mathematics Subject Classification 65M12 ·65M15 ·65M60 ·65M70 1 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 cf. Hughes 1995; Hughes and Stewart 1995; Hughes et al. 1998)). The VMS formulation is based upon the formulation of the Galerkin method as two variational problems, one satisfied by the resolved scales and another satisfied by the subgrid 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 elementwise 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 Codina (2000). Within these methods, the effect of the sub-grid scales is modelled by means of a dissipative interaction of operator terms acting on the resolved scales. The VMS methods have been successfully applied to many flow problems, and in particular to Large Eddy Simulation (LES) models of turbulent flows (cf. Hughes et al. 2000; John 2006; Chacón Rebollo and Lewandowski 2014). The application of VMS method to evolution PDEs dates back to the 1990s, when the results from Hughes (1995) were extended to nonsymmetric linear evolution operators, see (Hughes and Stewart 1995). The papers (Harari 2004; Harari and Hauke 2007) deal with the spurious oscillations generated in the Galerkin method for parabolic problems due to very small time steps. The series of articles (Hauke and Doweidar 2005a,b,2006) deal with transient Galerkin and SUPG methods, transient subgrid scale (SGS) stabilized methods and transient subgrid scale/gradient subgrid scale (SGS/GSGS), making a Fourier analysis for the one-dimensional advection-diffusion-reaction equation. A stabilised finite element method for the transient Navier–Stokes equations based on the decomposition of the unknowns into resolvable and subgrid scales is considered in Codina (2002); Codina et al. (2007). Further, (Asensio et al. 2007) compares the Rothe method with the so-called Method of Lines, which consists on first, discretise in space by means of a stabilized finite element method, and then use a finite difference scheme to approximate the solution. More recently, (Chacón Rebollo and Dia 2015) introduced the use of spectral techniques to model the sub-grid scales for 1D steady advection-diffusion equations. The basic observation is that the eigenpairs of the advection-diffusion operator may be calculated analytically on each grid element. A feasible VMS-spectral discretization is then built by truncation of this spectral expansion to a finite number of modes. An enhanced accuracy with respect to preceding VMS methods is achieved. In Chacón Rebollo et al. (2020), the spectral VMS method is extended to 2D steady advection-diffusion problems. It is cast for low-order elements as a standard VMS method with specific stabilised coefficients, that are anisotropic in the sense that they depend on two grid Péclet numbers. To reduce the computing time, the stabilised coefficients are precomputed at the nodes of a grid in an offline step, and then interpolated by a fast procedure in the online computation. The present paper deals with the building of the spectral VMS numerical approximation to evolution advection-diffusion equations. We construct an abstract spectral VMS discretisation of parabolic equations, that is particularised to 1D advection-diffusion equations. The 123 Spectral variational multi-scale method... Page 3 of 27 43 sub-grid components are exactly calculated in terms of spectral expansions when the driving 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 solution provided by the spectral VMS method coincides with the exact solution of the implicit Euler semi-discretisation at the Lagrange interpolation nodes. We also build a feasible method to solve the evolutive advection-diffusion problem by means of an offline/online strategy that pre-computes the action of the sub-grid scales on the resolved scales. This allows to dramatically reduce the computing times required by the method. We further perform some numerical tests for strongly advection dominated flows. The spectral VMS method is found to satisfy the discrete maximum principle, even for very small time steps. A remarkable increase of accuracy with respect to several stabilised methods is achieved. The outline of the paper is as follows. In Sect. 2, we describe the abstract spectral VMS discretisation to linear parabolic problems, which is applied to transient advection-diffusion problems in Sect. 3. A feasible method is built in Sect. 4, based upon an offline/online strategy. We present in Sect. 5our numerical results, and address some conclusions in Sect. 6. 2 Spectral VMS method In this section, we build the spectral VMS discretisation to abstract linear parabolic equation. Let a bounded domain in Rdand T>0 a final time. Let us consider two separable Hilbert spaces on ,Xand H,sothatX⊂Hwith dense and continuous embedding. We denote (·,·)the scalar product in X;Xand Hare the dual topological spaces of Xand H, respectively, and ·,· is the duality pairing between Xand X. We identify Hwith its topological dual Hso that X⊂H≡H⊂X.DenotebyL(X)the space of bilinear bounded forms on Xand consider b∈L1(0,T;L(X)) uniformly bounded and X-elliptic with respect to t∈(0,T). Given the data f∈L2(0,T;X)and u0∈H, we consider the following variational parabolic problem: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Find u∈L2((0,T);X)∩C0([0,T];H)such that, d dt (u(t), v) +b(t;u(t), v) =f(t), v∀v∈X,in D(0,T), u(0)=u0in H. (1) It is well known that this problem is well posed and, in particular, admits a unique solution (Dautray and Lions 1992). To discretize this problem, we proceed through the so-called Horizontal Method of Lines (Asensio et al. 2007;Bernardietal.2004; Harari and Hauke 2007). 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=t0<t1< ... < tN=T}, with time-step 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 u0=u0, ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Find un+1∈Xsuch that, un+1−un t,v+bn+1(un+1,v) =tfn+1,v∀v∈X, ∀n=0,1,...,N−1, (2) 123 43 Page 4 of 27 T.C. Rebollo et al. where bn+1and fn+1are some approximations of b(t;·,·)and f(t), respectively, at t=tn+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 , {Th}h>0, formed by simplicial elements, where the parameter hdenotes the largest diameter of the elements of the triangulation Th. The VMS method is based on the decomposition, X=Xh⊕˜ X, where Xhis a continuous finite element sub-space of Xconstructed on the grid Th,and ˜ X is a complementary, infinite-dimensional, sub-space of X.Notice that this is a multi-scale decomposition of the space X,beingXhthe large or resolved scale space and ˜ Xthe small or sub-grid scale space. This decomposition defines two projection operators Ph:X→ Xh and ˜ P:X→ ˜ X,by Ph(v) =vh,˜ P(v) =˜v, ∀v∈X,(3) where vhand ˜vare the unique elements belonging to Xhand ˜ X, respectively, such that v=vh+˜v. Hence, one can decompose the solution of problem (2)as un+1=un+1 h+˜un+1, where un+1 h=Ph(un+1)and ˜un+1=˜ P(un+1)satisfy the coupled problem, ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ un+1 h−un h t,v h+˜un+1−˜un t,v h+bn+1(un+1 h,v h)+bn+1(˜un+1,v h)=fn+1,v h(4.1) un+1 h−un h t,˜v+˜un+1−˜un t,˜v+bn+1(un+1 h,˜v) +bn+1(˜un+1,˜v) =fn+1,˜v(4.2) ∀vh∈Xh,∀˜v∈˜ X, for all n=0,1,...,N−1. The small scales component ˜un+1thus satisfies, (˜un+1,˜v) +tb n+1(˜un+1,˜v) =Rn+1(un+1 h), ˜v(4) where Rn+1(un+1 h), ˜vis the residual of the large scales component, defined as, Rn+1(un+1 h), ˜v:=(un h+˜un,˜v) +tfn+1,˜v−(un+1 h,˜v) −tb n+1(un+1 h,˜v), ∀˜v∈˜ X. In condensed notation, this may be written as, ˜un+1=n+1(Rn+1(un+1 h)), (5) where n+1:˜ X→˜ X g→ n+1(g)=˜ G is the static condensation operator on ˜ Xdefined as, (˜ G,˜v) +tb n+1(˜ G,˜v) =g,˜v∀˜v∈˜ X,for any g∈˜ X. Inserting expression (5) in the large scales equation (4.1), leads to the condensed VMS formulation of problem (2): ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Find un+1 h∈Xhsuch that (un+1 h,v h)+tb n+1(un+1 h,v h)+(n+1(Rn+1(un+1 h)), vh)+tb n+1(n+1(Rn+1(un+1 h)), vh) =tfn+1,v h+(un h+n(Rn(un h)), vh) ∀vh∈Xh,∀n=0,1,...,N−1, (6) 123 Spectral variational multi-scale method... Page 5 of 27 43 with u0 h=Ph(u0). This problem is an augmented Galerkin formulation, where the additional terms represent the effect of the small scales component of the solution (˜un+1)on the large scales component (un+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 ˜ Xis approximated by the “bubble” sub-spaces, ˜ X≃˜ Xh= K∈Th ˜ XK,with ˜ XK={˜v∈˜ X,such that supp(˜v) ⊂K}.(7) Hence, we approximate ˜un+1≃˜un+1 h= K∈Th˜un+1 K,with ˜un+1 K∈˜ XK,∀n=0,1,...,N−1.(8) Then, problem (4) is approximated by the following family of decoupled problems, (˜un+1 K,˜vK)+tb n+1(˜un+1 K,˜vK)=Rn+1(un+1 h), ˜vK,∀˜vK∈˜ XK,∀K∈Th.(9) Let Ln+1:X→ Xbe the operator defined by Ln+1w, v=bn+1(w, v), ∀w, v ∈X,(10) and let Ln+1 Kbe the restriction of this operator to ˜ XK. Let us also consider the weighted L2 space, L2 p(K)={w:K→Rmeasurable such that p|w|2∈L1(K)}, where pis some measurable real function defined on K,which is positive a.e. on K.Thisis a Hilbert space endowed with the inner product (w, v)p=K p(x)w(x)v(x)dx. We denote by ·pthe norm on L2 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 {˜zn,K j}j∈Non ˜ XKformed by eigenfunctions of the operator Ln K, which is an orthonormal system in L2 pn,K(K)for some weight function pn,K∈C1(¯ K). Then, ˜un K=∞  j=1 βn,K jrn,K j˜zn,K j,∀n=1,...,N,(11) where βn,K j=(n,K j)−1, with n,K j=1+tλn,K jbeing λn,K jthe eigenvalue of Ln K associated to ˜zn,K j, and rn,K j=Rn(un h), pn,K˜zn,K j. This is a rather straightforward application of Theorem 1 in Chacón Rebollo and Dia (2015), that we do not detail for brevity. 123 43 Page 6 of 27 T.C. Rebollo et al. Once the eigenpairs (˜zn+1,K j,λ n+1,K j)are known, the previous procedure allows us to directly compute un+1 hfrom problem (6), approximating the sub-grid component ˜un+1by expressions (8)and(11). This gives the spectral VMS method to fully discretize problem (1). Namely, ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Find un+1 h∈Xhsuch that (un+1 h,v h)+tb n+1(un+1 h,v h)+(˜un+1 h,v h)+tb n+1(˜un+1 h,v h) =tfn+1,v h+(un h,v h)+(˜un h,v h) ∀vh∈Xh,∀n=0,1,...,N−1, (12) where, ˜un+1 h= K∈Th ∞  j=1 βn+1,K jRn+1 h(un+1 h), pn+1,K˜zn+1,K j˜zn+1,K j,∀n=0,...,N−1,(13) with Rn+1 h(un+1 h), ˜v:=(un h+˜un h,˜v) +tfn+1,˜v−(un+1 h,˜v) −tb n+1(un+1 h,˜v), ∀˜v∈˜ X. u0 h=Ph(u0)and ˜u0 h∈˜ Xhsome approximation of ˜u0. 3 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, ⎧ ⎨ ⎩ ∂tu+a·∇u−μu=fin ×(0,T), u=0on∂ ×(0,T), u(0)=u0on , (14) where a∈L∞(0,T;W1,∞())dis the advection velocity field, μ>0 is the diffusion coefficient, f∈L2((0,T);L2()) is the source term and u0∈L2() 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, ⎧ ⎪ ⎨ ⎪ ⎩ Find u∈L2((0,T);H1 0()) ∩C0([0,T];L2()) such that, d dt (u(t), v) +(a·∇u(t), v) +μ(∇u(t), ∇v) =f(t), v∀v∈H1 0(), u(0)=u0. (15) Problem (15) admits the abstract formulation (1) with H=L2(),X=H1 0() and b(w, v) =(a·∇w, v) +μ(∇w, ∇v), ∀w, v ∈H1 0(). 123 Spectral variational multi-scale method... Page 7 of 27 43 In practice, we replace the velocity field aby ah, the piecewise constant function defined a. e. on such that ah=aKon the interior of each element K∈Th. Then, we apply the spectral VMS method to the approximated problem, ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Find Un+1∈H1 0() such that Un+1−Un t,v+(an+1 h·∇Un+1,v)+μ(∇Un+1,∇v) =fn+1,v,∀v∈H1 0(), ∀n=0,1,...,N−1, (16) with u0=u0. In this case, Lnw=an h·∇w−μw is the advection–diffusion operator. Proposition 1 in Chacón Rebollo and Dia (2015) proved that the eigenpairs (˜wn,K j,λ n,K j)of operator Ln K can be obtained from the eigenpairs (˜ WK j,σK j)of the Laplace operator in H1 0(K),inthe following way: ˜wn,K j=ψn,K˜ WK j,ψ n,K(x)=exp 1 2μan K·x λn,K j=μσK j+|an K|2 4μ2,∀j∈N. (17) Moreover, for the weight function pn,K(x)=(ψn,K)−2=exp −1 μaK·x(18) the sequence ˜zn,K j=˜wn,K j ˜wn,K jpn,K ,∀j∈N,(19) is a complete and orthonormal system in L2 pn,K(K)(see Theorem 2 in Chacón Rebollo and Dia (2015)). Then, Theorem 2.1 holds and it is possible to apply the method (12) to problem (16). 3.1 One dimensional 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 K∈Th are closed intervals, K=[a,b]. The eigenpairs (˜ WK j,σK j)are solutions of the problem −∂xx ˜ WK=σK˜ WKin K, ˜ WK(a)=˜ WK(b)=0. Solutions of this problem are ˜ WK j=sin σK j(x−a),σ K j=jπ hK2 ,with hK=b−a,for any j∈N. As the function pn,Kdefined in (18) is unique up to a constant factor, to express the eigenpairs in terms of non-dimensional parameters, we replace pn,Kby (we still denote it in the same 123 43 Page 8 of 27 T.C. Rebollo et al. way), pn,K(x)=exp −2Pn,K x−a hK,(20) where Pn,K=|an K|hK 2μis the element Péclet number. Then, from expressions (17)and(19), ˜zn,K j=2 hK exp Pn,K x−a hKsin jπx−a hK,λ K j=μjπ hK2 +|an K|2 4μ.(21) It follows βn,K j=1 1+SK(P2 n,K+π2j2)for any j∈N,(22) where SK=tμ h2 K 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 Xh.Letωi=(αi−1,α i), i=1,...,I. Setting ˜ Xi=H1 0(ωi), it holds, H1 0() =Xh⊕˜ X,with ˜ X= I  i=1˜ Xi. Observe that this decomposition generalises (7) with ˜ Xh=˜ X. Moreover, when operator in (10)isLnw=an h·∇w−μw, problem (4) can be exactly decoupled into the family of problems (9). In particular, if the projection operator Phin (3) is the Lagrange interpolate on Xh,thenUn h=Ph(Un),˜ Un h=Un−Un h∈˜ Xand consequently, Un h∈Xhsatisfies method (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 analyze the discretisation error we compare the solution of problem (16)to the solution of the implicit Euler time semi-discretisation of problem (15), ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Find un+1∈H1 0() such that un+1−un t,v+(an+1∂xun+1,v)+μ(∂ xun+1,∂ xv) =fn+1,v∀v∈H1 0(), ∀n=0,1,...,N−1, (23) with u0=u0. 123 Spectral variational multi-scale method... Page 9 of 27 43 We assume that ahrestricted to each Kis extended by continuity to ∂K. Given a sequence b={bn,n=1,...,N}of elements of a normed space Y, let us denote, blp(Y)=t N  n=1bnp Y1/p ,bl∞(Y)=max n=1,...,NbnY. 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 (1−σt)α n+1+βn+1≤(1+τt)α n+γn+1 for some σ≥0,τ≥0. Assume that σt≤1−δfor some δ>0. Then it holds αn≤eρtnα0+1 δ n  l=1 eρ(tn−tl)γl, and n  l=1 βl≤1+τ σ+(σ +τ)eρtn−1tn−1α0+1 δ1+(σ +τ)eρtn−1tn−1 n  l=1 γl, with ρ=(σ +τ)/δ. Let e={en,n=0,1,...,N}⊂H1 0() be the sequence of errors en=un−Un∈H1 0(), where we recall that Unis the solution of the discrete problem (16), and denote δten+1= en+1−en t. It holds the following result. Proposition 3.2 Assume that a∈L∞( ×(0,T))d,f∈L2( ×(0,T)),t≤(1− ε) μ a2 L∞(×(0,T)) for some ε∈(0,1)and ahL∞(×(0,T)) ≤DaL∞(×(0,T)) for some constant D >0. Then, δtel2(L2()) +μel∞(H1 0()) ≤Cah−al2(L∞()),(24) for some constant C >0independent of h, t and μ. Proof Let us subtract (16) from (23) with v=vh∈Xh. This yields en+1−en t,v h+(an+1 h∂xen+1,v h)+μ(∂ xen+1,∂ xvh)=((an+1 h−an+1)∂xun+1,v h). Setting vh=δten+1, and using the identity 2(b,b−a)=b2 L2() −a2 L2() +b−a2 L2() for any a,b∈L2()dyields tδten+12 L2() +t(an+1 h∂xen+1,δ ten+1)+μ 2∂xen+12 L2() −∂xen2 L2() ≤t((an+1 h−an+1)∂xun+1,δ ten+1). (25) It holds |(an+1 h∂xen+1,δ ten+1)|≤an+1 hL∞() ∂xen+1L2() δten+1L2() 123 43 Page 16 of 27 T.C. Rebollo et al. Fig. 6 Test 1. l∞(L2)and l2(H1)errors for the spectral VMS solution of problem (33) 5.2 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(considering that the length of is 1), 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 Fig. 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 ctimes-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 u0is exact, then u1 hshould coincide with the exact solution at grid nodes. This can already be observed in Fig. 7(a). We also test this result with different discretisation parameters. We actually set t=10−5and h=0.02 that corresponds to P=10 and S=0.025. The solution in the first time-step is represented in Fig. 8(a)and a zoom around x=0.7indepictedin(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 CFL < CFLbound =P/(3(1−P)) (see Harari and Hauke 2007). For that, we consider the same problem as in this section but with a=20, h=0.01 and the time-step tis chosen such that CFL/CFLbound =1/2. We obtain the results shown in Fig. 9, where we have represented the first five time-steps. As one can see the spectral solution does not present any oscillation. 123 Spectral variational multi-scale method... Page 17 of 27 43 Fig. 7 Solution of problem (14)fora=1000,μ =1,f=0andu0givenby(34) with t=10−3and h=0.02 (P=10, S=2.5). The spectral VMS solution is compared to the exact solution and the Galerkin solution. The results for time-steps numbers 1 to 4, 5 to 7 and 8 to 9 are respectively represented in figures (a),(b)and (c) Fig. 8 Solution of problem (14)fora=1000,μ =1,f=0andu0givenby(34) with t=10−5and 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 123 43 Page 18 of 27 T.C. Rebollo et al. Fig. 9 Solution of problem (14)fora=20,μ=1,f=0andu0givenby(34) with h=0.01 and tsuch that CFL/CFLbound =1/2(P=0.1, S=0.0926). Red lines represent Galerkin solution and cyan lines represent spectral VMS solution in each step-time 5.3 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, (M+tR n+ta 2τMs)un+1=Mun, where Mand Rnare, respectively, mass and stiffness matrices, while Msis a tridiagonal matrix defined by (Ms)i,i=2 h,(Ms)i+1,i=(Ms)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 (Christie et al. 1976; John and Novo 2011), τ1D=μ |a|2(Pcoth(P)−1). (35) 2. The stabilisation coefficient based on orthogonal sub-scales proposed by Codina in Codina (2000), τC=4μ h22+2|a| h2−1/2 .(36) 3. The stabilisation coefficient based on L2proposed by Hauke et al. (2008), τH=min h √3|a|,h2 24.24μ,t.(37) 123 Spectral variational multi-scale method... Page 19 of 27 43 Fig. 10 Comparison of different stabilised methods to solve problem (14)whenP=3, S=25 with t=10−2and h=0.02. Solutions in the three first time-steps Table 1 l∞(L2)and l2(H1) errors for the solutions represented in Fig. 10 Method l∞(L2)l2(H1) P=3 Galerkin 1.1784e-02 4.7505e-02 S=25 Spectral 8.7889e-06 5.4716e-05 Codina 3.2285e-03 1.4329e-02 1D 1.3805e-03 1.3446e-03 Hauke 2.1713e-03 1.1124e-02 Franca 9.9020e-03 5.0380e-02 4. The stabilisation coefficient separating the diffusion-dominated from the convectiondominated regimes proposed by Franca in Chacón Rebollo et al. (2015), τF=h |a|min{P,˜ P},(38) where ˜ P>0 is a threshold separating the diffusion dominated (P≤˜ P) to the advection dominated (Pe >˜ P)regimes. In Figs. 10,11 and 12, we show the solutions of each method for different values of Pand S, always for advection-dominated regime P>1. We also display the errors in l∞(L2)and l2(H1)norms for the solutions of these problems in Tables 1,2and 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. Next, we consider the same tests performed in Sect. 5.2, but applying the feasible spectral VMS method. Firstly, we check the behaviour of the feasible spectral VMS method (29)forverylarge Péclet numbers. In Fig. 13 we represent the solution of same problem as in Fig. 7obtained with this method. We show solutions in time-steps 1 to 4 in (a), times-steps 5 to 7 in (b) 123 43 Page 20 of 27 T.C. Rebollo et al. Fig. 11 Comparison of different stabilised methods to solve problem (14)whenP=1andS=5 with t=10−3and h=10−2. Solutions in the three first time-steps. Right: zoom around x=0.7 Table 2 l∞(L2)and l2(H1) errors for the solutions represented in Fig. 11 Method l∞(L2)l2(H1) P=1 Galerkin 9.6551e-03 7.7424e-02 S=5 Spectral 7.2887e-05 5.2396e-04 Codina 1.3580e-02 6.4992e-02 1D 3.7524e-03 5.3902e-03 Hauke 4.2353e-03 3.3330e-02 Franca 4.4200e-02 3.1419e-01 Fig. 12 Comparison of different stabilised methods to solve problem (14)whenP=3.5andS=100 with t=10−2and h=10−2. Solutions in the three first time-steps. Right: zoom around x=0.99 Table 3 l∞(L2)and l2(H1) errors for the solutions represented in Fig. 12 Method l∞(L2)l2(H1) P=3.5 Galerkin 4.5006e−03 3.3305e−02 S=100 Spectral 1.6381e−06 2.0138e−05 Codina 8.752e−04 8.4455e−03 1D 3.3968e−04 4.5556e−04 Hauke 5.6656e−04 5.6930e−03 Franca 3.0238e−03 3.0336e−02 123 Spectral variational multi-scale method... Page 21 of 27 43 Fig. 13 Solution of problem (14)fora=1000,μ =1,f=0andu0givenby(34) with t=10−3and h=0.02 (P=10, S=2.5). The feasible spectral VMS is compared with different stabilised methods. The results for time-steps numbers 1 to 4, 5 to 7 and 8 to 9 are respectively represented in figures (a),(b)and (c) 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, Fig. 14 is the analogous to Fig. 8, but comparing the feasible spectral VMS with different stabilised methods. Although Hauke’s solution is closer to the exact solution than 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 CFL < CFLbound .InFig.15, we can see the first five time-steps solutions obtained with each method using a time-step that verifies CFL/CFLbound =1/2. In the following subsection we compare the computing costs of the Spectral method with and without offline/online strategy with the other considered methods. 5.4 Computing cost We are going compare the computing time required by the methods introduced in expressions (35)–(38) with the spectral method with and without offline/online strategy. We have 123 43 Page 22 of 27 T.C. Rebollo et al. Fig. 14 First time-step solution of problem (14)fora=1000,μ =1,f=0andu0given by (34) with t=10−5and h=0.02 (P=10, S=0.025). The feasible spectral VMS is compared with different stabilised methods in the whole domain =(0,1)in (a)and in a zoom around x=0.7in(b) Fig. 15 Solution of problem (14)fora=20,μ=1,f=0andu0givenby(34) with h=0.01 and tsuch that CFL/CFLbound =1/2(P=0.1, S=0.0926). The feasible spectral VMS is compared with different stabilised methods 123 Spectral variational multi-scale method... Page 23 of 27 43 Table 4 Computing time required by each method to solve problem (14)whenP=3, S=25 (a=300,μ=1,h= 0.02,t=0.01) and u0given in (34)andT=1 in increasing times order Method Computing time Galerkin 0.015300 s Spectral offline/online 0.015966 s Hauke 0.028337 s Franca 0.039018 s 1D 0.042189 s Codina 0.034441 s Spectral 0.072913 s implemented all considered methods with the Matlab application R2022b (9.13.0) in a Mac Mini Model M1 2020 computer. First, let us note that the offline stage described in Subsection 4.1, which is performed once for all, has taken 3.15 hours. On the other hand, Table 4shows the computing times (in seconds) required for each method to solve problem (14)whenP=3, S=25 (a=300,μ= 1,h=0.02,t=0.01), u0given in (34)andT=1. In the case of the offline/online strategy, the time corresponds to the online stage. As we can see, Spectral method without offline/online strategy requires approximately two to three times the time required by the other considered methods. However, this disadvantage is solved when we use the offline-online strategy. 6 Conclusions In this paper we have extended to parabolic problems the spectral VMS method developed in Chacón Rebollo and Dia (2015) for elliptic problems. We have constructed a feasible method to solve the evolutive advection-diffusion problem by means of an offline/online strategy that pre-computes the effect of the sub-grid scales on the resolved scales. We have proved that when Lagrange finite element discretisations in space are used, the solution obtained by the fully spectral VMS method (12) coincides with the exact solution of the implicit Euler semi-discretisation of the advection-diffusion problem at the Lagrange interpolation nodes. We have performed some numerical tests that have confirmed this property for very large Péclet numbers and very small time steps, by the fully spectral VMS method. Also some additional tests show an improved accuracy with respect to several stabilised methods for the feasible spectral VMS method (29), with moderate gains of computing times. The methodology introduced here may be extended to multi-dimensional advectiondiffusion equations, by parameterising the sub-grid scales in an offline step. This research is at present in progress. Acknowledgements The research of D. Moreno has been funded by Programa Operativo FEDER Andalucía 2014-2020 grant US-1254587. The research of T. Chacón, S. Fernández and I. Sánchez has been partially funded by Plan Estatal 2021-2023- Proyectos Investigación Orientada Grant PID2021-123153OB-C21. Funding Funding for open access publishing: Universidad de Sevilla/CBUA. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, 123 43 Page 24 of 27 T.C. Rebollo et al. and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Appendix: Matrix formulation of the scheme Problem (29) is equivalent to a linear system with a particular structure that we describe next. If {ϕm}L+1 m=1is a basis of the space Xh, the solution un+1 his obtained as un+1 h= L+1  m=1 un+1 mϕm, where un+1=(un+1 1,...,un+1 L+1)t∈RL+1is the unknown vector. Taking vh=ϕl, with l=1...L, each term in (29) can be written in the following way: (un+1 h,ϕ l)= L  m=1 (ϕm,ϕ l)un+1 m=Mun+1l bn+1(un+1 h,ϕ l)= L  m=1 bn+1(ϕm,ϕ l)un+1 m=Rn+1un+1l (fn+1,ϕ l)=Fn+1l where (M)lm =(ϕm,ϕ l), (Rn+1)lm =bn+1(ϕm,ϕ l), (˜un+1 h,ϕ l)= L  m=1 K∈Th ∞  j=1 βn+1,K j(ϕm,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l)un m + K∈Th ∞  j=1 βn+1,K j(˜un h,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l) +t K∈Th ∞  j=1 βn+1,K j(fn+1,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l) − L  m=1 K∈Th ∞  j=1 βn+1,K j(ϕm,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l)un+1 m −t L  m=1 K∈Th ∞  j=1 βn+1,K jbn+1(ϕm,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l)un+1 m =An+1 1un+Gn+1 1+tFn+1 1−An+1 1un+1−tA n+1 2un+1l, with (An+1 1)lm = K∈Th ∞  j=1 βn+1,K j(ϕm,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l), (39) 123 Spectral variational multi-scale method... Page 25 of 27 43 (An+1 2)lm = K∈Th ∞  j=1 βn+1,K jbn+1(ϕm,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l), (Fn+1 1)l= K∈Th ∞  j=1 βn+1,K jfn+1,pn+1,K˜zn+1,K j(˜zn+1,K j,ϕ l), (Gn+1 1)l= K∈Th ∞  j=1 βn+1,K j(˜un h,pn+1,K˜zn+1,K j)(˜zn+1,K j,ϕ l), (40) and bn+1(˜un+1 h,v h)= L  m=1 K∈Th ∞  j=1 βn+1,K j(ϕm,pn+1,K˜zn+1,K j)bn+1(˜zn+1,K j,ϕ l)un m + K∈Th ∞  j=1 βn+1,K j(˜un h,pn+1,K˜zn+1,K j)bn+1(˜zn+1,K j,ϕ l) +t K∈Th ∞  j=1 βn+1,K j(fn+1,pn+1,K˜zn+1,K j)bn+1(˜zn+1,K j,ϕ l) − L  m=1 K∈Th ∞  j=1 βn+1,K j(ϕm,pn+1,K˜zn+1,K j)bn+1(˜zn+1,K j,ϕ l)un+1 m −t L  m=1 K∈Th ∞  j=1 βn+1,K jbn+1(ϕm,pn+1,K˜zn+1,K j)bn+1(˜zn+1,K j,ϕ l)un+1 m =An+1 3un+Gn+1 2+tFn+1 2−An+1 3un+1−tA n+1 4un+1l, where (An+1 3)lm = K∈Th ∞  j=1 βn+1,K j(ϕm,pn+1,K˜zn+1,K j)bn+1(˜zn+1,K j,ϕ l), (41) (An+1 4)lm = K∈Th ∞  j=1 βn+1,K jbn+1(ϕm,pn+1,K˜zn+1,K jbn+1(˜zn+1,K j,ϕ l), (Fn+1 2)l= K∈Th ∞  j=1 βn+1,K jfn+1,pn+1,K˜zn+1,K jbn+1(˜zn+1,K j,ϕ l), (Gn+1 2)l= K∈Th ∞  j=1 βn+1,K j(˜un h,pn+1,K˜zn+1,K j)bn+1(˜zn+1,K j,ϕ l). (42) Taking in account the definition of ˜un hin (30), the second terms Gn+1 1and Gn+1 2can be expressed in the following way: Gn+1 1l= K∈Th ∞  j=1 βn+1,K jβn,K jˆ Rn h(un h), pn,K˜zn,K j(˜zn+1,K j,ϕ l =Bn+1 1un−1+tFn+1 3−Bn+1 1un−tB n+1 2unl, 123