# Efficient Direct Space-Time Finite Element Solvers for Parabolic Initial-Boundary Value Problems in Anisotropic Sobolev Spaces

We consider a space-time variational formulation of parabolic initial-boundary value problems in anisotropic Sobolev spaces in combination with a Hilbert-type transformation. This variational setting is the starting point for the space-time Galerkin finite element discretization that leads to a large global linear system of algebraic equations. We propose and investigate new efficient direct solvers for this system. In particular, we use a tensor-product approach with piecewise polynomial, globally continuous ansatz and test functions. The developed solvers are based on the Bartels-Stewart method and on the Fast Diagonalization method, which result in solving a sequence of spatial subproblems. The solver based on the Fast Diagonalization method allows to solve these spatial subproblems in parallel leading to a full parallelization in time. We analyze the complexity of the proposed algorithms, and give numerical examples for a two-dimensional spatial domain, where sparse direct solvers for the spatial subproblems are used.

## Authors

• 12 publications
• 6 publications
03/07/2021

### Numerical results for an unconditionally stable space-time finite element method for the wave equation

In this work, we introduce a new space-time variational formulation of t...
02/15/2021

### Higher-Order Space-Time Continuous Galerkin Methods for the Wave Equation

We consider a space-time variational formulation of the second-order wav...
12/20/2019

### Matrix oriented reduction of space-time Petrov-Galerkin variational problems

Variational formulations of time-dependent PDEs in space and time yield ...
07/12/2021

### Parallel Element-based Algebraic Multigrid for H(curl) and H(div) Problems Using the ParELAG Library

This paper presents the use of element-based algebraic multigrid (AMGe) ...
02/11/2020

### Direct Domain Decomposition Method (D3M) for Finite Element Electromagnetic Computations

An exact arithmetic, memory efficient direct solution method for finite ...
01/07/2022

### A Direct Parallel-in-Time Quasi-Boundary Value Method for Inverse Space-Dependent Source Problems

Inverse source problems arise often in real-world applications, such as ...
07/26/2021

### A Parallel Boundary Element Method for the Electromagnetic Analysis of Large Structures With Lossy Conductors

In this paper, we propose an efficient parallelization strategy for boun...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Parabolic initial-boundary value problems are usually discretized by time-stepping schemes and spatial finite element methods. These methods treat the time and spatial variables differently, see, e.g., [45]

. In addition, the resulting approximation methods are sequential in time. In contrast to these approaches, space-time methods discretize time-dependent partial differential equations without separating the temporal and spatial directions. In particular, they are based on space-time variational formulations. There exist various space-time techniques for parabolic problems, which are based on variational formulations in Bochner-Sobolev spaces, see, e.g.,

[3, 4, 14, 20, 23, 28, 31, 35, 39, 43, 46], or on discontinuous Galerkin methods, see, e.g., [19, 32, 33], or on discontinuous Petrov-Galerkin methods, see, e.g., [13], and the references therein. We refer the reader to [15] and [40] for a comprehensive overview of parallel-in-time and space-time methods, respectively. An alternative is the discretization of space-time variational formulations in anisotropic Sobolev spaces, see, e.g., [8, 24, 36, 41, 48]. These variational formulations allow the complete analysis of inhomogeneous Dirichlet or Neumann conditions, and were used for the analysis of the resulting boundary integral operators, see [6, 9]. Hence, discretizations for variational formulations in anisotropic Sobolev spaces can be used for the interior problems of FEM-BEM couplings for transmission problems.

In this work, the approach in anisotropic Sobolev spaces is applied in combination with a novel Hilbert-type transformation operator , which has recently been introduced in [41, 48]. This transformation operator maps the ansatz space to the test space, and gives a symmetric and elliptic variational setting of the first-order time derivative. The homogeneous Dirichlet problem for the nonstationary diffusion respectively heat equation

 ∂tu(x,t)−Δxu(x,t)=f(x,t)for(x,t)∈Q=Ω×(0,T),u(x,t)=0for(x,t)∈Σ=∂Ω×[0,T],u(x,0)=0forx∈Ω,⎫⎪⎬⎪⎭ (1)

