Guaranteed error control bounds for the stabilised

www.oeaw.ac.at

www.ricam.oeaw.ac.at

Guaranteed error control bounds for the stabilised

space-time IgA

approximations to parabolic problems

U. Langer, S. Matculevich, S. Repin

RICAM-Report 2017-43

Guaranteed error control bounds for the stabilised space-time IgA approximations to parabolic problems

Ulrich Langer

^∗

, Svetlana Matculevich

^†

, and Sergey Repin

^‡

December 15, 2017

Abstract

The paper is concerned with space-time IgA approximations of parabolic initial-boundary value problems.

We deduce guaranteed and fully computable error bounds adapted to special features of IgA approximations and investigate their applicability. The derivation method is based on the analysis of respective integral iden- tities and purely functional arguments. Therefore, the estimates do not contain mesh-dependent constants and are valid for any approximation from the admissible (energy) class. In particular, they imply computable error bounds for norms associated with stabilised space–time IgA approximations. The last section of the paper contains a series of numerical examples where approximate solutions are recovered by IgA techniques.

They illustrate reliability and efficiency of the error estimates presented.

1 Introduction

Time dependent systems governed by parabolic partial differential equations (PDEs) are typical models in sci- entific and engineering applications. This triggers their active investigation in modelling, mathematical analysis and numerical solution. By virtue of fast development of parallel computers, treating time as yet another dimen- sion in space in evolutionary equations became quite natural. The space-time approach is not affected by the disadvantages of time-marching schemes. Its various versions can be useful in combination with parallelisation methods, e.g., those discussed in [12, 29].

Investigation of effective adaptive refinement methods is crucial for constructing fast and efficient solvers for PDEs. At the same time, scheme localisation is strongly linked with reliable and quantitatively efficient a posteriori error estimation tools. These tools are intended to identify the areas of a computational domain with relatively high discretisation errors and by that provide a fully automated refinement strategy in order to reach the desired accuracy level for the current approximation. Local refinement tools in IgA such as T-splines, THB-splines, and LR-splines were combined with various a posteriori error estimation techniques, e.g., error estimates using hierarchical bases [9, 40], residual-based [18, 41, 5, 23, 13], and goal-oriented error estimates [39, 8, 24, 25]. Below we use a different (functional) method providing fully guaranteed error estimates in various weighted norms equivalent to the global energy norm. These estimates include only global constants (independent of the mesh characteristich) and are valid for any approximation from the admissible functional space. Functional type error estimates (so-called majorants and minorants of deviation from the exact solution) were introduced in [36, 37] and later applied to different mathematical models [33, 30]. They provide guaranteed, sharp, and fully computable upper and lower bounds of errors. This approach, in combination with the IgA approximations generated by tensor-product splines, was proposed and investigated in [22] for elliptic boundary value problems (BVP).

In this paper, we derive functional type a posteriori error estimates for time-dependent problems (see also [35]) in the context of the space-time IgA scheme introduced in [29]. This scheme exploits the time-upwind test function based on the space-time streamline diffusion method (see, e.g., [16, 19, 20]) and the approximations provided by the IgA framework. By exploiting the universality and efficiency of considered error estimates as well as the smoothness of IgA approximations, we aim to construct fully-adaptive, fast and efficient parallel space-time methods that could tackle complicated problems inspired by industrial applications. We also study the numerical properties of newly derived error bounds and compare their performance to the behaviour of the bounds known from [35] on an extensive set of examples.

∗RICAM, Austrian Academy of Sciences, Linz, Austria,[email protected]

†RICAM, Austrian Academy of Sciences, Linz, Austria,[email protected]

‡University of Jyvaskyla, Finland; St. Petersburg Department of V.A. Steklov Institute of Mathematics RAS,[email protected], [email protected]

The outline of the paper is as follows: Section 2 states the problem, discusses its solvability and provides an overview of the existing error control tools for considered initial BVPs (I–BVPs). In Section 3, we deduce new functional type a posteriori error estimates using a stabilised formulation of parabolic I–BVPs. Our anal- ysis is based on a series of transformations performed on a stabilised variational setting; the result of these transformations defines respective generalised solutions. Section 4 presents a stabilised space-time IgA scheme with its main properties along with an overview of main ideas and definitions of the IgA framework. Section 5 is devoted to the algorithmic realisation of an adaptive procedure based on the a posteriori error estimates dis- cussed above. Finally, in Section 6 we present and discuss our numerical results that demonstrate the efficiency of several majorants and the error identity for a wide range of examples.

2 Model problem

Let Q:= Q∪∂Q, Q:= Ω×(0, T), denote a space-time cylinder, where Ω ⊂R^d, d∈ {1,2,3}, is a bounded Lipschitz domain with a boundary∂Ω, and (0, T) is a given time interval with the final timeT, 0< T <+∞.

Here, the boundary∂Qof the space-time cylinderQis defined as∂Q:= Σ∪Σ₀∪Σ_T with Σ =∂Ω×(0, T), Σ₀= Ω× {0}, and Σ_T = Ω× {T}. We discuss our approach to guaranteed error control of space-time approximations within the paradigm of a classical linear parabolic I–BVP: findu:Q→Rsatisfying the parabolic PDE, the boundary condition, and the initial condition

∂tu−∆xu=f in Q, u= 0 on Σ, and u=u0 on Σ0, (1) respectively, where ∂t is a time derivative, ∆x denotes the Laplace operator in space, f ∈ L²(Q) is a given source function, andu0∈H₀¹(Σ0) is prescribed initial data. Here,L²(Q) is a space of square-integrable functions overQequipped with the usual norm and the scalar product denoted respectively by

