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.



There are no comments yet.


page 1

page 2

page 3

page 4


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...

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

We consider a space-time variational formulation of the second-order wav...

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

Variational formulations of time-dependent PDEs in space and time yield ...

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) ...

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

An exact arithmetic, memory efficient direct solution method for finite ...

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 ...

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


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


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


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

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

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

are endowed with the Hilbertian norms

with the usual Bochner-Sobolev norms

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

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

[25, 26], see also [21].

Theorem 2.1.

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

with a constant Furthermore, the solution operator

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

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

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

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


of piecewise multilinear, continuous functions, i.e.

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


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

where the given function is represented by


with the eigenfunctions

and eigenvalues

, satisfying

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

is norm preserving, bijective and fulfills the coercivity property


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


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


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

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


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

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

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

fulfills the space-time error estimates


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

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


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


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


with the system matrix

where and denote spatial mass and stiffness matrices given by


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

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

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

holds true with where

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

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

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

  • For regular matrices , we have

  • The mixed-product property

    is valid.

  • It holds true that


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

For a given vector , define the matrix

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


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.

Ordering and Factorization Steps Solving Step Memory
Table 1: Summary of the complexity of sparse direct solver for sparse matrices coming from structured grids.

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


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

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

gives the representations

Hence, the global linear system is equivalent to solving

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


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


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