serves as model problem for a parabolic initial-boundary value problem, where is a bounded Lipschitz domain with boundary , is a given terminal time, and is a given right-hand side. With the help of the Hilbert-type transformation operator , a Galerkin finite element method is derived, which results in one global linear system

 Khu––=f––. (2)

When using a tensor-product approach, the system matrix can be represented as a sum of Kronecker products. The purpose of this paper is the development of efficient direct space-time solvers for the global linear system 2, exploiting the Kronecker structure of . Therefore, we apply the Bartels-Stewart method [5] and the Fast Diagonalization method [29] to solve (2), see also [16, 37]

. For both methods, we derive complexity estimates of the resulting algorithms.

The rest of this paper is organized as follows: In Section 2, we consider the space-time variational formulation in anisotropic Sobolev spaces and the Hilbert-type transformation operator with its main properties. In Section 3, we rephrase properties of the Kronecker product and of sparse direct solvers, which are needed for the new space-time solver. Section 4 is devoted to the construction of efficient space-time solvers. Numerical examples for a two-dimensional spatial domain are presented in Section 5. Finally, we draw some conclusions in Section 6.

## 2 Space-Time Method in Anisotropic Sobolev Spaces

In this section, we give the variational setting for the parabolic model problem (1), which is studied in greater detail in [6, 21, 25, 26, 41, 48]. We consider the space-time variational formulation of (1) in anisotropic Sobolev spaces to find such that

 a(u,v)=⟨f,v⟩Q (3)

for all where is a given right-hand side. Here, the bilinear form

 a(u,v):=⟨∂tu,v⟩Q+⟨∇xu,∇xv⟩L2(Q)

for , , is bounded, i.e. there exists a constant such that

 ∀u∈H1,1/20;0,(Q):∀v∈H1,1/20;,0(Q):|a(u,v)|≤C∥u∥H1,1/20;0,(Q)∥v∥H1,1/20;,0(Q),

see [6, Lemma 2.6, p. 505]. The anisotropic Sobolev spaces

 H1,1/20;0,(Q) :=H1/20,(0,T;L2(Ω))∩L2(0,T;H10(Ω)), H1,1/20;,0(Q) :=H1/2,0(0,T;L2(Ω))∩L2(0,T;H10(Ω))

are endowed with the Hilbertian norms

 ∥v∥H1,1/20;0,(Q) :=√∥v∥2H1/20,(0,T;L2(Ω))+∥∇xv∥2L2(Q), ∥w∥H1,1/20;,0(Q) :=√∥w∥2H1/2,0(0,T;L2(Ω))+∥∇xw∥2L2(Q)

with the usual Bochner-Sobolev norms

 ∥v∥H1/20,(0,T;L2(Ω)) := ⎷∥v∥2H1/2(0,T;L2(Ω))+∫T0∥v(⋅,t)∥2L2(Ω)tdt, ∥w∥H1/2,0(0,T;L2(Ω)) := ⎷∥w∥2H1/2(0,T;L2(Ω))+∫T0∥w(⋅,t)∥2L2(Ω)T−tdt,

see [25, 26, 41, 48] for more details. The dual space is characterized as completion of with respect to the Hilbertian norm

 ∥f∥[H1,1/20;,0(Q)]′:=sup0≠w∈H1,1/20;,0(Q)∣∣⟨f,w⟩Q∣∣∥w∥H1,1/20;,0(Q),

where denotes the duality pairing as extension of the inner product in In [6]

, the following existence and uniqueness theorem is proven by a transposition and interpolation argument as in

###### Theorem 2.1.

Let the right-hand side be given. Then, the variational formulation (3) has a unique solution satisfying

 ∥u∥H1,1/20;0,(Q)≤C∥f∥[H1,1/20;,0(Q)]′

with a constant Furthermore, the solution operator

 L:[H1,1/20;,0(Q)]′→H1,1/20;0,(Q),Lf:=u,

is an isomorphism.

For simplicity, we only consider homogeneous Dirichlet conditions, where inhomogeneous Dirichlet conditions can be treated via homogenization as for the elliptic case, see [38, p. 61-62], since for any Dirichlet data , an extension with exists, see [6, 9] for more details.

For a discretization scheme, let the bounded Lipschitz domain be an interval for or polygonal for or polyhedral for For a tensor-product ansatz, we consider admissible decompositions

 ¯¯¯¯Q=¯¯¯¯Ω×[0,T]=Nx⋃i=1¯¯¯¯¯ωi×Nt⋃ℓ=1[tℓ−1,tℓ]