kvkQ:=kvkL²(Q)= (v, v)¹_L^/₂²_(Q) and (v, w)_Q= (v, w)_L2(Q):=

Z

Q

v(x, t)w(x, t) dxdt, ∀v, w∈L²(Q).

By H^k(Q), k≥1, we denote the spaces of functions that have generalised square-summable derivatives of the orderkwith respect to (w.r.t.) space and time. Next, we introduce the Sobolev spaces

H₀^1,0(Q) :=

u∈L²(Q) :∇_xu∈[L²(Ω)]^d, u _Σ= 0 , V0:=H₀¹(Q) :=

u∈H¹(Q) :u _Σ= 0 , H_0,0¹ (Q) :=

u∈V0 :u _Σ

T = 0 , V_0,0:=H_0,0¹ (Q) :=

u∈V₀ :u _Σ

0= 0 , V₀^∆x,1:=H₀^∆x,1(Q) :=

u∈V₀ : ∆_xu∈L²(Q) , V_0,0^∆x,1:=H₀^∆x,1(Q) :=

u∈V0 : ∆xu∈L²(Q), u _Σ

0 = 0 , and

V_0,0^∇x∂t =H_0,0^∇x∂t(Q) :={w∈V_0,0^∆x :∇_x∂_tw∈L²(Q)}.

(2)

Moreover, later we use the auxiliary Hilbert spaces H^divx,0(Q) :=

y∈[L²(Q)]^d : div_xy∈L²(Q) and H^divx,1(Q) :=

y∈H^divx,0(Q) : ∂_ty∈[L²(Q)]^d for vector-valued functions equipped with the semi-norms defined respectively

kyk²_Hdivx ,0 :=kdivxyk²_Q and kyk²_Hdivx ,1:=kdivxyk²_Q+k∂tyk²_Q. Throughout the paper,CFstands for the constant in the Friedrichs inequality [11]

kwkQ≤CFk∇xwkQ, ∀w∈H₀^1,0(Q). (3) It follows from [27] that, if f ∈ L²(Q) and u₀ ∈ H₀¹(Σ₀), then problem (1) is uniquely solvable in V₀^∆x, and the solutionudepends continuously ontin theH₀¹(Ω)-norm. Moreover, according to [27, Remark 2.2], the normk ∇xu(·, t)k²_Ω is an absolutely continuous function oft∈[0, T] for anyu∈V₀^∆x. Ifu0∈L²(Σ0), then the problem is uniquely solvable in a wider classH₀^1,0(Q), and meets the modified variational formulation

(∇xu,∇xw)Q−(u, ∂tw)Q=:a(u, w) =l(w) := (f, w)Q+ (u0, w)Σ0 (4)

for allw∈H_0,0¹ (Q), where (u0, w)Σ₀ :=R

Σ₀u0(x)w(x,0)dx=R

Ωu0(x)w(x,0)dx. According to well-established arguments (see [26, 27, 42, 6, 7]), without loss of generality, we can ‘homogenise’ the problem, i.e., consider (4) withu0= 0.

Our main goal is to establish fully computable estimates for space-time approximations of this class of problems. For this purpose, we use a functional approach to derive a posteriori error estimates. The first and the simplest forms of such estimates have been derived in [35] for (1). The paper [35] provides the upper bound of the norm

|||e|||²_(ν1,ν2):=ν₁k∇xek²_Q+ν₂kek²_ΣT, ν_i≥0, (5) wheree=u−vis the difference between the exact solutionuand any approximationvin the respective energy classV₀¹. Assuming for simplicity that the initial condition is satisfied exactly, it is shown that for anyv∈V₀¹ approximatingu∈V₀¹and any y∈H^divx,1(Q), we have the following inequality:

|||e|||²_(2−ν,1):= (2−ν)kek²_Q+ kek²_ΣT

≤M^I(v,y;β) := ¹_ν

(1 +β)ky− ∇xvk²_Q+ (1 +¹_β)C_F²kdivxy+f−∂tvk²_Q

, (6)

whereν ∈(0,2] is an auxiliary parameter. The numerical properties of (6) w.r.t. the time-marching and space- time methods are discussed in [32, 31, 17] in the framework of finite-difference and finite-element schemes. The advanced upper error-bound M^II(v,y, η) (valid for the same error norm) introduced in [35] contains an additional auxiliary function η ∈ V0. For the same v and y, as well as any η ∈ V0, any fixed parameters ν ∈(0,2] and γ∈[1,+∞), an alternative majorant has the form

(2−ν)k∇xek²_Q+ (1−_γ¹)kekΣ_T ≤M^II(v,y, η) :=γkηk²_ΣT +ku0−vk²_Σ0−2 (η, u0−uh)Σ₀+ 2F(v, η) +_ν¹n

(1 +β)ky− ∇_xv+∇_xηk²_Q+C_F²(1 +_β¹)kdiv_xy+f−∂_tv−∂_twk²_Qo

, (7)

where

F(v, η) := (∇xv,∇xη) + (∂tv−f, η).

The optimal values of parameters ν, β, and γ are defined uniquely for each of the majorants M^I and M^II, however, in order to keep the notation simpler, we omit the indices I and II.

Furthermore, we note that for the case where u, v∈V₀^∆x, the heat equation (1) imposes the error identity (see [1]):

k∆xek²_Q+k∂tek²_Q+k∇xek²_Σ

