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

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 $H_{T}$ , which has recently been introduced in [41, 48]. This transformation operator $H_{T}$ 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

\begin{matrix} \partial_{t} u (x, t) - Δ_{x} u (x, t) & = & f (x, t) & for (x, t) \in Q = Ω \times (0, T), u (x, t) & = & 0 & for (x, t) \in Σ = \partial Ω \times [0, T], u (x, 0) & = & 0 & for x \in Ω, \end{matrix} ⎫ ⎪ ⎬ ⎪ ⎭

(1)

serves as model problem for a parabolic initial-boundary value problem, where $Ω \subset R^{d},$ $d = 1, 2, 3,$ is a bounded Lipschitz domain with boundary $\partial Ω$ , $T > 0$ is a given terminal time, and $f$ is a given right-hand side. With the help of the Hilbert-type transformation operator $H_{T}$ , a Galerkin finite element method is derived, which results in one global linear system

K_{h} u - - = f - - .

(2)

When using a tensor-product approach, the system matrix $K_{h}$ 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 $K_{h}$ . 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 $H_{T}$ 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 $u \in H_{0; 0,}^{1, 1 / 2} (Q)$ such that

a (u, v) = {⟨ f, v ⟩}_{Q}

(3)

for all $v \in H_{0;, 0}^{1, 1 / 2} (Q),$ where $f \in [H_{0;, 0}^{1, 1 / 2} (Q)]^{'}$ is a given right-hand side. Here, the bilinear form $a (\cdot, \cdot) : H_{0; 0,}^{1, 1 / 2} (Q) \times H_{0;, 0}^{1, 1 / 2} (Q) \to R,$

a (u, v) := {⟨ \partial_{t} u, v ⟩}_{Q} + {⟨ \nabla_{x} u, \nabla_{x} v ⟩}_{L^{2} (Q)}

for $u \in H_{0; 0,}^{1, 1 / 2} (Q)$ , $v \in H_{0;, 0}^{1, 1 / 2} (Q)$ , is bounded, i.e. there exists a constant $C > 0$ such that

\forall u \in H_{0; 0,}^{1, 1 / 2} (Q) : \forall v \in H_{0;, 0}^{1, 1 / 2} (Q) : | a (u, v) | \leq C {∥ u ∥}_{H_{0; 0,}^{1, 1 / 2} (Q)} {∥ v ∥}_{H_{0;, 0}^{1, 1 / 2} (Q)},

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

	$H_{0; 0,}^{1, 1 / 2} (Q)$	$:= H_{0,}^{1 / 2} (0, T; L^{2} (Ω)) \cap L^{2} (0, T; H_{0}^{1} (Ω)),$
	$H_{0;, 0}^{1, 1 / 2} (Q)$	$:= H_{, 0}^{1 / 2} (0, T; L^{2} (Ω)) \cap L^{2} (0, T; H_{0}^{1} (Ω))$

are endowed with the Hilbertian norms

	$∥ v ∥_{H_{0; 0,}^{1, 1 / 2} (Q)}$	$:= \sqrt{∥ v ∥_{H_{0,}^{1 / 2} (0, T; L^{2} (Ω))}^{2} + {∥ \nabla_{x} v ∥}_{L^{2} (Q)}^{2}},$
	$∥ w ∥_{H_{0;, 0}^{1, 1 / 2} (Q)}$	$:= \sqrt{∥ w ∥_{H_{, 0}^{1 / 2} (0, T; L^{2} (Ω))}^{2} + {∥ \nabla_{x} w ∥}_{L^{2} (Q)}^{2}}$

with the usual Bochner-Sobolev norms

	$∥ v ∥_{H_{0,}^{1 / 2} (0, T; L^{2} (Ω))}$	$:= \sqrt{{∥ v ∥}_{H^{1 / 2} (0, T; L^{2} (Ω))}^{2} + \int_{0}^{T} \frac{{∥ v (\cdot, t) ∥}_{L^{2} (Ω)}^{2}}{t} d t},$
	$∥ w ∥_{H_{, 0}^{1 / 2} (0, T; L^{2} (Ω))}$	$:= \sqrt{{∥ w ∥}_{H^{1 / 2} (0, T; L^{2} (Ω))}^{2} + \int_{0}^{T} \frac{{∥ w (\cdot, t) ∥}_{L^{2} (Ω)}^{2}}{T - t} d t},$

see [25, 26, 41, 48] for more details. The dual space $[H_{0;, 0}^{1, 1 / 2} (Q)]^{'}$ is characterized as completion of $L^{2} (Q)$ with respect to the Hilbertian norm

∥ f ∥_{[H_{0;, 0}^{1, 1 / 2} (Q)]^{'}} := sup 0 \neq w \in H_{0;, 0}^{1, 1 / 2} (Q) \frac{∣ ∣ ⟨ f, w ⟩_{Q} ∣ ∣}{∥ w ∥_{H_{0;, 0}^{1, 1 / 2} (Q)}},

where $⟨ \cdot, \cdot ⟩_{Q}$ denotes the duality pairing as extension of the inner product in $L^{2} (Q) .$ 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 $f \in [H_{0;, 0}^{1, 1 / 2} (Q)]^{'}$ be given. Then, the variational formulation (3) has a unique solution $u \in H_{0; 0,}^{1, 1 / 2} (Q),$ satisfying

{∥ u ∥}_{H_{0; 0,}^{1, 1 / 2} (Q)} \leq C {∥ f ∥}_{[H_{0;, 0}^{1, 1 / 2} (Q)]^{'}}

with a constant $C > 0.$ Furthermore, the solution operator

L : [H_{0;, 0}^{1, 1 / 2} (Q)]^{'} \to H_{0; 0,}^{1, 1 / 2} (Q), L f := 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 $g \in H^{1 / 2, 1 / 4} (Σ)$ , an extension $u_{g} \in H_{; 0,}^{1, 1 / 2} (Q)$ with $u_{g | Σ} = g$ exists, see [6, 9] for more details.

For a discretization scheme, let the bounded Lipschitz domain $Ω \subset R^{d}$ be an interval $Ω = (0, L)$ for $d = 1,$ or polygonal for $d = 2,$ or polyhedral for $d = 3.$ For a tensor-product ansatz, we consider admissible decompositions

¯ ¯¯ ¯ Q = ¯ ¯¯ ¯ Ω \times [0, T] = N_{x} ⋃ i = 1_{i} \times N_{t} ⋃ ℓ = 1 [t_{ℓ - 1}, t_{ℓ}]

with $N := N_{x} \cdot N_{t}$ space-time elements, where the time intervals $(t_{ℓ - 1}, t_{ℓ})$ with mesh sizes $h_{t, ℓ} = t_{ℓ} - t_{ℓ - 1}$ are defined via the decomposition

0 = t_{0} < t_{1} < t_{2} < \dots < t_{N_{t} - 1} < t_{N_{t}} = T

of the time interval $(0, T)$ . The maximal and the minimal time mesh sizes are denoted by $h_{t} := h_{t, max} := {max}_{ℓ} h_{t, ℓ}$ and $h_{t, min} := {min}_{ℓ} h_{t, ℓ}$ , respectively. For the spatial domain $Ω$ , we consider a shape-regular sequence $(T_{ν})_{ν \in N}$ of admissible decompositions

T_{ν} := {ω_{i} \subset R^{d} : i = 1, \dots, N_{x}}

of $Ω$ into finite elements $ω_{i} \subset R^{d}$ with mesh sizes $h_{x, i}$ and the maximal mesh size $h_{x} := {max}_{i} h_{x, i}$ . The spatial elements $ω_{i}$ are intervals for $d = 1$ , triangles or quadrilaterals for $d = 2$ , and tetrahedra or hexahedra for $d = 3$ . Next, we introduce the finite element space

Q_{h}^{1} (Q) := V_{h_{x}, 0}^{1} (Ω) \otimes S_{h_{t}}^{1} (0, T)

(4)

of piecewise multilinear, continuous functions, i.e.

V_{h_{x}, 0}^{1} (Ω) = span {ψ_{j}^{1}}_{j = 1}^{M_{x}} \subset H_{0}^{1} (Ω), S_{h_{t}}^{1} (0, T) = span {φ_{ℓ}^{1}}_{ℓ = 0}^{N_{t}} \subset H^{1} (0, T) .

In fact, $V_{h_{x}, 0}^{1} (Ω)$ is either the space $S_{h_{x}}^{1} (Ω) \cap H_{0}^{1} (Ω)$ of piecewise linear, continuous functions on intervals ( $d = 1$ ), triangles ( $d = 2$ ), and tetrahedra ( $d = 3$ ), or $V_{h_{x}, 0}^{1} (Ω)$ is the space $Q_{h_{x}}^{1} (Ω) \cap H_{0}^{1} (Ω)$ of piecewise linear/bilinear/trilinear, continuous functions on intervals ( $d = 1$ ), quadrilaterals ( $d = 2$ ), and hexahedra ( $d = 3$ ). Analogously, for a fixed polynomial degree $p \in N$ , we consider the space of piecewise polynomial, continuous functions

Q_{h}^{p} (Q) := V_{h_{x}, 0}^{p} (Ω) \otimes S_{h_{t}}^{p} (0, T) .

(5)

Using the finite element space (4), it turns out that a discretization of (3) with the conforming ansatz space $Q_{h}^{1} (Q) \cap H_{0; 0,}^{1, 1 / 2} (Q)$ and the conforming test space $Q_{h}^{1} (Q) \cap H_{0;, 0}^{1, 1 / 2} (Q)$ is not stable, see [48, Section 3.3]. A possible way out is the modified Hilbert transformation $H_{T}$ defined by

(H_{T} u) (x, t) := \infty \sum i = 1 \infty \sum k = 0 u_{i, k} cos ((\frac{π}{2} + k π) \frac{t}{T}) ϕ_{i} (x), (x, t) \in Q,

where the given function $u \in L^{2} (Q)$ is represented by

u (x, t) = \infty \sum i = 1 \infty \sum k = 0 u_{i, k} sin ((\frac{π}{2} + k π) \frac{t}{T}) ϕ_{i} (x), (x, t) \in Q,

(6)

with the eigenfunctions

ϕ_{i} \in H_{0}^{1} (Ω)

and eigenvalues

μ_{i} \in R

, satisfying

- Δ ϕ_{i} = μ_{i} ϕ_{i} in Ω, ϕ_{i} = 0 on \partial Ω, {∥ ϕ_{i} ∥}_{L^{2} (Ω)} = 1, i \in N .

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

H_{T} : H_{0; 0,}^{1, 1 / 2} (Q) \to H_{0;, 0}^{1, 1 / 2} (Q)

is norm preserving, bijective and fulfills the coercivity property

{⟨ \partial_{t} u, H_{T} v ⟩}_{Q} = \frac{1}{2} \infty \sum i = 1 \infty \sum k = 0 (\frac{π}{2} + k π) u_{i, k} \cdot v_{i, k} =: {⟨ u, v ⟩}_{H_{0,}^{1 / 2} (0, T; L^{2} (Ω)), F}

(7)

for functions $u, v \in H_{0,}^{1 / 2} (0, T; L^{2} (Ω))$ with expansion coefficients $u_{i, k}, v_{i, k}$ as in (6). Note that the norm induced by the inner product ${⟨ \cdot, \cdot ⟩}_{H_{0,}^{1 / 2} (0, T; L^{2} (Ω)), F}$ is equivalent to the norm ${∥ \cdot ∥}_{H_{0,}^{1 / 2} (0, T; L^{2} (Ω))} .$ Moreover, the relations

	$\forall v \in L^{2} (Q) :$	${⟨ v, H_{T} v ⟩}_{L^{2} (Q)} \geq 0,$
	$\forall s > 0 : \forall v \in H_{0,}^{s} (0, T; L^{2} (Ω)), v \neq 0 :$	${⟨ v, H_{T} v ⟩}_{L^{2} (Q)} > 0$		(8)

hold true. With the modified Hilbert transformation $H_{T}$ , the variational formulation (3) is equivalent to find $u \in H_{0; 0,}^{1, 1 / 2} (Q)$ such that

\forall v \in H_{0; 0,}^{1, 1 / 2} (Q) : a (u, H_{T} v) = {⟨ f, H_{T} v ⟩}_{Q} .

(9)

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

\forall u \in H_{0; 0,}^{1, 1 / 2} (Q) : c ∥ u ∥_{H_{0; 0,}^{1, 1 / 2} (Q)} \leq sup 0 \neq v \in H_{0; 0,}^{1, 1 / 2} (Q) \frac{| a (u, H_{T} v) |}{∥ v ∥_{H_{0; 0,}^{1, 1 / 2} (Q)}}

with a constant $c > 0.$ When using some conforming space-time finite element space $V_{h} \subset H_{0; 0,}^{1, 1 / 2} (Q),$ the Galerkin variational formulation of (9) is to find $u_{h} \in V_{h}$ such that

(10)

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

Theorem 2.2.

Let $V_{h} \subset H_{0; 0,}^{1, 1 / 2} (Q)$ be a conforming space-time finite element space and let $f \in [H_{0;, 0}^{1, 1 / 2} (Q)]^{'}$ be a given right-hand side. Then, a unique solution $u_{h} \in V_{h}$ of the Galerkin variational formulation (10) exists. If, in addition, the right-hand side fulfills $f \in [H_{, 0}^{1 / 2} (0, T; L^{2} (Ω))]^{'} \subset [H_{0;, 0}^{1, 1 / 2} (Q)]^{'},$ then the stability estimate

∥ u_{h} ∥_{H_{0,}^{1 / 2} (0, T; L^{2} (Ω))} \leq c ∥ f ∥_{[H_{, 0}^{1 / 2} (0, T; L^{2} (Ω))]^{'}}

is true with a constant $c > 0.$

Theorem 2.2 states that, under the assumption $f \in [H_{, 0}^{1 / 2} (0, T; L^{2} (Ω))]^{'}$ , any conforming space-time finite element space $V_{h} \subset H_{0; 0,}^{1, 1 / 2} (Q)$ 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

V_{h} = Q_{h}^{p} (Q) \cap H_{0; 0,}^{1, 1 / 2} (Q)

from (5), the Galerkin variational formulation (10) to find $u_{h} \in Q_{h}^{p} (Q) \cap H_{0; 0,}^{1, 1 / 2} (Q)$ such that

\forall v_{h} \in Q_{h}^{p} (Q) \cap H_{0; 0,}^{1, 1 / 2} (Q) : a (u_{h}, H_{T} v_{h}) = {⟨ f, H_{T} v_{h} ⟩}_{Q}

fulfills the space-time error estimates

${∥ u - u_{h} ∥}_{H_{0,}^{1 / 2} (0, T; L^{2} (Ω))}$	$\leq c h^{p + 1 / 2},$	(11)
${∥ u - u_{h} ∥}_{L^{2} (Q)}$	$\leq c h^{p + 1},$	(12)
${\| u - u_{h} \|}_{H^{1} (Q)}$	$\leq c h^{p}$	(13)

with $h = max {h_{t}, h_{x}}$ and with a constant $c > 0$ for a sufficiently smooth solution $u \in H_{0; 0,}^{1, 1 / 2} (Q)$ of (3) and a sufficiently regular boundary $\partial Ω,$ where for the $H^{1} (Q)$ error estimate (13), the sequence $(T_{ν})_{ν \in N}$ of decompositions of $Ω$ is additionally assumed to be globally quasi-uniform, see [41, 48] for details.

In the remainder of this work, we consider $p = 1$ , i.e. the tensor-product space of piecewise linear, continuous functions $V_{h} = Q_{h}^{1} (Q) \cap H_{0; 0,}^{1, 1 / 2} (Q)$ , where analogous results hold true for an arbitrary polynomial degree $p > 1.$

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

d o f = N_{t} \cdot M_{x} .

For an easier implementation, we approximate the right-hand side $f \in L^{2} (Q)$ by

f \approx Q_{h}^{0} f = N_{x} \sum j = 1 N_{t} \sum ℓ = 1 f_{j, ℓ} ψ_{j}^{0} φ_{ℓ}^{0} \in S_{h_{x}}^{0} (Ω) \otimes S_{h_{t}}^{0} (0, T)

(14)

with coefficients $f_{j, ℓ} \in R$ , where $Q_{h}^{0} : L^{2} (Q) \to S_{h_{x}}^{0} (Ω) \otimes S_{h_{t}}^{0} (0, T)$ is the $L^{2} (Q)$ projection on the piecewise constant functions $S_{h_{x}}^{0} (Ω) \otimes S_{h_{t}}^{0} (0, T)$ with $S_{h_{x}}^{0} (Ω) = span {ψ_{j}^{0}}_{j = 1}^{N_{x}}$ and $S_{h_{t}}^{0} (0, T) = span {φ_{ℓ}^{0}}_{ℓ = 1}^{N_{t}}$ . So, we consider the perturbed variational formulation to find ${˜ u}_{h} \in Q_{h}^{1} (Q) \cap H_{0; 0,}^{1, 1 / 2} (Q)$ such that

\forall v_{h} \in Q_{h}^{1} (Q) \cap H_{0; 0,}^{1, 1 / 2} (Q) : a ({˜ u}_{h}, H_{T} v_{h}) = {⟨ Q_{h}^{0} f, H_{T} v_{h} ⟩}_{L^{2} (Q)} .

(15)

Note that, for piecewise linear functions, i.e. $p = 1$ , the space-time error estimates (11), (12), (13) are not spoilt. Note additionally that, for $p > 1$ , a projection on polynomials of degree $p - 1$ should be used instead of $Q_{h}^{0}$ 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

K_{h} u - - = {˜ F - -}^{H_{T}}

(16)

with the system matrix

K_{h} = A_{h_{t}}^{H_{T}} \otimes M_{h_{x}} + M_{h_{t}}^{H_{T}} \otimes A_{h_{x}} \in R^{N_{t} \cdot M_{x} \times N_{t} \cdot M_{x}},

where $M_{h_{x}} \in R^{M_{x} \times M_{x}}$ and $A_{h_{x}} \in R^{M_{x} \times M_{x}}$ denote spatial mass and stiffness matrices given by

M_{h_{x}} [i, j] = ⟨ ψ_{j}^{1}, ψ_{i}^{1} ⟩_{L^{2} (Ω)}, A_{h_{x}} [i, j] = ⟨ \nabla_{x} ψ_{j}^{1}, \nabla_{x} ψ_{i}^{1} ⟩_{L^{2} (Ω)}, i, j = 1, \dots, M_{x},

(17)

and $M_{h_{t}}^{H_{T}} \in R^{N_{t} \times N_{t}}$ and $A_{h_{t}}^{H_{T}} \in R^{N_{t} \times N_{t}}$ are defined by

for $ℓ, k = 1, \dots, N_{t}$ . Note that the matrix $A_{h_{t}}^{H_{T}}$ is dense, symmetric and positive definite, see (7), whereas the matrix $M_{h_{t}}^{H_{T}}$ is dense, nonsymmetric and positive definite, see (8

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

16) is given by

{˜ F - -}^{H_{T}} := {({˜ f - -}_{1}, \dots, {˜ f - -}_{M_{x}})}_{1}^{⊤} \in R^{N_{t} \cdot M_{x}}

with the vectors ${˜ f - -}_{i} \in R^{N_{t}}$ , $i = 1, \dots, M_{x}$ , where, with the help of (14),

{˜ f - -}_{i} [k] := {⟨ Q_{h}^{0} f, ψ_{i}^{1} H_{T} φ_{k}^{1} ⟩}_{L^{2} (Q)} = N_{x} \sum j = 1 N_{t} \sum ℓ = 1 f_{j, ℓ} {⟨ ψ_{j}^{0}, ψ_{i}^{1} ⟩}_{L^{2} (Ω)} {⟨ φ_{ℓ}^{0}, H_{T} φ_{k}^{1} ⟩}_{L^{2} (0, T)},

$k = 1, \dots, N_{t} .$ To assemble the vector of the right-hand side in (16), the relation

{˜ f - -}_{i} [k] := ˜ F [i, k], i = 1, \dots, M_{x}, k = 1, \dots, N_{t},

holds true with $˜ F := M_{h_{x}}^{1, 0} F (C_{h_{t}}^{H_{T}})^{⊤} \in R^{M_{x} \times N_{t}},$ where

$M_{h_{x}}^{1, 0} [i, j]$	$:= {⟨ ψ_{j}^{0}, ψ_{i}^{1} ⟩}_{L^{2} (Ω)},$	$i = 1, \dots, M_{x}, j = 1, \dots, N_{x},$
$F [j, ℓ]$	$:= f_{j, ℓ},$	$j = 1, \dots, N_{x}, ℓ = 1, \dots, N_{t},$
$C_{h_{t}}^{H_{T}} [k, ℓ]$	$:= {⟨ φ_{ℓ}^{0}, H_{T} φ_{k}^{1} ⟩}_{L^{2} (0, T)},$	$k = 1, \dots, N_{t}, ℓ = 1, \dots, N_{t} .$

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 $A, C \in C^{N_{A} \times N_{A}}$ , $B, D \in C^{N_{B} \times N_{B}}$ and $X \in C^{N_{A} \times N_{B}}$ be given matrices for $N_{A}, N_{B} \in N .$ The Kronecker product is defined as the matrix

A \otimes B := ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} A [1, 1] B & A [1, 2] B & \dots & A [1, N_{A}] B A [2, 1] B & A [2, 2] B & \dots & A [2, N_{A}] B ⋮ & ⋮ & ⋱ & ⋮ A [N_{A}, 1] B & A [N_{A}, 2] B & \dots & A [N_{A}, N_{A}] B \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ \in C^{N_{A} \cdot N_{B} \times N_{A} \cdot N_{B}} .

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

v e c (X) := (X [1, 1], X [2, 1], \dots, X [N_{A}, 1], X [1, 2], \dots, X [N_{A}, N_{B}])^{⊤} \in C^{N_{A} \cdot N_{B} \times 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 \otimes B)^{*} = A^{*} \otimes B^{*} and (A \otimes B)^{⊤} = A^{⊤} \otimes B^{⊤} .$
For regular matrices $A, B$ , we have

$(A \otimes B)^{- 1} = A^{- 1} \otimes B^{- 1} .$
The mixed-product property

$(A \otimes B) (C \otimes D) = (A C) \otimes (B D)$

is valid.
It holds true that

$v e c (A X B) = (B^{⊤} \otimes A) v e c (X),$ (18)

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

For a given vector $v - \in C^{N_{A} \cdot N_{B}}$ , define the matrix

V := (\begin{matrix} {v -}_{1} & {v -}_{2} & \dots & {v -}_{N_{B}} \end{matrix}) \in C^{N_{A} \times N_{B}}

with ${v -}_{i} \in C^{N_{A}}$ given by ${v -}_{i} [k] = v - [(i - 1) N_{A} + k]$ for $k = 1, \dots, N_{A},$ $i = 1, \dots, N_{B}$ , i.e. $v e c (V) = v - .$ Then, the equality (18) yields

(B^{⊤} \otimes A) v - = (B^{⊤} \otimes A) v e c (V) = v e c (A V B) .

(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 $K_{h_{x}} z - = g -$ , see, e.g., [7, 10, 17, 22, 27, 30, 34]. Here, $K_{h_{x}} \in C^{n \times n}$ is a sparse matrix coming from finite element/difference discretizations of a physical domain $Ω \subset R^{d}$ , $d = 1, 2, 3$ , like the spatial mass or stiffness matrix (17). Sparse direct solvers exploit the sparsity pattern of the system matrix $K_{h_{x}}$ , 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:

Ordering Step, e.g., minimum degree or nested dissection methods,
Symbolic Factorization Step,
Numerical Factorization Step,
Solving Step.

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

$d$	Ordering and Factorization Steps	Solving Step	Memory
1	$O (n)$	$O (n)$	$O (n)$
2	$O (n^{3 / 2})$	$O (n ln n)$	$O (n ln n)$
3	$O (n^{2})$	$O (n^{4 / 3})$	$O (n^{4 / 3})$

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

(A_{h_{t}}^{H_{T}} \otimes M_{h_{x}} + M_{h_{t}}^{H_{T}} \otimes A_{h_{x}}) u - - = {˜ F - -}^{H_{T}}

(20)

given in (16) with the symmetric, positive definite matrices $M_{h_{x}} \in R^{M_{x} \times M_{x}}$ , $A_{h_{x}} \in R^{M_{x} \times M_{x}}$ , $A_{h_{t}}^{H_{T}} \in R^{N_{t} \times N_{t}}$ and the nonsymmetric, positive definite matrix $M_{h_{t}}^{H_{T}} \in R^{N_{t} \times N_{t}}$ . 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 $(M_{h_{t}}^{H_{T}}, A_{h_{t}}^{H_{T}})$ is decomposed in the form

(A_{h_{t}}^{H_{T}})^{- 1} M_{h_{t}}^{H_{T}} = X_{t} Z_{t} X_{t}^{- 1}

with real, regular matrices $X_{t}, Z_{t} \in R^{N_{t} \times N_{t}}$ , where $Z_{t}$ is an upper (quasi-)triangular matrix, or complex, regular matrices $X_{t}, Z_{t} \in C^{N_{t} \times N_{t}}$ , where $Z_{t}$ is an upper triangular or diagonal matrix. Defining

Y_{t} := (A_{h_{t}}^{H_{T}} X_{t})^{- 1}

gives the representations

A_{h_{t}}^{H_{T}} = Y_{t}^{- 1} X_{t}^{- 1} and M_{h_{t}}^{H_{T}} = A_{h_{t}}^{H_{T}}    = Y_{t}^{- 1} X_{t}^{- 1} X_{t} Z_{t} X_{t}^{- 1} = Y_{t}^{- 1} Z_{t} X_{t}^{- 1} .

Hence, the global linear system is equivalent to solving

(Y_{t}^{- 1} \otimes I_{M_{x}}) (I_{N_{t}} \otimes M_{h_{x}} + Z_{t} \otimes A_{h_{x}}) (X_{t}^{- 1} \otimes I_{M_{x}}) u - - = {˜ F - -}^{H_{T}}

with the identity matrices $I_{N_{t}} \in R^{N_{t} \times N_{t}}$ and $I_{M_{x}} \in R^{M_{x} \times M_{x}}$ . Thus, the solution of (20) is given by

u - - = (X_{t} \otimes I_{M_{x}}) (I_{N_{t}} \otimes M_{h_{x}} + Z_{t} \otimes A_{h_{x}})^{- 1} (Y_{t} \otimes I_{M_{x}}) {˜ F - -}^{H_{T}} .

(21)

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

g - := ({g -}_{1}, {g -}_{2}, \dots, {g -}_{N_{t}})^{⊤} := (Y_{t} \otimes I_{M_{x}}) {˜ F - -}^{H_{T}} = v e c (^F (A_{h_{t}}^{H_{T}})^{- 1} X_{t}^{- ⊤}) \in C^{N_{t} \cdot M_{x}}

(22)

with a matrix $^F \in R^{M_{x} \times N_{t}}$ corresponding to the relation (19), satisfying $v e c (^F) = {˜ F - -}^{H_{T}},$ where ${g -}_{i} \in C^{M_{x}} .$ The second step in (21) is to solve the linear system

(I_{N_{t}} \otimes M_{h_{x}} + Z_{t} \otimes A_{h}

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

Authors

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

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

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

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

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

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

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

1 Introduction

2 Space-Time Method in Anisotropic Sobolev Spaces

Theorem 2.1.

Theorem 2.2.

3 Preliminaries for the Space-Time Solvers

3.1 Kronecker Product

3.2 Sparse Direct Solver

4 Space-Time Solvers

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

Authors

Related Research

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

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

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

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

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

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

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

This week in AI

1 Introduction

2 Space-Time Method in Anisotropic Sobolev Spaces

Theorem 2.1.

Theorem 2.2.

3 Preliminaries for the Space-Time Solvers

3.1 Kronecker Product

3.2 Sparse Direct Solver

4 Space-Time Solvers