with space-time elements, where the time intervals with mesh sizes are defined via the decomposition

 0=t0

of the time interval . The maximal and the minimal time mesh sizes are denoted by and , respectively. For the spatial domain , we consider a shape-regular sequence of admissible decompositions

 Tν:={ωi⊂Rd:i=1,…,Nx}

of into finite elements with mesh sizes and the maximal mesh size . The spatial elements are intervals for , triangles or quadrilaterals for , and tetrahedra or hexahedra for . Next, we introduce the finite element space

 Q1h(Q):=V1hx,0(Ω)⊗S1ht(0,T) (4)

of piecewise multilinear, continuous functions, i.e.

 V1hx,0(Ω)=span{ψ1j}Mxj=1⊂H10(Ω),S1ht(0,T)=span{φ1ℓ}Ntℓ=0⊂H1(0,T).

In fact, is either the space of piecewise linear, continuous functions on intervals (), triangles (), and tetrahedra (), or is the space of piecewise linear/bilinear/trilinear, continuous functions on intervals (), quadrilaterals (), and hexahedra (). Analogously, for a fixed polynomial degree , we consider the space of piecewise polynomial, continuous functions

 Qph(Q):=Vphx,0(Ω)⊗Spht(0,T). (5)

Using the finite element space (4), it turns out that a discretization of (3) with the conforming ansatz space and the conforming test space is not stable, see [48, Section 3.3]. A possible way out is the modified Hilbert transformation defined by

 (HTu)(x,t):=∞∑i=1∞∑k=0ui,kcos((π2+kπ)tT)ϕi(x),(x,t)∈Q,

where the given function is represented by

 u(x,t)=∞∑i=1∞∑k=0ui,ksin((π2+kπ)tT)ϕi(x),(x,t)∈Q, (6)

with the eigenfunctions

and eigenvalues

, satisfying

 −Δϕi=μiϕi in Ω,ϕi=0 on ∂Ω,∥ϕi∥L2(Ω)=1,i∈N.

This approach was introduced recently in [41] and [48, Section 3.4]. The novel transformation acts on the finite terminal , whereas analogous considerations of an infinite time interval with the classical Hilbert transformation are investigated in [8, 11, 12, 24]. The most important properties of are summarized in the following, see [41, 42, 47, 48]. The map

 HT:H1,1/20;0,(Q)→H1,1/20;,0(Q)

is norm preserving, bijective and fulfills the coercivity property

 ⟨∂tu,HTv⟩Q=12∞∑i=1∞∑k=0(π2+kπ)ui,k⋅vi,k=:⟨u,v⟩H1/20,(0,T;L2(Ω)),F (7)

for functions with expansion coefficients as in (6). Note that the norm induced by the inner product is equivalent to the norm Moreover, the relations

 ∀v∈L2(Q): ⟨v,HTv⟩L2(Q)≥0, ∀s>0:∀v∈Hs0,(0,T;L2(Ω)),v≠0: ⟨v,HTv⟩L2(Q)>0 (8)

hold true. With the modified Hilbert transformation , the variational formulation (3) is equivalent to find such that

 ∀v∈H1,1/20;0,(Q):a(u,HTv)=⟨f,HTv⟩Q. (9)

Hence, unique solvability of the variational formulation (9) follows from the unique solvability of (3), which implies the stability estimate

 ∀u∈H1,1/20;0,(Q):c∥u∥H1,1/20;0,(Q)≤sup0≠v∈H1,1/20;0,(Q)|a(u,HTv)|∥v∥H1,1/20;0,(Q)

with a constant When using some conforming space-time finite element space the Galerkin variational formulation of (9) is to find such that

 (10)

Note that ansatz and test spaces are equal. In [48], the following theorem is proven.

###### Theorem 2.2.

Let be a conforming space-time finite element space and let be a given right-hand side. Then, a unique solution of the Galerkin variational formulation (10) exists. If, in addition, the right-hand side fulfills then the stability estimate

 ∥uh∥H1/20,(0,T;L2(Ω))≤c∥f∥[H1/2,0(0,T;L2(Ω))]′