T =:|||e|||²_L≡EId²(v) :=k∇x(u0−v)k²_Σ

0+k∆xv+f−∂tvk²_Q. (8) The numerical performance of estimates M^I(v,y) and M^II(v,y, η), and the error identity EId(v) is studied in Section 6.

3 Error majorants

In this section, we derive error majorants of the functional type for a stabilised weak formulation of parabolic I–

BVPs. They provide guaranteed and fully computable upper bounds of the distance between the exact solution uand some approximation v. The functional nature of these majorants allows us to obtain a posteriori error estimates foru∈V_0,0^∆x and any conforming approximationv∈V_0,0^∆x.

We begin by testing (1) with the time-upwind test function

λ w+µ ∂tw, w∈V_0,0^∇x∂t, λ, µ≥0, (9) and arrive at the stabilised weak formulation foru∈V_0,0^∇x∂t, i.e.,

∂_tu, λ w+µ ∂_tw

Q+ ∇xu,∇x(λ w+µ ∂_tw)

Q=:a_s(u, w) =l_s(w) := (f, λ w+µ ∂_tw)_Q, ∀w∈V_0,0^∇x∂t. (10) Then, the errore=u−v is measured in terms of the norm generated by the bilinear formas(u, w), i.e.,

|||e|||²_s,νi:=ν1k∇xek²_Q+ν2k∂tek²_Q+ν3k∇xek²_ΣT +ν4kek²_ΣT, (11) where{νi}i=1,...,4are the positive weights introduced in the derivation process.

To obtain guaranteed error bounds of|||e|||²_s,νi, we apply the method similar to the one developed in [35, 32]

for parabolic I–BVPs. As a starting point, we consider the space of functions V_0,0^∇x∂t (see (2)) equipped with the norm

kwk_V∇x ∂t 0,0

:= sup

t∈[0,T]

k∇xw(·, t)k²_Q+kwk²_V∆x 0,0

,

wherekwk²

V_0,0^∆x :=k∆_xwk²_Q+k∂_twk²_Q, which is dense inV_0,0^∆x. According to [27, Remark 2.2], the normsk·k_V∇x ∂t 0,0

andk · k_V∆x 0,0

are equivalent.

Let{un}^∞_n=1 be a sequence inV_0,0^∇x∂t, which is substituted into the identity (10), i.e.,

as(un, w) = (fn, λ w+µ ∂tw)Q, where fn= (un)_t−∆xun∈L²(Q). (12) Next, we consider the sequence {vn}^∞_n=1 ∈V_0,0^∇x∂t approximating{un}^∞_n=1. By subtractingas(vn, w) from the left- and right-hand side (LHS and RHS, respectively) of (12) and by setting the test function w to be the following differenceen=un−vn ∈V_0,0^∇x∂t, we arrive at the error identity

λk∇xenk²_Q+µk∂tenk²_Q+¹₂(µk∇xenk²_ΣT +λkenk²_ΣT)

=λ

(fn−∂tvn, en)Q−(∇xvn,∇xen)Q

+µ

(fn−∂tvn, ∂ten)Q−(∇xvn,∇x∂ten)Q

. (13) This identity is used to derive the majorants (11) in Theorems 1 and 2.

Theorem 1 For any v∈V_0,0^∆x andy∈H^divx,0(Q), the following estimate holds:

(2−¹_γ) (λk∇xek²_Q+µk∂tek²_Q) +λkek²_ΣT +µk∇xek²_ΣT

=:|||e|||^I,2_s ≤M^I_s,h(v,y;γ, β, α) :=γn

λM^I(v,y;β) +µMe^I(v,y;α)o

, (14) whereM^I(v,y;β) is the majorant defined in (6) withν = 1, i.e.,

M^I(v,y;β) := (1 +β)kr^I_dk²_Q+ (1 +_β¹)C_F²kr^I_eqk²_Q and

Me^I(v,y;α) := (1 +α)kdivxr^I_dk²_Q+ (1 + _α¹)kr^I_eqk²_Q. Here,CF is the Friedrichs constant (3),req andrdare the residuals defined by relations

r^I_eq(v,y) :=f −∂_tv+ div_xy and r^I_d(v,y) :=y− ∇xv, (15) λ, µ >0 are the weights introduced in (9),γ∈₁

2,+∞), andα, β >0.

Proof: First, we modify the RHS of (13) by means of the relation

(divxy, λ en+µ ∂ten)Q+ (y,∇x(λ en+µ ∂ten))Q= 0.

The obtained result can be presented as follows:

λk∇_xe_nk²_Q+µk∂_te_nk²_Q+¹₂(µk∇_xe_nk²_Σ

T +λke_nk²_Σ

T)

=λ (fn−∂tvn+ divxy, en)_Q+ (y− ∇xvn,∇xen)Q

+µ (fn−∂tvn+ divxy, ∂ten)_Q+ (y− ∇xv,∇x∂ten)Q

=λ r^I_eq(vn,y_n), en

Q+ (r^I_d(vn,y_n),∇xen)Q

+µ r^I_eq(vn,y_n), ∂ten

Q+ (r^I_d(vn,y_n),∇x∂ten)Q

. (16) Integrating by parts of the term (r^I_d,∇x∂_te_n)_Q leads to the identity

µ r^I_d,∇x(∂ten)

Q=−µ(divx(y− ∇xvn), ∂ten)Q=−µ(divxy−∆xvn, ∂ten)Q.

Using density arguments, i.e., u_n → u, v_n → v ∈ V_0,0^∆x, and f_n → f ∈ L²(Q), for n → ∞, we arrive at the identity formulated fore=u−v withu, v∈V_0,0^∆x, i.e.,

λk∇_xek²_Q+µk∂_tek²_Q+¹₂(µk∇_xek²_Σ

T +λkek²_Σ

T)

=λ (r_eq, e)_Q+ (r^I_d,∇xe)_Q

+µ (r_eq, ∂_te)_Q−µ(div_xr^I_d, ∂_te)_Q

. (17)

a b Σ₀

ΣT

T

0

Q Σ

x t

ˆ x ˆt

Qb Φ⁻¹

Φ 0 1

1

Figure 1: Mapping of the single-patch reference (parameter) domainQb to the physical single-patch space-time cylinderQ.

By means of the H¨older, Friedrichs, and Young inequalities with positive scalar-valued parametersγ,β, andα,

we deduce estimate (14).

The next theorem assumes higher regularity on the approximation vand the auxiliary functiony.

Theorem 2 For any v∈V_0,0^∇x∂t andy∈H^divx,1(Q), we have the inequality

(2−¹_ζ)(λk∇xek²_Q+µk∂tek²_Q) +µ(1−¹)k∇xek²_ΣT +λkek²_ΣT =:|||e|||^II,2_s ≤M^II_s,h(v,y;ζ, α, , β) := µkrdk²_ΣT +ζ

λ (1 +α) M^I(v,y;β) + (1 +_α¹)^µ_λ²2k∂tr^I_dk²_Q

+µkr^I_eqk²_Q ,

whereM^I(v,y;β)is the majorant defined in (6) withν = 1,CF is the Friedrichs constant in (3),r^I_eq(v, y)and r^I_d(v, y)are the residuals defined in (15),λ, µ >0 are the parameters from (9),ζ∈₁

2,+∞),∈[1,+∞), and β, α >0.

4 Stabilised formulation and its IgA discretisation

For the reader convenience, we recall the general concept of the IgA approach, the definitions of B-splines (NURBS), and their use in geometrical representation of the space-time cylinderQ as well as in construction of the IgA trial spaces, which are used to approximate the solutions that satisfy variational formulation (4).

Throughout the paper,p≥2 denotes the degree of polynomials used for the IgA approximations, whereasn denotes the number of basis functions used to construct aB-spline curve. Aknot-vectoris a non-decreasing set of coordinates in the parameter domain, written as Ξ ={ξ₁, ..., ξ_n+p+1},ξ_i∈R, whereξ₁= 0 andξ_n+p+1= 1.

The knots can be repeated, and multiplicity of thei-th knot is indicated by mi. In what follows, we consider only open knot vectors, i.e., m1 =mn+p+1 = p+ 1. For the one-dimensional parametric domain Qb := (0,1), Kbh := {K}b denotes a locally quasi-uniform mesh, where each element Kb ∈ Kbh is constructed by distinct neighbouring knots. The global size ofKb_his denoted by ˆh:= max

K∈b Kbh{ˆh

Kb}, where ˆh

Kb := diam(K).b

Theunivariate B-spline basis functionsBbi,p:Qb→Rare defined by means of Cox-de Boor recursion formula

Bbi,p(ξ) := _ξ^ξ−ξi

i+p−ξ_iBb_i,p−1(ξ) +_ξ^ξi+p+1−ξ

i+p+1−ξ_i+1Bb_i+1,p−1(ξ), Bbi,0(ξ) :=