is true with a constant

Theorem 2.2 states that, under the assumption , any conforming space-time finite element space leads to an unconditionally stable method, i.e. no CFL condition is required. For the choice of the tensor-product space-time finite element space

 Vh=Qph(Q)∩H1,1/20;0,(Q)

from (5), the Galerkin variational formulation (10) to find such that

 ∀vh∈Qph(Q)∩H1,1/20;0,(Q):a(uh,HTvh)=⟨f,HTvh⟩Q

fulfills the space-time error estimates

 ∥u−uh∥H1/20,(0,T;L2(Ω)) ≤chp+1/2, (11) ∥u−uh∥L2(Q) ≤chp+1, (12) |u−uh|H1(Q) ≤chp (13)

with and with a constant for a sufficiently smooth solution of (3) and a sufficiently regular boundary where for the error estimate (13), the sequence of decompositions of is additionally assumed to be globally quasi-uniform, see [41, 48] for details.

In the remainder of this work, we consider , i.e. the tensor-product space of piecewise linear, continuous functions , where analogous results hold true for an arbitrary polynomial degree

So, the number of the degrees of freedom is given by

 dof=Nt⋅Mx.

For an easier implementation, we approximate the right-hand side by

 f≈Q0hf=Nx∑j=1Nt∑ℓ=1fj,ℓψ0jφ0ℓ∈S0hx(Ω)⊗S0ht(0,T) (14)

with coefficients , where is the projection on the piecewise constant functions with and . So, we consider the perturbed variational formulation to find such that

 ∀vh∈Q1h(Q)∩H1,1/20;0,(Q):a(˜uh,HTvh)=⟨Q0hf,HTvh⟩L2(Q). (15)

Note that, for piecewise linear functions, i.e. , the space-time error estimates (11), (12), (13) are not spoilt. Note additionally that, for , a projection on polynomials of degree should be used instead of for preserving the space-time error estimates (11), (12), (13). After an appropriate ordering of the degrees of freedom, the discrete variational formulation (15) is equivalent to the global linear system

 Khu––=˜F––HT (16)

with the system matrix

 Kh=AHTht⊗Mhx+MHTht⊗Ahx∈RNt⋅Mx×Nt⋅Mx,

where and denote spatial mass and stiffness matrices given by

 Mhx[i,j]=⟨ψ1j,ψ1i⟩L2(Ω),Ahx[i,j]=⟨∇xψ1j,∇xψ1i⟩L2(Ω),i,j=1,…,Mx, (17)

and and are defined by

for . Note that the matrix is dense, symmetric and positive definite, see (7), whereas the matrix is dense, nonsymmetric and positive definite, see (8

). Additionally, the vector of the right-hand side in (

16) is given by

 ˜F––HT:=(˜f––1,…,˜f––Mx)⊤∈RNt⋅Mx

with the vectors ,   , where, with the help of (14),

 ˜f––i[k]:=⟨Q0hf,ψ1iHTφ1k⟩L2(Q)=Nx∑j=1Nt∑ℓ=1fj,ℓ⟨ψ0j,ψ1i⟩L2(Ω)⟨φ0ℓ,HTφ1k⟩L2(0,T),

To assemble the vector of the right-hand side in (16), the relation

 ˜f––i[k]:=˜F[i,k],i=1,…,Mx,k=1,…,Nt,

holds true with where

 M1,0hx[i,j] :=⟨ψ0j,ψ1i⟩L2(Ω), i=1,…,Mx,j=1,…,Nx, F[j,ℓ] :=fj,ℓ, j=1,…,Nx,ℓ=1,…,Nt, CHTht[k,ℓ] :=⟨φ0ℓ,HTφ1k⟩L2(0,T), k=1,…,Nt,ℓ=1,…,Nt.

## 3 Preliminaries for the Space-Time Solvers

In this section, some properties of the Kronecker product and direct solvers, which are needed in Section 4, are summarized.

### 3.1 Kronecker Product

In this subsection, some basic properties of the Kronecker product are stated, see, e.g., [18, 37]. Let ,   and be given matrices for The Kronecker product is defined as the matrix

 A⊗B:=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝A[1,1]BA[1,2]B⋯A[1,NA]BA[2,1]BA[2,2]B⋯A[2,NA]B⋮⋮⋱⋮A[NA,1]BA[NA,2]B⋯A[NA,NA]B⎞⎟ ⎟ ⎟ ⎟ ⎟⎠∈CNA⋅NB×NA⋅NB.