(1 if ξi ≤ξ≤ξi+1

0 otherwise , (18)

and are (p−mi)-times continuously differentiable across thei-th knot with multiplicitymi. The multivariate B-splines on the parameter domainQb := (0,1)^d+1, d={1,2,3}, are defined as a tensor-product of the corre- sponding univariate ones. In the multidimensional case, we define the knot-vector dependent on the coordinate direction Ξ^α={ξ₁^α, ..., ξ_n^αα+p^α+1}, ξ_i^α∈R, whereα= 1, ..., d+ 1 indicates the direction (in space or in time).

Furthermore, we introduce a set of multi-indicesI =

i= (i1, ..., id+1) :iα = 1, ..., nα, α= 1, ..., d+ 1 and a multi-indexp:= (p1, ..., pd+1) indicating the order of polynomials. The tensor-product of univariate B-spline basis functions generates multivariate B-spline basis functions

Bb_i,p(ξ) :=

d+1

Y

α=1

Bb_iα,pα(ξ^α), ξ= (ξ¹, ..., ξ^d+1)∈Q.b

The univariate and multivariate NURBS basis functions are defined in the parametric domain by means of B-spine basis functions, i.e., for the givenpand anyi∈ I, the NURBS basis functionsRbi,p:Qb→Rare defined as

Rbi,p(ξ) :=^wi_W^Bbi,p_(ξ)^(ξ). (19) Here, W(ξ) is the weighting functionW(ξ) :=P

i∈Iw_iBb_i,p(ξ), wherew_i ∈R⁺.

The physical space-time domain Q ⊂ R^d+1 is defined by the geometrical mapping Φ : Qb → Q of the parametric domainQb:= (0,1)^d+1:

Q:= Φ(Q)b ⊂R^d+1, Φ(ξ) :=X

i∈I

Rbi,p(ξ)Pi, (20)

where {Pi}i∈I ∈ R^d+1 are the control points. For simplicity, we assume the same polynomial degree for all coordinate directions, i.e.,p_α=pfor allα= 1, ..., d+ 1. By means of geometrical mapping (20), the meshKh

discretisingQis defined asK_h:=

K= Φ(K) :b Kb ∈Kb_h . The global mesh size is denoted by h:= max

K∈Kh

{h_K}, h_K :=k∇ΦkL^∞(K)hˆ

Kb. (21)

Moreover, we assume thatK_h is a quasi-uniform mesh, i.e., there exists a positive constantC_u independent of h, such thath_K ≤h≤C_uh_K.

The finite dimensional spaces V_h on Q are constructed by a push-forward of the NURBS basis functions, i.e.,

Vh:= span

φh,i:=Rbi,p◦Φ^ư1

i∈I, (22)

where the geometrical mapping Φ is invertible inQ, with smooth inverse on each elementK∈ Kh, see, e.g., see [38, 3]. The subspace

V_0h:=V_h∩V_0,0^∆x,1

is introduced for the functions satisfying homogeneous boundary and initial conditions.

A stable space-time IgA scheme for (1) was presented and analysed in [29], where the authors proved its efficiency for fixed and moving spatial computational domains. In our analysis, we use spline bases of sufficiently high order, so thatvh∈V0h ⊂V_0,0^∆x,1. In order to provide an efficient discretisation method, we consider (10), whereλ= 1 andµ=δh=θhin (9) with some positive parameterθand the global mesh-sizehdefined in (21).

It implies the stabilised space-time IgA scheme: finduh∈V0h satisfying (∂_tu_h, v_h+δ_h∂_tv_h)_Q+ ∇_xu_h,∇_x(v_h+δ_h∂_tv_h)

Q =:a_s,h(u_h, v_h) =l_s,h(v_h) := (f, v_h+δ_h∂_tv_h)_Q. (23) for allvh∈V0h. TheV0h-coercivity ofah(·,·) :V0h×V0h→Rw.r.t. the norm

|||v_h|||²_s,h:=k∇_xv_hk²_Q+δ_hk∂_tv_hk²_Q+kv_hk²_Σ

T +δ_hk∇_xv_hk²_Σ

T (24)

follows from [29, Lemma 1]. As was noted in [29], coercivity implies that the IgA solutionu_h ∈ V_0h of (23) is unique. Moreover, since the IgA scheme (23) is posed in the finite dimensional spaceV_0h, uniqueness yields existence of the solution inV_0h. Moreover, following [29, 28], we can show boundedness of the bilinear form in (23) with respect to appropriately chosen norms. Combining coercivity and boundedness properties ofa_s,h(·,·) with the consistency of the scheme and approximation results for IgA spaces implies a corresponding a priori error estimate presented in Theorem 3 below.

Theorem 3 Let u∈H₀^s(Q) :=H^s(Q)∩H₀^1,0(Q),s∈N,s≥2, be the exact solution to (4), and letu_h∈V_0h be the solution to (23)with some fixed parameterθ. Then, the following a priori error estimate

kuưu_hks,h≤C h^rư1kukH^r(Q) (25) holds, wherer= min{s, p+ 1}, andC >0 is a generic constant independent ofh.

Proof: See, e.g., [29, Theorem 8].

Corollary 1 presents a posteriori error majorants forλ= 1 andµ=δ_h, whereδ_h=θ h,θ >0.

Corollary 1 (i) Ifv∈V_0,0^∆x andy∈H^divx,0(Q), Theorem 1 yields the estimate (2ư¹_γ) (k∇xek²_Q+δ_hk∂tek²_Q) +kek²_ΣT +δ_hk∇xek²_ΣT =:|||e|||^I,2_s,h

≤M^I_s,h(v,y;γ, β, α) :=γ

M^I(v,y;γ, β) +δhMe^I(v,y;α)

, (26)

whereM^I andMe^I are defined in Theorem 1 and the best β andαare given by relations β= ^CFkr

I eqk_Q

kr^I_dkQ and α= ^kr

I eqk_Q kdivxr^I_dkQ. A particularly useful form of (26)follows γ= 1, i.e.,

k∇xek²_Q+δ_hk∂tek²_Q+kek²_ΣT +δ_hk∇xek²_ΣT =:|||e|||²_s,h≤M^I_s,h(v,y;α, β) := M^I(v,y;β) +δ_hMe^I(v,y;α).

(ii) Ifv∈V_0,0^∇x∂t andy∈H^divx,1(Q), then Theorem 2 yields

(2ư¹_ζ)(k∇xek²_Q+δhk∂tek²_Q) + (1ư¹)k∇xek²_ΣT +δhkek²_ΣT =:|||e|||^II,2_s,h≤M^II_s,h(v,y;ζ, β, α, ^I) := δhkr^I_dk²_Σ

T +ζ

(1 +α) M^I(v,y;γ, β) + (1 + _α¹)δ²_hk∂tr^I_dk²_Q+δhkr^I_eqk²_Q

, (27) where the optimal parameters are given by relations

β= ^CFkr

I eqkQ

kr^I_dkQ and α=r ^{δhk∂trIdkQ} (1+β)kr^I_dk²_Q+ 1+1 β

C_F²kr^I_eqk²_Q

.

Forζ= 1and= 2, we obtain

k∇xek²_Q+δhk∂tek²_Q+kek²_ΣT +^δ₂^hk∇xek²_ΣT ≤M^II_s,h(v,y;β, α) := δhkr^I_dk²_ΣT +ζ

(1 +α) M^I(v,y;β) + (1 +_α¹) +δ²_hk∂tr^I_dk²_Q+δhkr^I_eqk²_Q

. (28)

In both (i) and (ii),r^I_d andr^I_eq are defined in (15),CF is the Friedrichs constant in (3),δh is the discretisation parameter,γ, ζ∈1

2,+∞),∈[1,+∞), andβ, α >0.

5 Algorithmic realisation

In this section, we discuss the IgA discretisation of the variational formulation presented above as well as the estimates that control the reconstructed approximations quality. We also suggest efficient algorithms for the reconstruction of a posteriori error bounds. The numerical examples presented in Section 6 demonstrate computational properties of the majorants that follow from [35], of the error identity EId, and of the error bounds exposed in Section 3.

5.1 Computation of the majorants in the IgA framework

We consider the approximations

uh∈V0h:=Vh∩V_0,0^∆x,1, where Vh ≡ S_h^p :=

φh,i := ˆS_h^p◦Φ^ư1 are generated with the NURBS of degreep= 2. Due to the restriction on knots-multiplicity of ˆS_h^p in the framework of one-patch domains,uh∈C^pư1 is automatically provided. It is important to note that the scope of this paper is limited to a single-patch domain since it is important to first fully analyse the behaviour of the error-estimation tool in a simplified setting. The extension of this simpler setting to a widely used in practice multi-patch case, in which the physical domain is decomposed into several simple patches, will be a focus of the subsequent paper.

Then approximation has the form

uh(x, t) =uh(x1, ..., xd+1) :=X

i∈I

u_h,iφh,i(x1, ..., xd+1), where u_h :=

u_h,i

i∈I ∈ R^|I| is the vector of degrees of freedom (also called control points in the IgA community) defined by the linear system

Khu_h= fh, Kh:=

as,h(φh,i, φh,j)

i,j∈I, fh:=

ls,h(φh,i)

i∈I. (29)

In the numerical tests presented in Section 6, we analyse the approximation properties ofuhby looking at the convergence of the errore=uưuh measured in terms of the following three norms earlier defined in (5) (with ν = 1), (8), and (24), i.e.,

|||e|||²_(1,1)=|||e|||²:=k∇xek²_Q+kek²_ΣT,

|||e|||²_L:=k∆xek²_Q+k∂tek²_Q+k∇xek²_ΣT, and

|||e|||²_s,h:=k∇_xek²_Q+δ_hk∂_tek²_Q+kek²_Σ

T +δ_hk∇_xek²_Σ

T.

(30)

The majorant for|||e|||² (defined in (6) withν= 1) has the form

M^I(uh,y_h) := (1 +β)ky_hư ∇xuhk²_Q+ (1 +_β¹)C_F²kdivxy_h+fư∂tuhk²_Q = (1 +β) m^I_d+ (1 +_β¹)C_F²m^I_eq, (31) where β >0 and y_h ∈Yh ⊂H^divx,0(Q). The space Yh ≡ S_h^q :=

ψ_h,i :=⊕^d+1Sˆ_h^q◦Φ^ư1 is generated by the push-forward of⊕^d+1Sˆ_h^q, where ˆS_h^q is the space of NURBS of degreeq approximating each ofd+ 1 components ofy_h= y_h⁽¹⁾, . . . , y_h^(d+1)T

. The best estimate follows from the minimisation of M^I(uh,y_h) w.r.t.

y_h(x, t) =y_h(x1, ..., xd+1) = X

i∈I×(d+1)

y_h,iψ_h,i(x1, ..., xd+1).

Here, ψ_h,i are the basis functions generating the spaceY_h, and y

h :=

yh,i

i∈I ∈ R^(d+1)|I| is defined by the linear system

C_F²Divh+βMh

y_h=ưC_F²zh+βgh, (32) where

Divh:=

(divxψ_h,i,divxψ_h,j)Q

(d+1)|I|

i,j=1 , zh:=

fưvt,divxψ_h,j

Q

(d+1)|I|

j=1 , M_h:=

(ψ_h,i,ψ_h,j)_Q(d+1)|I|

i,j=1 , g_h:=

∇_xv,ψ_h,j

Q

(d+1)|I|

j=1 .

(33)

Next, we consider a discretisation of the second form of the majorant M^II(uh,y_h, ηh). For simplicity of exposition, we assume that the initial condition on Σ0 is satisfied exactly, and parametersδ^IIand γare set to 1. In order to make the reconstruction ofηhtransparent and overcome minimisation of M^II(uh,y_h, ηh) w.r.t.ηh, we representηh asηh=whưuh. Here,uhis the approximation at hand obtained by solving (29) andwhis the solution to the same variational problem (12) using wider approximation space

W0h:=Wh∩H₀¹(Ω), Wh≡ S_h^r:=

χh,i:= ˆS_h^r◦χ^ư1 ,

whereS_h^r is the space of NURBS of degreer. As a result, the functionwh can be represented in the form wh(x, t) =wh(x1, ..., xd+1) :=X

i∈I

w_h,iχh,i. Here, w_h:=

w_h,i

i∈I∈R^|I|is the vector of control points ofwhdefined by the linear system

K^(r)_h w_h= f_h^(r), (34)

where K^(r)_h :=

as,h(χh,i, χh,j)

i,j∈I, f_h^(r):=

ls,h(χh,i)

i∈I, andrindicates the degree of splines used to construct the basisχ_h,i. Taking the new representation ofη_h into account, (7) can be reformulated as follows:

k∇xek²_Q≤M^II(uh, wh) =kwhưuhk²_ΣT + 2F(uh, whưuh) + (1 +β) r^II_d

2

Q+C_F²(1 +_β¹) r^II_eq