Furthermore, the vectorization of a matrix converts the matrix into a column vector, i.e. we define

 vec(X):=(X[1,1],X[2,1],…,X[NA,1],X[1,2],…,X[NA,NB])⊤∈CNA⋅NB×1.

In the remainder of this work, we use the following properties of the Kronecker product and the vectorization of a matrix:

• For the conjugate transposition and transposition, it holds true that

 (A⊗B)∗=A∗⊗B∗ and (A⊗B)⊤=A⊤⊗B⊤.
• For regular matrices , we have

 (A⊗B)−1=A−1⊗B−1.
• The mixed-product property

 (A⊗B)(C⊗D)=(AC)⊗(BD)

is valid.

• It holds true that

 vec(AXB)=(B⊤⊗A)vec(X), (18)

where also in the case of complex matrices, only the transposition is applied.

For a given vector , define the matrix

 V:=(v–1v–2⋯v–NB)∈CNA×NB

with given by for , i.e. Then, the equality (18) yields

 (B⊤⊗A)v–=(B⊤⊗A)vec(V)=vec(AVB). (19)

### 3.2 Sparse Direct Solver

In this subsection, we repeat some properties of sparse direct solver, like left-looking/right-looking/multifrontal methods, for solving linear systems , see, e.g., [7, 10, 17, 22, 27, 30, 34]. Here, is a sparse matrix coming from finite element/difference discretizations of a physical domain , , like the spatial mass or stiffness matrix (17). Sparse direct solvers exploit the sparsity pattern of the system matrix , and are based on a divide-and-conquer technique, which can be interpreted as procedure of subdividing the physical domain , which leads also to a subdivision of the degrees of freedom. Usually, the following steps have to be applied in such a method:

1. Ordering Step, e.g., minimum degree or nested dissection methods,

2. Symbolic Factorization Step,

3. Numerical Factorization Step,

4. Solving Step.

For structured grids, the complexity of these methods is summarized in Table 1.

There are several open-source software packages for sparse direct solvers. In this paper, we use only the sparse direct solver MUMPS 5.3.3

[1, 2].

## 4 Space-Time Solvers

In this section, efficient solvers for the large-scale space-time system (16) are developed. Our new solver is based on [19, Section 3] and [44, Section 4], where analogous results are derived for methods in isogeometric analysis. In greater detail, we state solvers for the global linear system

 (AHTht⊗Mhx+MHTht⊗Ahx)u––=˜F––HT (20)

given in (16) with the symmetric, positive definite matrices , , and the nonsymmetric, positive definite matrix . Since (20) is a (generalized) Sylvester equation, we can apply the Bartels-Stewart method [5] with real- or complex-Schur decomposition and the Fast Diagonalization method [29] to solve (20), see also [16, 37]. In all three cases, the matrix pencil is decomposed in the form

 (AHTht)−1MHTht=XtZtX−1t

with real, regular matrices , where is an upper (quasi-)triangular matrix, or complex, regular matrices , where is an upper triangular or diagonal matrix. Defining

 Yt:=(AHThtXt)−1

gives the representations

 AHTht=Y−1tX−1t and MHTht=AHTht=Y−1tX−1tXtZtX−1t=Y−1tZtX−1t.

Hence, the global linear system is equivalent to solving

 (Y−1t⊗IMx)(INt⊗Mhx+Zt⊗Ahx)(X−1t⊗IMx)u––=˜F––HT

with the identity matrices and . Thus, the solution of (20) is given by

 u––=(Xt⊗IMx)(INt⊗Mhx+Zt⊗Ahx)−1(Yt⊗IMx)˜F––HT. (21)

The first step in (21) is the calculation of the vector

 g–:=(g–1,g–2,…,g–Nt)⊤:=(Yt⊗IMx)˜F––HT=vec(^F(AHTht)−1X−⊤t)∈CNt⋅Mx (22)

with a matrix corresponding to the relation (19), satisfying where The second step in (21) is to solve the linear system

 (INt⊗Mhx+Zt⊗Ah