2

Q, (35) where

r^II_d(uh,y_h, wh) :=y_h+∇xwhư2∇xuh and r^II_eq(y_h, wh) := divxy_h+fư∂twh. Since∂_tw_his approximated by a richer space, the term

r^II_eq(y_h, w_h)

2

Qis expected to be smaller thankreq(y_h, u_h)k²_Q. Therefore, the value of the error bound M^II must be improved. The optimal parameter β is calculated by β:=CFkr^II_eqkQ/kr^II_dkQ.

In [35], it was shown that if wh=uand y_h=∇xu, inequality (35) can be reformulated as follows:

k∇xek²_Q≤M^II(uh, u) :=kuưuhk²_ΣT + 2F(uh, uưuh) + 4 (1 +β)k∇x(uưuh)k²_Q. Moreover, after rearranging the terms of

F(uh, uưuh) =

(∇xu,∇x(uưuh)) + (∂tuưf, uưuh)

+ (∇x(uhưu),∇x(uưuh)) + (∂t(uhưu), uưuh), it is easy to see that the first scalar product on the RHS ofF(uh, uưuh) vanishes. As a result, we obtain

k∇_xek²_Q≤M^II(u_h, u) :=kuưu_hk²_Σ

T + (4 (1 +β)ư2)k∇_x(uưu_h)k²_Qư2(∂_t(uưu_h), uưu_h)

= (4 (1 +β)ư2)k∇x(uưuh)k²_Q. (36)

Thus, we have the following double inequality k∇_xek²_Q ≤M^II(u_h, u)≤C_MII

gapk∇_xek²_Q, C_MII

gap:= (4 (1 +β)ư2), and thereforeC^ư1

M^IIgapM^II(uh, u) can be used for more efficient error indication.

Algorithm 1Reliable reconstruction ofuh(a single refinement step) Input:Kh{discretisation ofQ}

span

φh,i(x1, ..., xd+1) ,i= 1, ...,|I| {Vh-basis}

APPROXIMATE:

• ASSEMBLE the matrix Khand RHS fh :tas(uh)

• SOLVE Khu_h= fh :tsol(uh)

• Reconstructuh(x, t) =uh(x1, ..., xd+1) :=P

i∈Iu_h,iφh,i

Compute the errore=uưuhmeasured in terms of|||e|||,|||e|||s,h, and|||e|||L :te/w(|||e|||) +te/w(|||e|||s,h) +te/w(|||e|||L) ESTIMATE:

• compute M^I(uh,y_h) :tas(y_h) +tsol(y_h) +te/w(M^I)

• compute M^II(uh, wh) :tas(wh) +tsol(wh) +te/w(M^II)

• compute M^I_s,h(uh,y_h) :t_e/w(M^I_s,h)

• compute EId(uh) :te/w(EId)

MARK: Using the marking criteriaM^∗(σ), select elementsKof the meshKhthat must be refined REFINE: Execute the refinement strategy:Kh_ref =R(Kh)

Output:Kh_ref {refined discretisation of Ω}

5.2 Algorithms

Next, we concentrate on the algorithms providing an adaptive procedure based on the a posteriori error estimates presented above. A reliableuh-approximation procedure is summarised in Algorithm 1. We assume thatf,u0, andQin (1) are given. As an input to Algorithm 1, the initial (or obtained on a previous refinement step) mesh Kh discretising the space-time cylinderQis provided. As an output, Algorithm 1 returns a refined version of the mesh denoted byKhref. Overall, the algorithm is structured according to the classic block-chain

APPROXIMATE→ESTIMATE→MARK→REFINE.

The APPROXIMATE step involves assembling of the system of the IgA solution u_h, i.e., the matrix K_h and RHS fh in (29), and solving it with sparse direct LU factorisations (like Eigen SparseLU [14] that is used in our numerical example). Such a choice of a solver is made in order to provide a fair comparison of time spent on solving (29), (32), and (34). On coarser grids, sparse direct solvers are quite efficient. However, on finer grids, iterative solvers like multigrid become more and more efficient in a nested iteration setting, where one can use the interpolated coarse grid solution as an initial guess on the next, adaptively refined grid, see, e.g., [2, 4, 15]. The time spent on assembling and solving sub-procedures is tracked and saved in the vectors tas(uh) and tsol(uh), respectively. This notation is used in the upcoming examples to analyse the efficiency of Algorithm 1 and to compare the computational costs for its subroutines. After the APPROXIMATE step, the error contained in uh is evaluated in terms of several norms defined in (30), i.e., |||e|||, |||e|||s,h, and |||e|||L. To measure the time for element-wise (e/w) assembling of the latter quantities, we uset_e/w(|||e|||),t_e/w(|||e|||s,h), and t_e/w(|||e|||L), respectively.

The next ESTIMATE step focuses on the reconstruction of the global estimates M^I(uh,y_h), M^II(uh, wh), and Ms,h(uh,y_h), as well as the error identity EId. The time spent on each of the error estimators is measured in the same way, for instance,t_as(y_h),t_sol(y_h), andt_e/w(M^I) correspond to the times required to assemble system (32), solve it, and evaluate e/w contributions of M^I(uh,y_h). Analogously, since M^II(uh,y_h) depends onwh, we store in tas(wh) the time corresponding to the assembling of system (34) and intsol(y_h) the time spent on (34). Element- wise evaluation costs are tracked inte/w(M^II(uh,y_h)). The reconstruction of Ms,h(uh,y_h) as well as EId narrows down to their e/w assembly since they do not have to be optimised and can be directly computed. Therefore, the time-expenses are saved in t_e/w(M_s,h(u_h,y_h)) as well ast_e/w(EId). A detailed description of the majorant M^I(uh,y_h) calculation procedure is presented in Algorithm 2, whereas the steps of M^II(wh)-reconstruction are described in Algorithm 3.

In the third chain-block MARK, we use a marking criterion denoted byM^∗(σ). It provides an algorithm for defining the thresholdS∗ for selecting thoseK∈ Kh for further refinement that satisfies the criterion

m^I,2_d,K ≥S∗(M^∗(σ)), K∈ Kh.

Having reconstructed EId(uh) in addition to m^I_d(uh,y_h), which is defined by one term of M^I(uh,y_h), we have a variety of different error indicators to base the mesh refinement strategy on. In the open source C++ library G+Smo [21] used for carrying out the numerical examples presented further, several marking strategies are considered. In particular, the marking based on ‘absolute threshold’ is denoted as GARU (an abbreviation for ‘greatest appearing residual utilisation’), the ‘relative threshold’ is denoted as PUCA (which stands for

‘percent-utilising cutoff ascertainment’), and the most widely used bulk marking (also known as the D¨orfler marking [10]) is denoted by BULK. In further examples, we mainly use the latter marking criterion. In the case of uniform refinement, all elements ofKh are marked for refinement (i.e,σ= 0). If the numerical IgA scheme is implemented correctly, the error is supposed to decrease at least as O(h^p), which is verified throughout the numerical tests in Section 6.

Finally, on the last REFINE step, we apply the refinement algorithm Rto those elements that have been selected on the MARK level. Since the THB-splines are based on subdomains of different hierarchical levels, the procedureRincreases the level of subdomains by applying the dyadic cell refinement.

In the following, we concentrate on the structure of Algorithm 2, which clarifies the ESTIMATE step of Algorithm 1 in the context of functional type error estimates. On the Input step, the algorithm receives the approximate solution uh reconstructed by the IgA scheme. Moreover, since the majorant is minimised w.r.t.

the vector-valued variable y_h ∈ Y_h, the collection of basis functions generating the space Y_h := span ψ_h,i , i = 1, ...,(d+ 1)|I| is provided. The last input parameter N_maj^it defines the number of the optimisation loops executed to obtain a good enough minimiser of M^I. According to the tests performed in [34, 32, 17], one or two iterations is usually sufficient to achieve the reasonable accuracy of error majorant. Another criterion to exit the cycle earlier and, therefore, minimise the computational costs of the error control, can be the condition that the ratio (1 + ¹_β)C_F²m^I,2_eq/(1 +β) m^I,2_d is small enough. In this case, the efficiency index is automatically close to one. When the calculation of M^I is followed up by the reconstruction of M^II, we consider onlyN_maj^it = 1 iteration. In addition to M^I and M^II, we evaluate the majorant M^I_s,h specifically derived in Theorem 1 for the stabilised scheme (12) and the control of the error|||e|||s,h.

We emphasise that both matrices Divhand Mh, as well as vectors zhand gh, are assembled only once. The loop iteratesN_maj^it times such that each time the optimaly⁽ⁿ⁾_h andβ^I,(n)are reconstructed. In our implementa- tion, the optimality system for the flux (see (32)) is solved by the sparse direct LDL^T Cholesky factorisations.

The time spent on ASSEMBLE and SOLVE steps w.r.t. the system (32) is measured by tas(y_h) andtsol(y_h) respectively and compared tot_as(u_h) andt_sol(u_h) in forthcoming numerical examples. It is crucial to note that the matrices Div_h and M_h have block structure (of (d+ 1)×(d+ 1) blocks) due to the properties of the ap- proximation spacesV_h andY_h. Moreover, since Div_h, M_h, r_h and g_h are generated by the scalar product of the derivatives or divergence w.r.t. spatial coordinates only, (d+ 1, d+ 1)-th block of Div_h is zero as well as the (d+ 1)-th block of the RHS of (32), i.e.,

C_F²

Div^(d)_h 0

0 0

+β

"

M^(d)_h 0 0 M_h⁽¹⁾

#!

·

"

y^(d)_h y⁽¹⁾_h

#

=−C_F² z^(d)_h

0

+β g^(d)_h

0

,

where (1)-block corresponds to the time variable. This resolves into the vectorYh with zero (d+ 1)-th block, which in turn allows us to solve the system composed only of spatial blocks. Besides the computational costs related to the assembling and solving of (29) and (32), we measure the time spent on the e/w evaluation of all the majorants.

Analogously to the selection ofqfor the spaceYh, we letr=p+l,l∈N⁺. At the same time, we use a coarser mesh KLh,L∈N⁺ in order to recoverwh. For the reader convenience, we collect the notation related to the spaces parametrisation in Table 1. The sequence of steps of thewh-approximation, as well as M^II-reconstruction corresponding to it, are presented in Algorithm 3. Its structure is similar to the structure of Algorithm 2 with the exception that the free variable of M^II(vh,y_h, wh) is a scalar function and we solve system (34) to reconstruct the degrees of freedom (d.o.f.) ofwh only once.

Evaluation of the error identity does not require any optimisation techniques. Therefore, it can be computed straightforwardly by using

EId²(u_h) :=k∇_x(u₀−u_h)k²_Σ

0+k∆_xu_h+f−∂_tu_hk²_Q

without any overhead in time performance. Time spent on the element-wise assembly of EId is tracked int_e/w(EId).