# Parallelization in time by diagonalization

Andrea Bressan<EMAIL_ADDRESS>Alen Kushova
<EMAIL_ADDRESS>Gabriele Loli<EMAIL_ADDRESS>Monica Montardini<EMAIL_ADDRESS>Giancarlo Sangalli
<EMAIL_ADDRESS>Mattia Tani<EMAIL_ADDRESS>Dipartimento di
Matematica “F. Casorati”, Università di Pavia, Via A. Ferrata, 5, 27100 Pavia,
Italy.  Istituto di Matematica Applicata e Tecnologie Informatiche, “E.
Magenes” del CNR, Via A. Ferrata, 1, 27100 Pavia, Italy.

###### Abstract

This is a review of preconditioning techniques based on fast-diagonalization
methods for space-time isogeometric discretization of the heat equation. Three
formulation are considered: the Galerkin approach, a discrete least-square and
a continuous least square. For each formulation the heat differential operator
is written as a sum of terms that are kronecker products of uni-variate
operators. These are used to speed-up the application of the operator in
iterative solvers and to construct a suitable preconditioner. Contrary to the
fast-diagonalization technique for the Laplace equation where all uni-variate
operators acting on the same direction can be simultaneously diagonalized in
the case of the heat equation this is not possible. Luckily this can be done
up to an additional term that has low rank allowing for the utilization of
arrow-head like factorization or inversion by Sherman-Morrison formula. The
proposed preconditioners work extremely well on the parametric domain and,
when the domain is parametrized or when the equation coefficients are not
constant, they can be adapted and retain good performance characteristics.

###### keywords:

Isogeometric Analysis , splines , heat equation , space-time Galerkin
formulation , Fast Diagonalization.

##  1 Introduction

This work is a review of preconditioners for isogeometric space-time
discretizations of the heat equation.

Isogeometric Analysis (IgA), introduced in the seminal paper [1] (see also the
book [2]), is an evolution of the classical finite element methods. IgA uses
spline functions, or their generalizations, both to represent the
computational domain and to approximate the solution of the partial
differential equation that models the problem of interest. This is meant to
simplify the interoperability between computer aided design and numerical
simulations. IgA also benefits from the approximation properties of splines,
whose high-continuity yields higher accuracy when compared to $C^{0}$
piecewise polynomials, see e.g., [3, 4, 5].

Space-time finite element methods originated in the papers [6, 7, 8], where
standard finite elements are assigned an extra dimension for the time and,
typically, adopt a discontinuous approximation in time, since this produces a
time marching algorithm with a traditional step-by-step format (see e.g. [9]).

One of the first work concerning space-time isogeometric discretization is
[10], in which a stabilized variational formulation produces a discrete
bilinear form that is elliptic with respect to a discrete energy norm. The
resulting linear system is then solved through a standard parallel AMG
preconditioned GMRES solver. Other papers in literature propose isogeometric
space-time Galerkin methods, favoring a step-by-step structure in time. In
[11], the same variational formulation of [10] is used in combination with a
space-time domain decomposition into space-time slabs that are sequentially
coupled in time by a stabilized discontinuous Galerkin method. In [12], two
different methods, called ST-C-SPT and ST-C-DCT, are outlined.

The first one, analysed in [13], is a way to project a previously computed
solution, possibly discontinuous, into isogeometric spaces in order to get a
more regular solution and to save memory for its storage.

In the ST-C-DCT method, the solution with continuous temporal representation
is computed sequentially from the space-time variational formulation
associated with each slab.

Related multigrid solvers have been proposed in [14, 15] and low-rank
approximations in [16]. In [17] the authors consider $C^{0}$ coupling between
the space-time slabs with a suitable stabilized formulation that also yields a
sequential scheme. Finally, the interest in space-time isogeometric analysis
for complex real-world simulations is attested by the recent papers [18, 19,
20], where, again, a sequential (discontinuous) approximation in time is
adopted.

In this review paper we focus on the heat equation and on its space-time
isogeometric discretizations, that allows smooth approximation in both space
and time. This originated in [21], focusing on a $L^{2}$ least-squares
formulation, while in [22], the author focused on the plain Galerkin space-
time formulation, whose well-posedness has been studied, for finite element
discretizations and for the heat equation, in the recent papers [23] and [24].
Here we consider both this formulations, focusing on the study of stable
preconditioning strategies. Indeed, when adopting smooth approximation in
space and in time, the major issue is its computational cost and the key
ingredient is an efficient solver for the linear system, which is global in
time.

The preconditioners proposed by the authors in [21] and [22], exploit the
Kronecker structure of the arising linear systems. For the plain Galerkin
formulation, and assuming that the spatial domain does not change with time,
the linear system has the structure

$\mathbf{A}:=\gamma\mathbf{W}_{t}\otimes\mathbf{M}_{s}+\nu\mathbf{M}_{t}\otimes\mathbf{L}_{s},$
(1.1)

where $\mathbf{W}_{t}$ is given by the discretization of the time derivative,
$\mathbf{L}_{s}$ is given by the discretization of the Laplacian in the
spatial variables, $\mathbf{M}_{t}$ and $\mathbf{M}_{s}$ are “mass matrices”
in time and space, respectively, and $\gamma,\nu>0$ are constants of the
problem. The construction of a preconditioner for (1.1) is based on a
generalization of the classical Fast Diagonalization (FD) method [25]. Indeed,
the FD method cannot be directly applied to (1.1), as this would require to
compute the eigendecomposition of the pencil $(\mathbf{W}_{t},\mathbf{M}_{t})$
which is numerically unstable. We circumvent this difficulty by introducing an
ad-hoc factorization of the time matrices which allows to design a solver
conceptually similar to the FD method.

For the Least Squares formulation, instead, the linear system has the
structure

$\mathbf{B}:=\gamma^{2}\mathbf{L}_{t}\otimes\mathbf{M}_{s}+\nu^{2}\mathbf{M}_{t}\otimes\mathbf{J}_{s}+\gamma\nu\mathbf{R}_{t}\otimes\mathbf{L}_{s},$
(1.2)

where $\mathbf{L}_{t}$ is the stiffness matrix w.r.t. the time derivative,
$\mathbf{L}_{s}$ is given by the discretization of the Laplacian in the
spatial variables, $\mathbf{R}_{t}$ is a rank 1 matrix associated to the final
time, and $\mathbf{J}_{s}$ is given by the discretization of the Bi-Laplacian
in the spatial variables. Again, $\mathbf{M}_{t}$ and $\mathbf{M}_{s}$ are
“mass matrices” in time and space, respectively, and $\gamma,\nu>0$ are
constants of the problem. Thus, the problem becomes elliptic and a
preconditioner for the linear system associated to (1.2) can be easily
designed as in [26], again based on the FD method.

Alongside the preconditioners mentioned above, we propose a third approach
that is suitable for both introduced formulations. In the case of the plain
Galerkin formulation, we can consider the discrete least-square form
associated to the system (1.1), that is a linear system with matrix

$\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{A}:=\gamma^{2}\mathbf{W}_{t}^{T}\mathbf{M}_{t}^{-1}\mathbf{W}_{t}\otimes\mathbf{M}_{s}+\nu^{2}\mathbf{M}_{t}\otimes\mathbf{J}_{s}\mathbf{M}_{s}^{-1}\mathbf{J}_{s}+\gamma\nu\mathbf{R}_{t}\otimes\mathbf{L}_{s}.$

This has a more similar structure to (1.2). In particular both discrete
problems are now elliptic, and the third term appearing is the Kronecker
product between a rank 1 perturbation and the stiffness matrix in space. Thus
we designed two new preconditioners, one for the discrete least square
formulation above and one for (1.2), both of them relying on FD method and
Sherman-Morrison formula.

The computational cost of the setup of all mentioned preconditioners is
$O(N_{dof})$ FLOating-Point operations (FLOPs), while the application cost is
$O(N_{dof}^{1+1/d})$ FLOPs, where $d$ is the number of spatial dimensions and
$N_{dof}$ denotes the total number of degrees-of-freedom (assuming, for
simplicity, to have the same number of degrees-of-freedom in time and in each
spatial direction). Our numerical benchmarks show that the computing time
(serial and single-core execution) is close to optimality, that is,
proportional to $N_{dof}$. The preconditioner are also robust with respect to
the polynomial degree and number of elements. Furthermore, our approaches are
optimal in terms of memory requirement: denoting by $N_{s}$ the total number
of degrees-of-freedom in space, the storage cost is $O(p^{d}N_{s}+N_{dof})$.
We also remark that global space-time methods in principle facilitate the full
parallelization of the solver, see [27, 28, 29].

The outline of the paper is as follows. In Section 2 we present the basics of
B-splines based IgA and the main properties of the Kronecker product
operation. The plain Galerkin formulation and its isogeometric discretization
are introduced in Section 3, while in Section 4 we recall the preconditioner
introduced in [22] and we discuss its application. In Section 5 we introduce
the new preconditioner for the discrete least square form plain Galerkin
discretization. The Least Squares formulation and its isogeometric
discretization are introduced in Section 6, while in Section 7 we recall the
preconditioner introduced in [21] and we discuss its application. In Section 8
we introduce the new preconditioner for the least square formulation. Section
9 is devoted to the computational costs of the proposed preconditioners and to
memory requirements. We present the numerical results assessing the
performance of the proposed preconditioner in Section 10. Finally, in the last
section we draw some conclusions and we highlight some future research
directions.

##  2 Preliminaries

###  2.1 B-Splines

Given $m$ and $p$ two positive integers, a knot vector in $[0,1]$ is a
sequence of non-decreasing points
$\Xi:=\left\\{0=\xi_{1}\leq\dots\leq\xi_{m+p+1}=1\right\\}$. We consider open
knot vectors, i.e. we set $\xi_{1}=\dots=\xi_{p+1}=0$ and
$\xi_{m}=\dots=\xi_{m+p+1}=1$. Then, according to Cox-De Boor recursion
formulas (see [30]), univariate B-splines
$\widehat{b}_{i,p}:(0,1)\rightarrow\mathbb{R}$ are piecewise polynomials
defined for $i=1,\dots,m$ as
for $p=0$:

$\displaystyle\widehat{b}_{i,0}(\eta)=\begin{cases}1&{\textrm{if
}}\xi_{i}\leq\eta<\xi_{i+1},\\\ 0&\textrm{otherwise,}\end{cases}$

for $p\geq 1$:

$\displaystyle\widehat{b}_{i,p}(\eta)=\\!\begin{cases}\dfrac{\eta-\xi_{i}}{\xi_{i+p}-\xi_{i}}\widehat{b}_{i,p-1}(\eta)+\dfrac{\xi_{i+p+1}-\eta}{\xi_{i+p+1}-\xi_{i+1}}\widehat{b}_{i+1,p-1}(\eta)&{\textrm{if
}}\xi_{i}\leq\eta<\xi_{i+p+1},\\\\[8.0pt] 0&\textrm{otherwise,}\end{cases}$

where we adopt the convention $0/0=0$. The univariate spline space is defined
as

$\widehat{\mathcal{S}}_{h}^{p}:=\mathrm{span}\\{\widehat{b}_{i,p}\\}_{i=1}^{m},$

where $h$ denotes the mesh-size, i.e. $h:=\max\\{|\xi_{i+1}-\xi_{i}|\ |\
i=1,\dots,m+p\\}$. The interior knot multiplicity influences the smoothness of
the B-splines at the knots (see [30]). For more details on B-splines
properties and their use in IgA we refer to [2].

Multivariate B-splines are defined as tensor product of univariate B-splines.
We consider functions that depend on $d$ spatial variables and the time
variable. Therefore, given positive integers $m_{l},p_{l}$ for $l=1,\dots,d$
and $m_{t},p_{t}$, we introduce $d+1$ univariate knot vectors
$\Xi_{l}:=\left\\{\xi_{l,1}\leq\dots\leq\xi_{l,m_{l}+p_{l}+1}\right\\}$ for
$l=1,\ldots,d$ and
$\Xi_{t}:=\left\\{\xi_{t,1}\leq\dots\leq\xi_{t,m_{t}+p_{t}+1}\right\\}$. Let
$h_{l}$ be the mesh-size associated to the knot vector $\Xi_{l}$ for
$l=1,\dots,d$, let $h_{s}:=\max\\{h_{l}\ |\ l=1,\dots,d\\}$ be the maximal
mesh-size in all spatial knot vectors and let $h_{t}$ be the mesh-size of the
time knot vector. Let also $\boldsymbol{p}$ be the vector that contains the
degree indexes, i.e. $\boldsymbol{p}:=(\boldsymbol{p}_{s},p_{t})$, where
$\boldsymbol{p}_{s}:=(p_{1},\dots,p_{d})$. For simplicity, we assume to have
the same polynomial degree in all spatial directions, i.e., with abuse of
notations, we set $p_{1}=\dots=p_{d}=:p_{s}$, but the general case is similar.

We assume that the following quasi-uniformity of the knot vectors holds.

######  Assumption 1.

There exists $0<\alpha\leq 1$, independent of $h_{s}$ and $h_{t}$, such that
each non-empty knot span $(\xi_{l,i},\xi_{l,i+1})$ of $\Xi_{l}$ fulfils
$\alpha h_{s}\leq\xi_{l,i+1}-\xi_{l,i}\leq h_{s}$ for $l=1,\dots,d$ and each
non-empty knot-span $(\xi_{t,i},\xi_{t,i+1})$ of $\Xi_{t}$ fulfils $\alpha
h_{t}\leq\xi_{t,i+1}-\xi_{t,i}\leq h_{t}$.

The multivariate B-splines are defined as

$\widehat{B}_{\boldsymbol{i},\boldsymbol{p}}(\boldsymbol{\eta},\tau):=\widehat{B}_{\boldsymbol{i_{s}},\boldsymbol{p}_{s}}(\boldsymbol{\eta})\widehat{b}_{i_{t},p_{t}}(\tau),$

where

$\widehat{B}_{\boldsymbol{i_{s}},\boldsymbol{p}_{s}}(\boldsymbol{\eta}):=\widehat{b}_{i_{1},p_{s}}(\eta_{1})\ldots\widehat{b}_{i_{d},p_{s}}(\eta_{d}),$
(2.1)

$\boldsymbol{i_{s}}:=(i_{1},\dots,i_{d})$,
$\boldsymbol{i}:=(\boldsymbol{i_{s}},i_{t})$ and
$\boldsymbol{\eta}=(\eta_{1},\ldots,\eta_{d})$. The corresponding spline space
is defined as

$\widehat{\mathcal{S}}^{\boldsymbol{p}}_{{h}}:=\mathrm{span}\left\\{\widehat{B}_{\boldsymbol{i},\boldsymbol{p}}\
\middle|\ i_{l}=1,\dots,m_{l}\text{ for
}l=1,\dots,d;i_{t}=1,\dots,m_{t}\right\\},$

where $h:=\max\\{h_{s},h_{t}\\}$. We have that
$\widehat{\mathcal{S}}^{\boldsymbol{p}}_{{h}}=\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{{h}_{s}}\otimes\widehat{\mathcal{S}}^{p_{t}}_{h_{t}},$
where

$\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{h_{s}}:=\mathrm{span}\left\\{\widehat{B}_{\boldsymbol{i_{s}},\boldsymbol{p}_{s}}\
\middle|\ i_{l}=1,\dots,m_{l};l=1,\dots,d\right\\}$

is the space of tensor-product splines on $\widehat{\Omega}:=(0,1)^{d}$.

######  Assumption 2.

We assume that $p_{t},p_{s}\geq 1$ and that
$\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{h_{s}}\subset
C^{0}(\widehat{\Omega})$ and $\widehat{\mathcal{S}}^{{p}_{t}}_{h_{t}}\subset
C^{0}\left((0,1)\right)$ .

###  2.2 Isogeometric spaces

The space-time computational domain that we consider is $\Omega\times(0,T)$,
where $\Omega\subset\mathbb{R}^{d}$ and $T>0$ is the final time. We make the
following assumptions.

######  Assumption 3.

We assume that $\Omega$ is parametrized by
$\boldsymbol{F}:\widehat{\Omega}\rightarrow{\Omega}$, with
$\boldsymbol{F}\in{\left[\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{{h}_{s}}\right]^{d}}$
on the closure of $\widehat{\Omega}$. Moreover, we assume that
$\boldsymbol{F}^{-1}$ has piecewise bounded derivatives of any order.

We define
$\boldsymbol{x}=(x_{1},\dots,x_{d}):=\boldsymbol{F}(\boldsymbol{\eta})$ and
$t:=T\tau$. Then the space-time domain is given by the parametrization
$\boldsymbol{G}:\widehat{\Omega}\times(0,1)\rightarrow\Omega\times(0,T)$, such
that
$\boldsymbol{G}(\boldsymbol{\eta},\tau):=(\boldsymbol{F}(\boldsymbol{\eta}),T\tau)=(\boldsymbol{x},t).$

We introduce the spline space with initial and boundary conditions, in
parametric coordinates, as

$\widehat{\mathcal{X}}_{h}:=\left\\{\widehat{v}_{h}\in\widehat{\mathcal{S}}^{\boldsymbol{p}}_{h}\
\middle|\ \widehat{v}_{h}=0\text{ on
}\partial\widehat{\Omega}\times(0,1)\text{ and }\widehat{v}_{h}=0\text{ on
}\widehat{\Omega}\times\\{0\\}\right\\}.$

We also have that
$\widehat{\mathcal{X}}_{h}=\widehat{\mathcal{X}}_{s,h_{s}}\otimes\widehat{\mathcal{X}}_{t,h_{t}}$,
where

$\displaystyle\widehat{\mathcal{X}}_{s,h_{s}}$
$\displaystyle:=\left\\{\widehat{w}_{h}\in\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{h_{s}}\
\middle|\ \widehat{w}_{h}=0\text{ on }\partial\widehat{\Omega}\right\\}\ $
$\displaystyle\ =\
\text{span}\left\\{\widehat{b}_{i_{1},p_{s}}\dots\widehat{b}_{i_{d},p_{s}}\
\middle|\ i_{l}=2,\dots,m_{l}-1;\ l=1,\dots,d\ \right\\},$
$\displaystyle\widehat{\mathcal{X}}_{t,h_{t}}$
$\displaystyle:=\left\\{\widehat{w}_{h}\in\widehat{\mathcal{S}}^{p_{t}}_{h_{t}}\
\middle|\ \widehat{w}_{h}(0)=0\right\\}\ =\
\text{span}\left\\{\widehat{b}_{i_{t},p_{t}}\ \middle|\ i_{t}=2,\dots,m_{t}\
\right\\}.$

By introducing a colexicographical reordering of the basis functions, we can
write

$\displaystyle\widehat{\mathcal{X}}_{s,h_{s}}$ $\displaystyle=\
\text{span}\left\\{\widehat{b}_{i_{1},p_{s}}\dots\widehat{b}_{i_{d},p_{s}}\
\middle|\ i_{l}=1,\dots,n_{s,l};\ l=1,\dots,d\ \right\\}$ $\displaystyle\
=\text{span}\left\\{\widehat{B}_{i,\boldsymbol{p}_{s}}\ \middle|\
i=1,\dots,N_{s}\ \right\\},$ $\displaystyle\ \widehat{\mathcal{X}}_{t,h_{t}}$
$\displaystyle=\ \text{span}\left\\{\widehat{b}_{i,p_{t}}\ \middle|\
i=1,\dots,n_{t}\ \right\\},$

and then

$\widehat{\mathcal{X}}_{h}=\text{span}\left\\{\widehat{B}_{{i},\boldsymbol{p}}\
\middle|\ i=1,\dots,N_{dof}\right\\},$ (2.3)

where we defined $n_{s,l}:=m_{l}-2$ for $l=1,\dots,d$,
$N_{s}:=\prod_{l=1}^{d}n_{s,l}$, $n_{t}:=m_{t}-1$ and $N_{dof}:=N_{s}n_{t}$.

Finally, the isogeometric space we consider is the isoparametric push-forward
of (2.3) through the geometric map $\boldsymbol{G}$, i.e.

$\mathcal{X}_{h}:=\text{span}\left\\{B_{i,\boldsymbol{p}}:=\widehat{B}_{i,\boldsymbol{p}}\circ\boldsymbol{G}^{-1}\
\middle|\ i=1,\dots,N_{dof}\right\\}.$ (2.4)

We also have that
$\mathcal{X}_{h}=\mathcal{X}_{s,h_{s}}\otimes\mathcal{X}_{t,h_{t}}$, where

$\mathcal{X}_{s,h_{s}}:=\text{span}\left\\{{B}_{i,\boldsymbol{p}_{s}}:=\widehat{B}_{i,\boldsymbol{p}_{s}}\circ\boldsymbol{F}^{-1}\
\middle|\ i=1,\dots,N_{s}\right\\}$

and

$\mathcal{X}_{t,h_{t}}:=\text{span}\left\\{{b}_{i,p_{t}}:=\widehat{b}_{i,p_{t}}(\cdot/T)\
\middle|\ i=1,\dots,n_{t}\right\\}.$

###  2.3 Kronecker product

The Kronecker product of two matrices $\mathbf{C}\in\mathbb{C}^{n_{1}\times
n_{2}}$ and $\mathbf{D}\in\mathbb{C}^{n_{3}\times n_{4}}$ is defined as

$\mathbf{C}\otimes\mathbf{D}:=\begin{bmatrix}[\mathbf{C}]_{1,1}\mathbf{D}&\dots&[\mathbf{C}]_{1,n_{2}}\mathbf{D}\\\
\vdots&\ddots&\vdots\\\
[\mathbf{C}]_{n_{1},1}\mathbf{D}&\dots&[\mathbf{C}]_{n_{1},n_{2}}\mathbf{D}\end{bmatrix}\in\mathbb{C}^{n_{1}n_{3}\times
n_{2}n_{4}},$

where $[\mathbf{C}]_{i,j}$ denotes the $ij$-th entry of the matrix
$\mathbf{C}$. For extensions and properties of the Kronecker product we refer
to [31]. In particular, when a matrix has a Kronecker product structure, the
matrix-vector product can be efficiently computed. For this purpose, for
$m=1,\dots,d+1$ we introduce the $m$-mode product of a tensor
$\mathfrak{X}\in\mathbb{C}^{n_{1}\times\dots\times n_{d+1}}$ with a matrix
$\mathbf{J}\in\mathbb{C}^{\ell\times n_{m}}$, that we denote by
$\mathfrak{X}\times_{m}\mathbf{J}$. This is a tensor of size
$n_{1}\times\dots\times n_{m-1}\times\ell\times n_{m+1}\times\dots n_{d+1}$,
whose elements are defined as

$\left[\mathfrak{X}\times_{m}\mathbf{J}\right]_{i_{1},\dots,i_{d+1}}:=\sum_{j=1}^{n_{m}}[\mathfrak{X}]_{i_{1},\dots,i_{m-1},j,i_{m+1},\dots,i_{d+1}}[\mathbf{J}]_{i_{m},j}.$

Then, given $\mathbf{J}_{i}\in\mathbb{C}^{\ell_{i}\times n_{i}}$ for
$i=1,\dots,d+1$, it holds

$\left(\mathbf{J}_{d+1}\otimes\dots\otimes\mathbf{J}_{1}\right)\text{vec}\left(\mathfrak{X}\right)=\text{vec}\left(\mathfrak{X}\times_{1}\mathbf{J}_{1}\times_{2}\dots\times_{d+1}\mathbf{J}_{d+1}\right),$
(2.5)

where the vectorization operator “vec” applied to a tensor stacks its entries
into a column vector as

$[\text{vec}(\mathfrak{X})]_{j}=[\mathfrak{X}]_{i_{1},\dots,i_{d+1}}\text{ for
}i_{l}=1,\dots,n_{l}\text{ and for }l=1,\dots,d+1,$

where $j:=i_{1}+\sum_{k=2}^{d+1}\left[(i_{k}-1)\Pi_{l=1}^{k-1}n_{l}\right]$.

##  3 Space-time heat equation

###  3.1 Space-time variational formulation

Our model problem is the heat equation with homogeneous boundary and initial
conditions: we look for a solution $u$ such that

$\left\\{\begin{array}[]{rcllrcl}\gamma\partial_{t}u-\nabla\cdot(\nu\nabla
u)&=&f&\mbox{in }&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!(0,T),\\\\[1.0pt]
u&=&0&\mbox{on
}&\partial\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\![0,T],\\\\[1.0pt]
u&=&0&\mbox{in
}&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!\\{0\\},\end{array}\right.$ (3.1)

where $\Omega\subset\mathbb{R}^{d}$, $T$ is the final time, $\gamma>0$ is the
heat capacity constant and $\nu>0$ is the thermal conductivity constant. We
assume that $f\in L^{2}(0,T;H^{-1}(\Omega))$ and we introduce the Hilbert
spaces

$\mathcal{X}:=\left\\{v\in L^{2}(0,T;H_{0}^{1}(\Omega))\cap
H^{1}(0,T;H^{-1}(\Omega))\mid v(\boldsymbol{x},0)=0\right\\},$
$\mathcal{Y}:=L^{2}(0,T;H_{0}^{1}(\Omega)),$

endowed with the following norms

$\|v\|_{\mathcal{X}}^{2}:=\frac{\gamma^{2}}{{\nu}}\|\partial_{t}v\|^{2}_{L^{2}(0,T;H^{-1}(\Omega))}+\nu\|v\|^{2}_{L^{2}(0,T;H^{1}_{0}(\Omega))}\
\text{ and }\
\|v\|_{\mathcal{Y}}^{2}:=\nu\|v\|^{2}_{L^{2}(0,T;H^{1}_{0}(\Omega))},$

respectively. Then, the variational formulation of (3.1) reads:

$\text{Find }u\in\mathcal{X}\text{ such that
}\mathcal{A}({u},v)=\mathcal{F}(v)\quad\,\forall v\in\mathcal{Y},$ (3.2)

where the bilinear form $\mathcal{A}(\cdot,\cdot)$ and the linear form
$\mathcal{F}(\cdot)$ are defined $\forall w\in\mathcal{X}\text{ and }\forall
v\in\mathcal{Y}$ as

$\mathcal{A}(w,v):=\int_{0}^{T}\int_{\Omega}\left(\gamma\partial_{t}{w}\,v+\nu\nabla
w\cdot\nabla
v\right)\,\mathrm{d}\Omega\,\mathrm{dt}\quad\text{and}\quad\mathcal{F}(v):=\int_{0}^{T}\int_{\Omega}f\,v\,\mathrm{d}\Omega\,\mathrm{dt}.$

The well-posedness of the variational formulation above is a classical result,
see for example [23].

The previous setting can be generalized to non-homogeneous initial and
boundary conditions. For example, suppose that in (3.1) we have the initial
condition $u=u_{0}$ in $\Omega\times\\{0\\}$ with $u_{0}\in L^{2}(\Omega).$
Then, we consider a lifting $\underline{u}_{0}$ of $u_{0}$ such that
$\underline{u}_{0}\in L^{2}(0,T;H_{0}^{1}(\Omega))\cap
H^{1}(0,T;H^{-1}(\Omega))$, see e.g. [32]. Finally, we split the solution $u$
as $u=\underline{u}+\underline{u}_{0}$, where $\underline{u}\in\mathcal{X}$ is
the solution of the following heat equation with homogeneous initial and
boundary conditions:

$\left\\{\begin{array}[]{rcllrcl}\gamma\partial_{t}\underline{u}-\nabla\cdot(\nu\nabla\underline{u})&=&\underline{f}&\mbox{in
}&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!(0,T),\\\\[1.0pt]
\underline{u}&=&0&\mbox{on
}&\partial\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\![0,T],\\\\[1.0pt]
\underline{u}&=&0&\mbox{in
}&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!\\{0\\},\end{array}\right.$

where
$\underline{f}:=f-\gamma\partial_{t}\underline{u}_{0}+\nabla\cdot(\nu\nabla\underline{u}_{0})$.

###  3.2 Space-time Galerkin method

Let $\mathcal{X}_{h}\subset\mathcal{X}$ be the isogeometric space defined in
(2.4). We consider the following Galerkin method for (3.2):

$\text{Find }u_{h}\in\mathcal{X}_{h}\text{ such that
}\mathcal{A}(u_{h},v_{h})=\mathcal{F}(v_{h})\quad\,\forall
v_{h}\in\mathcal{X}_{h}.$ (3.3)

Following [23], let
$N_{h}:L^{2}(0,T;H^{-1}(\Omega))\rightarrow\mathcal{X}_{h}$ be the discrete
Newton potential operator, defined as follows: given $\phi\in
L^{2}(0,T;H^{-1}(\Omega))$ then $N_{h}\phi\in\mathcal{X}_{h}$ fulfills

$\int_{0}^{T}\int_{\Omega}{\nu}\nabla(N_{h}\phi)\cdot\nabla
v_{h}\,\mathrm{d}\Omega\,\mathrm{dt}=\gamma\int_{0}^{T}\int_{\Omega}\phi\
v_{h}\,\mathrm{d}\Omega\,\mathrm{dt}\quad\forall v_{h}\in\mathcal{X}_{h}.$

Thus, we define the norm in $\mathcal{X}_{h}$ as

$\|w\|^{2}_{\mathcal{X}_{h}}:=\nu\|N_{h}(\partial_{t}w)\|^{2}_{L^{2}(0,T;H^{1}_{0}(\Omega))}+{\nu}\|w\|^{2}_{L^{2}(0,T;H^{1}_{0}(\Omega))}.$

The stability and the well-posedness of formulation (3.3) are guaranteed by
[23, Equation (2.7)] and by a straightforward extension to IgA of [23, Theorem
3.1] and [23, Theorem 3.2]. We summarize these results in the following
Proposition 1 and Theorem 1.

######  Proposition 1.

It holds

$\mathcal{A}(w,v)\leq\sqrt{2}\|w\|_{\mathcal{X}}\|v\|_{\mathcal{Y}}\quad\forall
w\in\mathcal{X}\text{ and }\forall v\in\mathcal{Y},$

and

$\|w_{h}\|_{\mathcal{X}_{h}}\leq{2\sqrt{2}}\sup_{v_{h}\in\mathcal{X}_{h}}\frac{\mathcal{A}(w_{h},v_{h})}{\|v_{h}\|_{\mathcal{Y}}}\quad\forall
w_{h}\in\mathcal{X}_{h}.$

######  Theorem 1.

There exists a unique solution $u_{h}\in\mathcal{X}_{h}$ to the discrete
problem (3.3). Moreover, it holds

$\|u-u_{h}\|_{\mathcal{X}_{h}}\leq
5\inf_{w_{h}\in\mathcal{X}_{h}}\|u-w_{h}\|_{\mathcal{X}},$

where $u\in\mathcal{X}$ is the solution of (3.2).

We have then the following a-priori estimate for $h$-refinement.

######  Theorem 2.

Let $q$ be an integer such that $1<q\leq\min\\{p_{s},p_{t}\\}+1$. If
$u\in\mathcal{X}\cap H^{q}\left(\Omega\times(0,T)\right)$ is the solution of
(3.2) and $u_{h}\in\mathcal{X}_{h}$ is the solution of (3.3), then it holds

$\|u-u_{h}\|_{\mathcal{X}_{h}}\leq
C\sqrt{\frac{\gamma^{2}}{\nu}+\nu}\left(h_{t}^{q-1}+h_{s}^{q-1}\right)\|u\|_{H^{q}\left(\Omega\times(0,T)\right)}$
(3.4)

where $C$ is independent of $h_{s},h_{t},\gamma$, $\nu$ and $u$.

The proof of (3.4) is given in [22, Theorem 2], we conclude with the following
remark.

######  Remark 1.

In Theorem 1, the degrees $p_{t}$, $p_{s}$ and the mesh-sizes $h_{t}$, $h_{s}$
play a similar role. This motivates our choice $p_{t}=p_{s}=:p$ and
$h_{t}=h_{s}=:h$ for the numerical tests in Section 10. In this case, and if
the solution $u$ is smooth, (3.4) yields $h$-convergence of order $p$. A
sharper error analysis is possible taking into account a different regularity
of the solution $u$ in space and time, in the line of the anisotropic
estimates of [33].

###  3.3 Discrete system

The linear system associated to (3.3) is

$\mathbf{A}\mathbf{u}=\mathbf{f},$ (3.5)

where
$[\mathbf{A}]_{i,j}=\mathcal{A}(B_{j,\boldsymbol{p}},B_{i,\boldsymbol{p}})$
and $[\mathbf{f}]_{i}=\mathcal{F}(B_{i,\boldsymbol{p}})$. The tensor-product
structure of the isogeometric space (2.4) allows to write the system matrix
$\mathbf{A}$ as sum of Kronecker products of matrices as

$\mathbf{A}=\gamma\mathbf{W}_{t}\otimes\mathbf{M}_{s}+\nu\mathbf{M}_{t}\otimes\mathbf{L}_{s},$
(3.6)

where for $i,j=1,\dots,n_{t}$

$[\mathbf{W}_{t}]_{i,j}=\int_{0}^{T}b^{\prime}_{j,{p}_{t}}(t)\,b_{i,{p}_{t}}(t)\,\mathrm{dt}\quad\text{and}\quad[\mathbf{M}_{t}]_{i,j}=\int_{0}^{T}\,b_{i,p_{t}}(t)\,b_{j,p_{t}}(t)\,\mathrm{dt},$
(3.7a) while for $i,j=1,\dots,N_{s}$
$[\mathbf{L}_{s}]_{i,j}=\int_{\Omega}\nabla
B_{i,p_{s}}(\boldsymbol{x})\cdot\nabla B_{j,p_{s}}(\boldsymbol{x})\
\mathrm{d}\Omega\quad\text{and}\quad[\mathbf{M}_{s}]_{i,j}=\int_{\Omega}B_{i,p_{s}}(\boldsymbol{x})\
B_{j,p_{s}}(\boldsymbol{x})\ \mathrm{d}\Omega.$ (3.7b)

In what follows we define and investigate the application of different
preconditioners for solving (3.5).

##  4 Galerkin preconditioner

As first attempt we introduce, for the system (3.5), the preconditioner

$[\widehat{\mathbf{A}}]_{i,j}:=\widehat{\mathcal{A}}(\widehat{B}_{j,\boldsymbol{p}},\widehat{B}_{i,\boldsymbol{p}}),$

where

$\widehat{\mathcal{A}}(\widehat{v},\widehat{w}):=\int_{0}^{1}\int_{\widehat{\Omega}}\left(\gamma\partial_{t}\widehat{v}\,\widehat{w}+T\nu\nabla\widehat{v}\cdot\,\nabla\widehat{w}\right)\,\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau\quad\forall\widehat{v},\widehat{w}\in\widehat{\mathcal{X}}_{h}.$

We have

$\widehat{\mathbf{A}}=\gamma{\mathbf{W}}_{t}\otimes\widehat{\mathbf{M}}_{s}+\nu{\mathbf{M}}_{t}\otimes{\widehat{\mathbf{L}}_{s}},$
(4.1)

where ${\widehat{\mathbf{L}}_{s}}$ and $\widehat{\mathbf{M}}_{s}$ are the
equivalent of (3.7b) in the parametric domain, i.e. we define for
$i,j=1,\dots,N_{s}$

$[{\widehat{\mathbf{L}}_{s}}]_{i,j}=\int_{\widehat{\Omega}}\nabla{\widehat{B}}_{i,p_{s}}(\boldsymbol{\eta})\cdot\nabla{\widehat{B}}_{j,p_{s}}(\boldsymbol{\eta})\
\mathrm{d}\widehat{\Omega}\quad\text{and}\quad[\widehat{\mathbf{M}}_{s}]_{i,j}=\int_{\widehat{\Omega}}\widehat{B}_{i,p_{s}}(\boldsymbol{\eta})\
\widehat{B}_{j,p_{s}}(\boldsymbol{\eta})\ \mathrm{d}\widehat{\Omega}.$ (4.2)

We emphasize that the time matrices appearing in (4.1) are the same ones
appearing in the system matrix (3.6). This is because for $i,j=1,\dots,n_{t}$
we have

$[{\mathbf{W}}_{t}]_{i,j}=\int_{0}^{T}b^{\prime}_{j,{p}_{t}}(t)\
b_{i,{p}_{t}}(t)\
\mathrm{dt}=\int_{0}^{1}\widehat{b}^{\prime}_{j,{p}_{t}}(\tau)\
\widehat{b}_{i,{p}_{t}}(\tau)\ \mathrm{d}\tau$

and

$[{\mathbf{M}}_{t}]_{i,j}=\int_{0}^{T}b_{j,{p}_{t}}(t)\ b_{i,{p}_{t}}(t)\
\mathrm{dt}=T\int_{0}^{1}\widehat{b}_{j,{p}_{t}}(\tau)\
\widehat{b}_{i,{p}_{t}}(\tau)\ \mathrm{d}\tau.$

Thanks to (2.1), the spatial matrices (4.2) have the following structure

${\widehat{\mathbf{L}}_{s}}=\sum_{l=1}^{d}\widehat{\mathbf{M}}_{d}\otimes\dots\otimes\widehat{\mathbf{M}}_{l+1}\otimes{\widehat{\mathbf{L}}_{l}}\otimes\widehat{\mathbf{M}}_{l-1}\otimes\dots\otimes\widehat{\mathbf{M}}_{1}\quad\text{and}\quad\widehat{\mathbf{M}}_{s}=\widehat{\mathbf{M}}_{d}\otimes\dots\otimes\widehat{\mathbf{M}}_{1},$
(4.3)

where for $l=1,\dots,d$ and for $i,j=1,\dots,n_{s,l}$ we define

$[{\widehat{\mathbf{L}}_{l}}]_{i,j}:=\int_{0}^{1}\widehat{b}^{\prime}_{i,p_{s}}(\eta_{k})\widehat{b}^{\prime}_{j,p_{s}}(\eta_{k})\mathrm{d}\eta_{k}\quad\text{and}\quad[\widehat{\mathbf{M}}_{l}]_{i,j}:=\int_{0}^{1}\widehat{b}_{i,p_{s}}(\eta_{k})\widehat{b}_{j,p_{s}}(\eta_{k})\mathrm{d}\eta_{k}.$

The efficient application of the proposed preconditioner, that is, the
solution of a linear system with matrix $\widehat{\mathbf{A}}$, should exploit
the structure highlighted above. When the pencils
$(\mathbf{W}_{t},\mathbf{M}_{t})$,
$({\widehat{\mathbf{L}}_{1}},\widehat{\mathbf{M}}_{1}),\ldots,({\widehat{\mathbf{L}}_{d}},\widehat{\mathbf{M}}_{d})$
admit a stable generalized eigendecomposition, a possible approach is the Fast
Diagonalization (FD) method, see [34] and [25] for details. We will see in
Section 4.1 that the spatial pencils
$({\widehat{\mathbf{L}}_{1}},\widehat{\mathbf{M}}_{1}),\ldots,({\widehat{\mathbf{L}}_{d}},\widehat{\mathbf{M}}_{d})$
admit a stable diagonalization, but this is not the case of
$(\mathbf{W}_{t},\mathbf{M}_{t})$, that needs a special treatment as explained
in Section 4.2.

###  4.1 Stable factorization of
${({\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l})}$ for ${l=1,\dots,d}$

The spatial stiffness and mass matrices ${\widehat{\mathbf{L}}_{l}}$ and
$\widehat{\mathbf{M}}_{l}$ are symmetric and positive definite for
$l=1,\dots,d$. Thus, the pencils
$({\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l})$ for $l=1,\dots,d$
admit the generalized eigendecomposition

${\widehat{\mathbf{L}}_{l}}\mathbf{U}_{l}=\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}\mathbf{\Lambda}_{l},$

where the matrices $\mathbf{U}_{l}$ contain in each column the
$\widehat{\mathbf{M}}_{l}$-orthonormal generalized eigenvectors and
$\mathbf{\Lambda}_{l}$ are diagonal matrices whose entries contain the
generalized eigenvalues. Therefore we have for $l=1,\dots,d$ the
factorizations

$\mathbf{U}^{T}_{l}{\widehat{\mathbf{L}}_{l}}\mathbf{U}_{l}=\mathbf{\Lambda}_{l}\quad\text{
and
}\quad\mathbf{U}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}=\mathbb{I}_{n_{s,l}},$
(4.4)

where $\mathbb{I}_{n_{s,l}}$ denotes the identity matrix of dimension
$n_{s,l}\times n_{s,l}$. Figure 1 shows the shape of the generalized
eigenvectors in $\mathbf{U}_{l}$, with associated eigenvalue in $\Lambda_{l}$,
for a fixed univariate direction $l=1,\dots,d$ discretized with degree
$p_{s}=3$ B-Splines and uniform partition with $n_{el}=32$. The stability of
the decomposition (4.4) is expressed by the condition number of the
eigenvector matrix. In particular
$\mathbf{U}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}=\mathbb{I}_{n_{s,l}}$
implies that

$\kappa_{2}(\mathbf{U}_{l}):=\|\mathbf{U}_{l}\|_{2}\|\mathbf{U}_{l}^{-1}\|_{2}=\sqrt{\kappa_{2}(\widehat{\mathbf{M}}_{l})},$

where $\|\cdot\|_{2}$ is the norm induced by the Euclidean vector norm. The
condition number $\kappa_{2}(\widehat{\mathbf{M}}_{l})$ has been studied in
[35] and it does not depend on the mesh-size, but it depends on the polynomial
degree. Indeed, we report in Table 1 the behavior of
$\kappa_{2}(\mathbf{U}_{l})$ for different values of spline degree $p_{s}$ and
for different uniform discretizations with number of elements denoted by
$n_{el}$. We observe that $\kappa_{2}(\mathbf{U}_{l})$ exhibits a dependence
only on $p_{s}$, but stays moderately low for all low polynomial degrees that
are in the range of interest.

Figure 1: Generalized eigenvectors in space with associated eigenvalues for $p_{s}=3$ and $n_{el}=32$ elements. The real part is expressed in blue, while the imaginary part (null) is in red. $n_{el}$ | $p_{s}=2$ | $p_{s}=3$ | $p_{s}=4$ | $p_{s}=5$ | $p_{s}=6$ | $p_{s}=7$ | $p_{s}=8$
---|---|---|---|---|---|---|---
0032 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0064 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0128 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0256 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0512 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
1024 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
Table 1: $\kappa_{2}(\mathbf{U}_{l})$ for different polynomial degrees $p_{s}$
and number of elements $n_{el}$.

###  4.2 Stable factorization of ${(\mathbf{W}_{t},\mathbf{M}_{t})}$

####  4.2.1 Numerical instability of the eigendecomposition

While $\mathbf{M}_{t}$ is symmetric, $\mathbf{W}_{t}$ is neither symmetric nor
skew-symmetric. Indeed

$\displaystyle[\mathbf{W}_{t}]_{i,j}+[\mathbf{W}_{t}]_{j,i}=\\!\\!\int_{0}^{T}b^{\prime}_{j,{p}_{t}}(t)\
b_{i,{p}_{t}}(t)\ \mathrm{dt}+\\!\\!\int_{0}^{T}b^{\prime}_{i,{p}_{t}}(t)\
b_{j,{p}_{t}}(t)\ \mathrm{dt}=b_{i,{p}_{t}}(T)\ b_{j,{p}_{t}}(T)$ (4.5)

where $b_{i,{p}_{t}}(T)\ b_{j,{p}_{t}}(T)$ vanishes for all
$i=1,\dots,n_{t}-1$ or $j=1,\dots,n_{t}-1$. A numerical computation of the
generalized eigendecomposition of the pencil
$(\mathbf{W}_{t},\mathbf{M}_{t})$, that is

$\mathbf{W}_{t}\mathbf{U}=\mathbf{M}_{t}\mathbf{U}\mathbf{\Lambda}_{t},$ (4.6)

where $\mathbf{\Lambda}_{t}$ is the diagonal matrix of the generalized complex
eigenvalues and ${\mathbf{U}}$ is the complex matrix whose columns are the
generalized eigenvectors normalized w.r.t. the norm induced by
${\mathbf{M}_{t}}$, reveals that the eigenvectors are far from
$\mathbf{M}_{t}$-orthogonality, i.e. the matrix
${\mathbf{U}}^{*}\mathbf{M}_{t}{\mathbf{U}}$ is not diagonal. In Figure 2 we
set $T=1$ and plot these generalized eigenvectors with associated eigenvalue,
for $p_{t}=3$ and uniform partition with $n_{el}=32$. We report in Table 2 the
condition number $\kappa_{2}(\mathbf{U})$ for different values of spline
degree $p_{t}$ and for different uniform discretizations with $n_{el}$ number
of elements. In contrast to the spatial case (see Section 4.1),
$\kappa_{2}(\mathbf{U})$ is large and grows exponentially with respect to the
spline degree $p_{t}$ and the level of mesh refinement. This test clearly
indicates a numerical instability when computing the generalized
eigendecomposition of $(\mathbf{W}_{t},\mathbf{M}_{t})$. A similar behavior
has also been highlighted in [15].

Figure 2: Generalized eigenvectors in time with associated eigenvalues for degree $3$ B-Splines with uniform partition of $32$ elements. The real part is expressed in blue, while the imaginary part is in red. $n_{el}$ | $p_{t}=2$ | $p_{t}=3$ | $p_{t}=4$ | $p_{t}=5$ | $p_{t}=6$ | $p_{t}=7$ | $p_{t}=8$
---|---|---|---|---|---|---|---
0032 | $8.9\cdot 10^{2}$ | $3.0\cdot 10^{4}$ | $5.0\cdot 10^{4\phantom{0}}$ | $3.4\cdot 10^{5\phantom{0}}$ | $3.1\cdot 10^{6\phantom{0}}$ | $4.2\cdot 10^{7\phantom{0}}$ | $7.0\cdot 10^{8\phantom{0}}$
0064 | $4.4\cdot 10^{3}$ | $2.6\cdot 10^{5}$ | $5.0\cdot 10^{5\phantom{0}}$ | $5.4\cdot 10^{6\phantom{0}}$ | $8.9\cdot 10^{7\phantom{0}}$ | $3.1\cdot 10^{9\phantom{0}}$ | $2.0\cdot 10^{10}$
0128 | $2.3\cdot 10^{4}$ | $1.2\cdot 10^{6}$ | $5.8\cdot 10^{6\phantom{0}}$ | $1.0\cdot 10^{8\phantom{0}}$ | $3.0\cdot 10^{9\phantom{0}}$ | $6.4\cdot 10^{11}$ | $1.3\cdot 10^{12}$
0256 | $1.2\cdot 10^{5}$ | $9.4\cdot 10^{6}$ | $7.6\cdot 10^{7\phantom{0}}$ | $2.1\cdot 10^{9\phantom{0}}$ | $1.2\cdot 10^{11}$ | $1.2\cdot 10^{13}$ | $2.1\cdot 10^{13}$
0512 | $7.0\cdot 10^{5}$ | $8.3\cdot 10^{7}$ | $1.1\cdot 10^{9\phantom{0}}$ | $4.9\cdot 10^{10}$ | $4.5\cdot 10^{12}$ | $3.6\cdot 10^{13}$ | $4.9\cdot 10^{12}$
1024 | $4.1\cdot 10^{6}$ | $8.0\cdot 10^{8}$ | $1.9\cdot 10^{10}$ | $1.3\cdot 10^{12}$ | $9.6\cdot 10^{12}$ | $1.4\cdot 10^{12}$ | $5.6\cdot 10^{12}$
Table 2:  $\kappa_{2}(\mathbf{U})$ for different degree $p_{t}$ and number of
elements $n_{el}$.

####  4.2.2 Construction of the stable factorization

The analysis above motivates the search of a different but stable
factorization of the pencil $(\mathbf{W}_{t},\mathbf{M}_{t})$. We look now for
a factorization of the form

$\mathbf{W}_{t}\mathbf{U}_{t}=\mathbf{M}_{t}\mathbf{U}_{t}\mathbf{\Delta}_{t},$
(4.7)

where $\mathbf{\Delta}_{t}$ is a complex arrowhead matrix, i.e. with non-zero
entries allowed on the diagonal, on the last row and on the last column only.
We also require that $\mathbf{U}_{t}$ fulfils the orthogonality condition

$\mathbf{U}_{t}^{*}\mathbf{M}_{t}{\mathbf{U}_{t}}=\mathbb{I}_{n_{t}}.$ (4.8)

From (4.7)–(4.8) we then obtain the factorizations

$\mathbf{U}_{t}^{*}\mathbf{W}_{t}\mathbf{U}_{t}=\mathbf{\Delta}_{t}\quad\text{
and }\quad\mathbf{U}_{t}^{*}\mathbf{M}_{t}\mathbf{U}_{t}=\mathbb{I}_{n_{t}}.$
(4.9)

With this aim, we look for $\mathbf{U}_{t}$ as follows:

$\mathbf{U}_{t}:=\begin{bmatrix}\overset{\circ}{\mathbf{U}}_{t}&\mathbf{r}\\\\[2.0pt]
\boldsymbol{0}^{T}&\rho\end{bmatrix}$ (4.10)

where
$\overset{\circ}{\mathbf{U}}_{t}\in\mathbb{C}^{(n_{t}-1)\times(n_{t}-1)}$,
$\mathbf{r}\in\mathbb{C}^{n_{t}-1}$, $\rho\in\mathbb{C}$ and where
$\boldsymbol{0}\in\mathbb{R}^{n_{t}-1}$ denotes the null vector. In order to
guarantee the non-singularity of $\mathbf{U}_{t}$, we further impose $\rho\neq
0$. Accordingly, we split the time matrices $\mathbf{W}_{t}$ and
$\mathbf{M}_{t}$ as

$\mathbf{W}_{t}=\begin{bmatrix}\overset{\circ}{\mathbf{W}}_{t}&\mathbf{w}\\\\[2.0pt]
-\mathbf{w}^{T}&\omega\end{bmatrix}\quad\text{and}\quad\mathbf{M}_{t}=\begin{bmatrix}\overset{\circ}{\mathbf{M}}_{t}&\mathbf{m}\\\\[2.0pt]
\mathbf{m}^{T}&\mu\end{bmatrix},$ (4.11)

where we have defined

$\omega:=[\mathbf{W}_{t}]_{n_{t},n_{t}},\qquad\mu:=[\mathbf{M}_{t}]_{n_{t},n_{t}},$
$[\mathbf{w}]_{i}=[\mathbf{W}_{t}]_{i,n_{t}}\quad\text{and}\quad[\mathbf{m}]_{i}=[\mathbf{M}_{t}]_{i,n_{t}}\quad\text{for}\quad
i=1,\dots,n_{t}-1,$
$\phantom{0}[\overset{\circ}{\mathbf{W}}_{t}]_{i,j}=[\mathbf{W}_{t}]_{i,j}\quad\text{and}\quad[\overset{\circ}{\mathbf{M}}_{t}]_{i,j}=[\mathbf{M}_{t}]_{i,j}\quad\text{for}\quad
i,j=1,\dots,n_{t}-1.$

Recalling (4.5), we observe that $\overset{\circ}{\mathbf{W}}_{t}$ is skew-
symmetric and, since $\overset{\circ}{\mathbf{M}}_{t}$ is symmetric, we can
write the eigendecomposition of the pencils
$(\overset{\circ}{\mathbf{W}}_{t},\overset{\circ}{\mathbf{M}}_{t})$:

$\overset{\circ}{\mathbf{W}}_{t}\overset{\circ}{\mathbf{U}}_{t}=\overset{\circ}{\mathbf{M}}_{t}\overset{\circ}{\mathbf{U}}_{t}\overset{\circ}{\mathbf{\Lambda}}_{t}\quad\text{
with
}\quad\overset{\circ\phantom{0}}{\mathbf{U}_{t}^{*}}\overset{\circ}{\mathbf{M}}_{t}\overset{\circ}{\mathbf{U}}_{t}=\mathbb{I}_{n_{t}-1},$
(4.12)

where $\overset{\circ}{\mathbf{U}}_{t}$ contains the complex generalized
eigenvectors and $\overset{\circ}{\mathbf{\Lambda}}_{t}$ is the diagonal
matrix of the generalized eigenvalues, that are pairs of complex conjugate
pure imaginary numbers plus, eventually, the eigenvalue zero. From
(4.10)–(4.11), it follows

$\mathbf{U}_{t}^{*}\mathbf{M}_{t}\mathbf{U}_{t}=\begin{bmatrix}\mathbb{I}_{n_{t}-1}&\overset{\circ\phantom{0}}{\mathbf{U}^{*}_{t}}\overset{\circ}{\mathbf{M}}_{t}\mathbf{r}+\overset{\circ\phantom{0}}{\mathbf{U}^{*}_{t}}\mathbf{m}\rho\\\
\mathbf{r}^{*}\overset{\circ}{\mathbf{M}}_{t}\overset{\circ}{\mathbf{U}}_{t}+\rho^{*}\mathbf{m}^{T}\overset{\circ}{\mathbf{U}}_{t}&\left[\mathbf{r}^{*}\rho^{*}\right]\mathbf{M}_{t}\begin{bmatrix}\mathbf{r}\\\
\rho\end{bmatrix}\end{bmatrix},$

where for the top-left block we have used (4.12).

The orthogonality condition in (4.8) holds if and only if $\mathbf{r}$ and
$\rho$ fulfil the two conditions:

$\displaystyle\overset{\circ\phantom{0}}{\mathbf{U}^{*}_{t}}\overset{\circ}{\mathbf{M}}_{t}\mathbf{r}+\overset{\circ\phantom{0}}{\mathbf{U}^{*}_{t}}\mathbf{m}\rho=\boldsymbol{0},$
(4.13a)
$\displaystyle\left[\mathbf{r}^{*}\rho^{*}\right]\mathbf{M}_{t}\begin{bmatrix}\mathbf{r}\\\
\rho\end{bmatrix}=1.$ (4.13b)

In order to compute $\mathbf{r}$ and $\rho$, we first find
$\mathbf{v}\in\mathbb{C}^{n_{t}-1}$ such that

$\overset{\circ}{\mathbf{M}}_{t}\mathbf{v}=-\mathbf{m};$ (4.14)

then we normalize the vector $\begin{bmatrix}\mathbf{v}\\\ 1\end{bmatrix}$
w.r.t. the $\|\cdot\|_{\mathbf{M}_{t}}$-norm to get

$\begin{bmatrix}\mathbf{r}\\\
\rho\end{bmatrix}:=\frac{\begin{bmatrix}\mathbf{v}\\\
1\end{bmatrix}}{\left([\mathbf{v}^{*}\
1]\mathbf{M}_{t}\begin{bmatrix}\mathbf{v}\\\
1\end{bmatrix}\right)^{\tfrac{1}{2}}}$

that fulfils (4.13a)–(4.13b). Finally, we get (4.7) by defining

$\mathbf{\Delta}_{t}:=\mathbf{U}_{t}^{*}\mathbf{W}_{t}\mathbf{U}_{t}=\begin{bmatrix}\overset{\circ}{\mathbf{\Lambda}}_{t}&\mathbf{g}\\\\[2.0pt]
-\mathbf{g}^{*}&\sigma\end{bmatrix},$ (4.15)

where
$\mathbf{g}:=\overset{\circ}{{\mathbf{U}}_{t}^{*}}\left[\overset{\circ}{\mathbf{W}}_{t}\
\mathbf{w}\right]\begin{bmatrix}\mathbf{r}\\\ \rho\end{bmatrix}$ and
$\sigma:=\left[{\mathbf{r}}^{*}{\rho}^{*}\right]\mathbf{W}_{t}\begin{bmatrix}\mathbf{r}\\\
\rho\end{bmatrix}$. Note that matrix (4.15) has an arrowhead structure. Figure
3 shows the plot of the columns of $\mathbf{U}_{t}$ with associated diagonal
entry of $\mathbf{\Delta}_{t}$, for $p_{t}=3$ and uniform partition with
$n_{el}=32$.

To assess the stability of the new decomposition (4.9), we set $T=1$ and we
compute the condition number $\kappa_{2}(\mathbf{U}_{t})$ for different values
of spline degree $p_{t}$ and for various uniform discretizations with number
of elements $n_{el}$. Thanks to (4.8), we have
$\kappa_{2}(\mathbf{U}_{t})=\sqrt{\kappa_{2}({\mathbf{M}_{t})}}$. The results,
reported in Table 3, show that the condition numbers
$\kappa_{2}(\mathbf{U}_{t})$ are uniformly bounded w.r.t. the mesh refinement,
they grow with respect to the polynomial degree but they are moderately small
for all the degrees of interest. We conclude that the factorization (4.9) for
the time pencil $(\mathbf{W}_{t},\mathbf{M}_{t})$ is stable.

Figure 3: Real part (blue) and imaginary part (red) of columns of $\mathbf{U}_{t}$ with associated diagonal entry in $\mathbf{\Delta}_{t}$. Discretization with $p_{t}=3$ and $n_{el}=32$. $n_{el}$ | $p_{t}=2$ | $p_{t}=3$ | $p_{t}=4$ | $p_{t}=5$ | $p_{t}=6$ | $p_{t}=7$ | $p_{t}=8$
---|---|---|---|---|---|---|---
0032 | $3.2\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0064 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0128 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0256 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0512 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
1024 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
Table 3: $\kappa_{2}(\mathbf{U}_{t})$ for different degrees $p_{t}$ and number
of elements $n_{el}$.

###  4.3 Preconditioner application

The application of the preconditioner involves the solution of the linear
system

$\widehat{\mathbf{A}}\mathbf{s}=\mathbf{r},$ (4.16)

where $\widehat{\mathbf{A}}$ has the structure (4.1). We are able to
efficiently solve system (4.16) by extending the FD method. The starting
points, that are involved in the setup of the preconditioner, are the
following ones:

  * 1.

for the pencils $({\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l})$ for
$l=1,\dots,d$ we have the factorizations (4.4);

  * 2.

for the pencil $(\mathbf{W}_{t},\mathbf{M}_{t})$ we have the factorization
(4.9).

Then, by defining
$\mathbf{U}_{s}:=\mathbf{U}_{d}\otimes\dots\otimes\mathbf{U}_{1}$ and
$\mathbf{\Lambda}_{s}:=\sum_{l=1}^{d}\mathbb{I}_{n_{s,d}}\otimes\dots\otimes\mathbb{I}_{n_{s,l+1}}\otimes\mathbf{\Lambda}_{l}\otimes\mathbb{I}_{n_{s,l-1}}\otimes\dots\otimes\mathbb{I}_{n_{s,1}}$,
we have for the matrix $\widehat{\mathbf{A}}$ the factorization

$\widehat{\mathbf{A}}=\left(\mathbf{U}_{t}^{*}\otimes\mathbf{U}_{s}^{T}\right)^{-1}\left({\gamma}\mathbf{\Delta}_{t}\otimes\mathbb{I}_{N_{s}}+{\nu}{\mathbb{I}}_{n_{t}}\otimes{\mathbf{\Lambda}}_{s}\right)\left(\mathbf{U}_{t}\otimes\mathbf{U}_{s}\right)^{-1}.$
(4.17)

Note that the second factor in (4.17) has the block-arrowhead structure

$\displaystyle{\gamma}\mathbf{\Delta}_{t}\otimes\mathbb{I}_{N_{s}}+{\nu}{\mathbb{I}}_{n_{t}}\otimes{\mathbf{\Lambda}}_{s}$
$\displaystyle=\begin{bmatrix}\mathbf{H}_{1}&&&\mathbf{B}_{1}\\\\[2.0pt]
&\ddots&&\vdots\\\\[2.0pt]
\quad&&\mathbf{H}_{n_{t}-1}&\mathbf{B}_{n_{t}-1}\\\\[2.0pt]
-\mathbf{B}^{*}_{1}&\ldots&-\mathbf{B}^{*}_{n_{t}-1}&\mathbf{H}_{n_{t}}\end{bmatrix}$
(4.18)

where $\mathbf{H}_{i}$ and $\mathbf{B}_{i}$are diagonal matrices defined as

$\mathbf{H}_{i}:={\gamma}[\mathbf{\Lambda}_{t}]_{i,i}\mathbb{I}_{N_{s}}+{\nu}\mathbf{\Lambda}_{s}\quad\text{
and
}\quad\mathbf{B}_{i}:={\gamma}[\mathbf{g}]_{i}\mathbb{I}_{N_{s}}\quad\text{
for }\quad i=1,\dots,n_{t}-1,$
$\mathbf{H}_{n_{t}}:={\gamma}\sigma{\mathbb{I}}_{N_{s}}+{\nu}\mathbf{\Lambda}_{s}.$

The matrix (4.18) has the following easy-to-invert block LU decomposition

$\displaystyle\gamma\mathbf{\Delta}_{t}\otimes$
$\displaystyle\mathbb{I}_{N_{s}}+{\nu}{\mathbb{I}}_{n_{t}}\otimes{\mathbf{\Lambda}}_{s}$
(4.19) $\displaystyle=\begin{bmatrix}\mathbb{I}_{N_{s}}&&&\\\ &\ddots&&\\\
&&\mathbb{I}_{N_{s}}&\\\
-\mathbf{B}^{*}_{1}\mathbf{H}_{1}^{-1}&\dots&-\mathbf{B}^{*}_{n_{t}-1}\mathbf{H}_{n_{t}-1}^{-1}&\mathbb{I}_{N_{s}}\end{bmatrix}\begin{bmatrix}\mathbf{H}_{1}&&&\mathbf{B}_{1}\\\
&\ddots&&\vdots\\\ &&\mathbf{H}_{n_{t}-1}&\mathbf{B}_{n_{t}-1}\\\
&&&\mathbf{S}\end{bmatrix}$

where
$\mathbf{S}:=\mathbf{H}_{n_{t}}+\sum_{i=1}^{n_{t}-1}\mathbf{B}^{*}_{i}\mathbf{H}_{i}^{-1}\mathbf{B}_{i}$
is a diagonal matrix.

Summarizing, the solution of (4.16) can be computed by the following
algorithm.

Algorithm 1 Extended FD

1:Compute the factorizations (4.4) and (4.9).

2:Compute
$\widetilde{\mathbf{q}}=(\mathbf{U}_{t}^{*}\otimes\mathbf{U}_{s}^{T})\mathbf{r}$.

3:Compute
$\widetilde{\mathbf{s}}=\left({\gamma}\mathbf{\Delta}_{t}\otimes\mathbb{I}_{N_{s}}+{\nu}{\mathbb{I}}_{n_{t}}\otimes{\mathbf{\Lambda}}_{s}\right)^{-1}\widetilde{\mathbf{q}}.$

4:Compute $\mathbf{s}=(\mathbf{U}_{t}\otimes\mathbf{U}_{s})\
\widetilde{\mathbf{s}}.$

###  4.4 Preconditioner robustness: partial inclusion of the geometry

The preconditioner (4.1) does not incorporate any information on the spatial
parametrization $\boldsymbol{F}$. Thus, the quality of the preconditioning
strategy may depend on the geometry map: we see this trend in the numerical
tests presented in the upper tables of [22, Tables 4,6]. However, we can
generalize (4.1) by including in the univariate spatial matrices
${\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l}$ for $l=1,\dots,d$ a
suitable approximation of $\boldsymbol{F}$, without increasing the asymptotic
computational cost. A similar approach has been used first in [36] for the
Stokes problem and in [21] for a least-squares formulation of the heat
equation. We briefly give an overview of this strategy.

Referring to Section 2.2 for the notation of the basis functions, we rewrite
the entries of the system matrix (3.5) in the parametric domain as

$\displaystyle[\mathbf{A}]_{i,j}$
$\displaystyle=\mathcal{A}({B}_{j,\boldsymbol{p}},{B}_{i,\boldsymbol{p}})$
$\displaystyle=\gamma\int_{0}^{1}\\!\int_{\widehat{\Omega}}\tfrac{1}{T}\partial_{\tau}{\widehat{B}}_{{j},\boldsymbol{p}}{\widehat{B}}_{{i},\boldsymbol{p}}|\text{det}(J_{\boldsymbol{G}})|\,\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau$
$\displaystyle+\int_{0}^{1}\\!\int_{\widehat{\Omega}}\nu(\nabla{\widehat{B}}_{{j},\boldsymbol{p}})^{T}J_{\boldsymbol{G}}^{-1}J_{\boldsymbol{G}}^{-T}\nabla{\widehat{B}}_{{i},\boldsymbol{p}}|\text{det}(J_{\boldsymbol{G}})|\,\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau$
$\displaystyle={\int_{0}^{1}\\!\int_{\widehat{\Omega}}\left[(\nabla{\widehat{B}}_{{j},\boldsymbol{p}})^{T}\
\ {\partial_{\tau}}{\widehat{B}}_{{j},\boldsymbol{p}}\right]\begin{bmatrix}\nu
T\ \mathbb{I}_{d}&\\\
&\gamma\end{bmatrix}\mathfrak{C}\left[(\nabla{\widehat{B}}_{{i},\boldsymbol{p}})^{T}\
\
{\widehat{B}}_{{i},\boldsymbol{p}}\right]^{T}\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau,}$

where

$\mathfrak{C}:=\begin{bmatrix}J_{\boldsymbol{F}}^{-1}J_{\boldsymbol{F}}^{-T}|\det(J_{\boldsymbol{F}})|&\\\
&|\text{det}(J_{\boldsymbol{F}})|\end{bmatrix}$

and where we used that
${B}_{i,\boldsymbol{p}}=\widehat{B}_{i,\boldsymbol{p}}\circ\boldsymbol{G}^{-1}$,
${B}_{j,\boldsymbol{p}}=\widehat{B}_{j,\boldsymbol{p}}\circ\boldsymbol{G}^{-1}$
and $|\text{det}(J_{\boldsymbol{G}})|=T|\text{det}(J_{\boldsymbol{F}})|$. The
construction of the preconditioners is based on the following approximation of
the diagonal entries only of $\mathfrak{C}$:

$\displaystyle[\mathfrak{C}(\boldsymbol{\eta})]_{l,l}\approx[\widetilde{\mathfrak{C}}(\boldsymbol{\eta})]_{l,l}:={\varphi_{1}(\eta_{1})\dots\varphi_{l-1}(\eta_{l-1})\Phi_{l}(\eta_{l})\varphi_{l+1}(\eta_{l+1})\dots\varphi_{d}(\eta_{d})}\,\text{for
}l=1,\dots,d,$ (4.20a)
$\displaystyle[\mathfrak{C}(\boldsymbol{\eta})]_{d+1,d+1}\approx[\widetilde{\mathfrak{C}}(\boldsymbol{\eta})]_{d+1,d+1}:={\varphi_{1}(\eta_{1})\dots\varphi_{d}(\eta_{d}).}\phantom{0}$
(4.20b)

In order to compute such an approximation, we interpolate the functions
$[\widetilde{\mathfrak{C}}(\boldsymbol{\eta})]_{l,l}$ in (4.20) by piecewise
constants in each element and we build the univariate factors $\varphi_{l}$
and $\Phi_{l}$ by using the separation of variables algorithm detailed in [21,
Appendix C]. The computational cost of the approximation above is proportional
to the number of elements in $\Omega$, that, when using smooth B-splines, is
almost equal to $N_{s}$, independent of $p_{s}$ and $p_{t}$ and thus
negligible in the whole iterative strategy.

Then we define

$[\widetilde{\mathbf{A}}]_{i,j}:=\int_{0}^{1}\\!\int_{\widehat{\Omega}}\left[(\nabla{\widehat{B}}_{{j},\boldsymbol{p}})^{T}\
\ \partial_{\tau}{\widehat{B}}_{{j},\boldsymbol{p}}\right]\begin{bmatrix}\nu
T\ \mathbb{I}_{d}&\\\
&\gamma\end{bmatrix}\widetilde{\mathfrak{C}}\left[(\nabla{\widehat{B}}_{{i},\boldsymbol{p}})^{T}\
\
{\widehat{B}}_{{i},\boldsymbol{p}}\right]^{T}\,\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau.$

The previous matrix maintains the same Kronecker structure as (4.1). Indeed we
have that

$\widetilde{\mathbf{A}}=\gamma\mathbf{W}_{t}\otimes\widetilde{\mathbf{M}}_{s}+\nu\mathbf{M}_{t}\otimes\widetilde{\mathbf{K}}_{s},$
(4.21)

where

$\widetilde{\mathbf{K}}_{s}:=\sum_{l=1}^{d}\widetilde{\mathbf{M}}_{d}\otimes\dots\otimes\widetilde{\mathbf{M}}_{l+1}\otimes\widetilde{\mathbf{K}}_{l}\otimes\widetilde{\mathbf{M}}_{l-1}\otimes\dots\otimes\widetilde{\mathbf{M}}_{1},\qquad\widetilde{\mathbf{M}}_{s}:=\widetilde{\mathbf{M}}_{d}\otimes\dots\otimes\widetilde{\mathbf{M}}_{1},$

and where for $l=1,\dots,d$ and for $i,j=1,\dots,n_{s,l}$ we define

$[\widetilde{\mathbf{K}}_{l}]_{i,j}:=\int_{0}^{1}{\Phi_{l}}(\eta_{l})\widehat{b}^{\prime}_{i,p_{s}}(\eta_{l})\widehat{b}^{\prime}_{j,p_{s}}(\eta_{l})\mathrm{d}\eta_{l}\
\text{ and }\
[\widetilde{\mathbf{M}}_{l}]_{i,j}:=\int_{0}^{1}{\varphi_{l}}(\eta_{l})\widehat{b}_{i,p_{s}}(\eta_{l})\widehat{b}_{j,p_{s}}(\eta_{l})\mathrm{d}\eta_{l}.$

We remark that the application of (4.21) can still be performed by Algorithm
1. Finally, we apply a diagonal scaling on $\widetilde{\mathbf{A}}$ and we
define the Galerkin preconditioner as

$\widehat{\mathbf{A}}^{\boldsymbol{G}}:=\mathbf{D}^{\tfrac{1}{2}}\widetilde{\mathbf{A}}\mathbf{D}^{\tfrac{1}{2}}$
(4.22)

where
$\displaystyle{[\mathbf{D}]_{i,i}:=\frac{[\mathbf{A}]_{i,i}}{[\widetilde{\mathbf{A}}]_{i,i}}}$
for $i=1,\dots,N_{dof}$.

####  4.4.1 The case of non-constant separable coefficients

We briefly discuss a generalization of the preconditioning strategy to the
case of non-constant equation coefficients $\gamma$ and $\nu$. We assume that
$\gamma$ and $\nu$ are positive functions defined over $\Omega\times[0,T]$ and
that they are separable in space and in time, i.e. we can write

$\gamma(\boldsymbol{x},t)=\gamma_{s}(\boldsymbol{x})\gamma_{t}(t),\qquad\nu(\boldsymbol{x},t)=\nu_{s}(\boldsymbol{x})\nu_{t}(t),$

with $\gamma_{s},\gamma_{t},\nu_{s}$ and $\nu_{t}$ positive functions.

Now, the first equation of (3.1) can be written as

$\gamma_{s}\partial_{t}u-\nabla\cdot\left(\frac{\nu_{t}}{\gamma_{t}}\nu_{s}\nabla
u\right)=\frac{f}{\gamma_{t}}.$

We discretize this equation as described in Section 3 and we generalize the
definition of the linear system (3.6) with

$\mathbf{A}:=\mathbf{W}_{t}\otimes\underline{\mathbf{M}}_{s}+\underline{\mathbf{M}}_{t}\otimes\underline{\mathbf{K}}_{s},$

where $\mathbf{W}_{t}$ is defined as in (3.7a), while for $i,j=1,\dots,n_{t}$

$[\underline{\mathbf{M}}_{t}]_{i,j}:=\int_{0}^{T}\,\frac{\nu_{t}(t)}{\gamma_{t}(t)}b_{i,p_{t}}(t)\,b_{j,p_{t}}(t)\,\mathrm{dt}$

and for $i,j=1,\dots,N_{s}$

$\displaystyle[\underline{\mathbf{M}}_{s}]_{i,j}:=\int_{\Omega}\gamma_{s}(\boldsymbol{x})B_{i,\boldsymbol{p}_{s}}(\boldsymbol{x})\
B_{j,\boldsymbol{p}_{s}}(\boldsymbol{x})\ \mathrm{d}\Omega\ \text{ and }\
[\underline{\mathbf{K}}_{s}]_{i,j}:=\int_{\Omega}\nu_{s}(\boldsymbol{x})\nabla
B_{i,\boldsymbol{p}_{s}}(\boldsymbol{x})\cdot\nabla
B_{j,\boldsymbol{p}_{s}}(\boldsymbol{x})\ \mathrm{d}\Omega.$

Then, the preconditioner that we propose is defined as in (4.22)

$\widehat{\mathbf{A}}^{\boldsymbol{G}}:=\mathbf{D}^{\tfrac{1}{2}}\widetilde{\mathbf{A}}\mathbf{D}^{\tfrac{1}{2}},$

but here we generalize (4.21) with
$\widetilde{\mathbf{A}}:=\mathbf{W}_{t}\otimes\
\breve{{\mathbf{M}}}_{s}+\underline{\mathbf{M}}_{t}\otimes\breve{{\mathbf{K}}}_{s}$,
where the matrices $\breve{{\mathbf{K}}}_{s}$ and $\breve{{\mathbf{M}}}_{s}$
are obtained by using an approximation technique analogous to the one
described previously in this section, with $\gamma_{s}$ and $\nu_{s}$ included
in the coefficient matrix $\mathfrak{C}$. The preconditioner
$\widetilde{\mathbf{A}}$ can still be applied as described in Section 4.3.
Note that, for this purpose, it is crucial that $\mathbf{W}_{t}$ does not
incorporate any time-dependent coefficient, since this would invalidate (4.5).

##  5 Discrete Least-Square preconditioner

The second attempt of solving (3.3) is by considering its discrete least-
square formulation, that is

$\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{A}\mathbf{u}=\mathbf{g},$ (5.1)

where we recall
$[\mathbf{A}]_{i,j}=\mathcal{A}(B_{j,\boldsymbol{p}},B_{i,\boldsymbol{p}})$,
$\mathbf{g}:=\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{f}$ with
$[\mathbf{f}]_{i}=\mathcal{F}(B_{i,\boldsymbol{p}})$ and
$\mathbf{M}=\mathbf{M}_{t}\otimes\mathbf{M}_{s}$ is the mass matrix with
$\mathbf{M}_{t}$ and $\mathbf{M}_{s}$ as in (3.7b). The computation of
$\mathbf{g}$ requires to invert the mass $\mathbf{M}$, which can be
efficiently performed with different methods, i.e., ad hoc sparse
approximations of the inverse [37, 38], or preconditioning with the parametric
mass or its extensions proposed in [39, 40, 41] or using low rank
approximations as in [42, 43, 44]. Here we iteratively invert the mass with
Conjugate Gradients and the preconditioner designed in [41].

The tensor-product structure of the isogeometric space (2.4) allows to write
the system matrix $\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{A}$ as sum of
Kronecker products of matrices as

$\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{A}\
=\gamma^{2}\mathbf{W}^{T}_{t}\mathbf{M}^{-1}_{t}\mathbf{W}_{t}\otimes\mathbf{M}_{s}+\nu^{2}\mathbf{M}_{t}\otimes\mathbf{L}_{s}\mathbf{M}^{-1}_{s}\mathbf{L}_{s}+\gamma\nu\left(\mathbf{W}_{t}+\mathbf{W}_{t}^{T}\right)\otimes\mathbf{L}_{s},$
(5.2)

where $\mathbf{W}_{t}$ and $\mathbf{L}_{s}$ are defined in (3.7b).

For the problem (5.1), we therefore introduce the preconditioner

$\widehat{\mathbf{P}}:=\gamma^{2}\mathbf{W}^{T}_{t}\mathbf{M}^{-1}_{t}\mathbf{W}_{t}\otimes\widehat{\mathbf{M}}_{s}+\nu^{2}\mathbf{M}_{t}\otimes\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}^{-1}_{s}\widehat{\mathbf{L}}_{s}+\gamma\nu\left(\mathbf{W}_{t}+\mathbf{W}_{t}^{T}\right)\otimes\widehat{\mathbf{L}}_{s},$
(5.3)

where $\widehat{\mathbf{L}}_{s}$ and $\widehat{\mathbf{M}}_{s}$ are defined in
(4.2). Again, the efficient application of the proposed preconditioner, that
is, the solution of a linear system with matrix $\widehat{\mathbf{P}}$, should
exploit the structure highlighted above.

Recall that the space pencils
$(\widehat{\mathbf{L}}_{l},\widehat{\mathbf{M}}_{l})$ for $l=1,\dots,d$ admit
the stable factorizations described in Section 4.1, that is

$\mathbf{U}^{T}_{l}\widehat{\mathbf{L}}_{l}\mathbf{U}_{l}=\mathbf{\Lambda}_{l}\quad\text{
and
}\quad\mathbf{U}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}=\mathbb{I}_{n_{s,l}},$

where $\mathbb{I}_{n_{s,l}}$ denotes the identity matrix of dimension
$n_{s,l}\times n_{s,l}$.

###  5.1 Stable factorization in time

The time matrices $\mathbf{W}^{T}_{t}\mathbf{M}^{-1}_{t}\mathbf{W}_{t}$ and
$\mathbf{M}_{t}$ are symmetric positive definite, therefore they admit a
generalized eigendecomposition of the kind

$\mathbf{W}^{T}_{t}\mathbf{M}^{-1}_{t}\mathbf{W}_{t}\widetilde{\mathbf{U}}_{t}=\widetilde{\mathbf{\Lambda}}_{t}\mathbf{M}_{t}\widetilde{\mathbf{U}}_{t},$
(5.4)

where the matrix $\widetilde{\mathbf{U}}_{t}$ contains in each column the
${\mathbf{M}}_{t}$-orthonormal generalized eigenvectors and
$\widetilde{\mathbf{\Lambda}}_{t}$ is the diagonal matrix whose entries
contain the generalized eigenvalues. We have the following factorization

$\widetilde{\mathbf{U}}_{t}^{T}\mathbf{W}^{T}_{t}\mathbf{M}^{-1}_{t}\mathbf{W}_{t}\widetilde{\mathbf{U}}_{t}=\widetilde{\mathbf{\Lambda}}_{t}\quad\text{
and
}\quad\widetilde{\mathbf{U}}_{t}^{T}\mathbf{M}_{t}\widetilde{\mathbf{U}}_{t}=\mathbb{I}_{n_{t}}.$
(5.5)

Figure 4 shows the shape of the generalized eigenvectors in
$\widetilde{\mathbf{U}}_{t}$, with associated eigenvalue in
$\widetilde{\mathbf{\Lambda}}_{t}$, for a fixed univariate direction
$l=1,\dots,d$ discretized with degree $p_{s}=3$ B-Splines and uniform
partition. The stability of the decomposition (5.4) is again expressed by the
condition number of the eigenvector matrix, and since
$\widetilde{\mathbf{U}}^{T}_{t}\mathbf{M}_{t}\widetilde{\mathbf{U}}_{t}=\mathbb{I}_{n_{t}}$,
it holds
$\kappa_{2}(\widetilde{\mathbf{U}}_{t})=\sqrt{\kappa_{2}(\mathbf{M}_{t})}$.

We investigate the stability of the diagonalization (5.4) by setting $T=1$ and
computing the condition number $\kappa_{2}(\widetilde{\mathbf{U}}_{t})$ for
different values of spline degree $p_{t}$ and for various uniform
discretizations with number of elements $n_{el}$. The results, reported in
Table 4, show that the condition numbers
$\kappa_{2}(\widetilde{\mathbf{U}_{t}})$ are uniformly bounded w.r.t. the mesh
refinement, they grow with respect to the polynomial degree but they are
moderately small for all the degrees of interest. We conclude that the
generalized diagonalization (5.4) for the time pencil
$(\mathbf{W}^{T}_{t}\mathbf{M}^{-1}_{t}\mathbf{W}_{t},\mathbf{M}_{t})$ is
stable.

Figure 4: Real part (blue) and imaginary part (null) of columns of $\widetilde{\mathbf{U}}_{t}$ with associated diagonal entry in $\widetilde{\mathbf{\Lambda}}_{t}$. Discretization with $p_{t}=3$ and $n_{el}=32$. $n_{el}$ | $p_{t}=2$ | $p_{t}=3$ | $p_{t}=4$ | $p_{t}=5$ | $p_{t}=6$ | $p_{t}=7$ | $p_{t}=8$
---|---|---|---|---|---|---|---
0032 | $3.2\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0064 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0128 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0256 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0512 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
1024 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
Table 4: $\kappa_{2}(\widetilde{\mathbf{U}}_{t})$ for different degrees
$p_{t}$ and number of elements $n_{el}$.

###  5.2 Preconditioner application

The application of the preconditioner involves the solution of the linear
system

$\widehat{\mathbf{P}}\mathbf{s}=\mathbf{r},$ (5.6)

where $\widehat{\mathbf{P}}$ has the structure (5.3). We are able to
efficiently solve system (5.6) by the FD method and the Sherman-Morrison
formula. The starting points, that are involved in the setup of the
preconditioner, are the following ones:

  * 1.

for the pencils $({\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l})$ with
$l=1,\dots,d$ we have the factorizations (4.4);

  * 2.

for the pencil
$(\mathbf{W}_{t}^{T}\mathbf{M}_{t}^{-1}\mathbf{W}_{t},\mathbf{M}_{t})$ we have
the factorization (5.4).

Then, define $\mathbf{U}_{s}:=\mathbf{U}_{d}\otimes\dots\otimes\mathbf{U}_{1}$
and
$\mathbf{\Lambda}_{s}:=\sum_{l=1}^{d}\mathbb{I}_{n_{s,d}}\otimes\dots\otimes\mathbb{I}_{n_{s,l+1}}\otimes\mathbf{\Lambda}_{l}\otimes\mathbb{I}_{n_{s,l-1}}\otimes\dots\otimes\mathbb{I}_{n_{s,1}}$.
Notice that
$\widehat{\mathbf{M}}_{s}^{-1}={\mathbf{U}}_{s}{\mathbf{U}}_{s}^{T}$,
therefore the matrix
$\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}$
admits the stable factorization

$\mathbf{U}^{T}_{s}\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}\mathbf{U}_{s}=\mathbf{\Lambda}_{s}^{2}.$
(5.7)

The preconditioner $\widehat{\mathbf{P}}$ admits the following factorization

$\widehat{\mathbf{P}}=\left(\mathbb{I}_{n_{t}}\otimes\mathbf{U}_{s}^{T}\right)^{-1}\\\
\left(\left(\widetilde{\mathbf{U}}_{t}^{T}\otimes\mathbb{I}_{N_{s}}\right)^{-1}\left(\gamma^{2}\mathbf{\Lambda}_{t}\otimes\mathbb{I}_{N_{s}}+\nu^{2}\mathbb{I}_{t}\otimes\mathbf{\Lambda}_{s}^{2}\right)\left(\widetilde{\mathbf{U}}_{t}\otimes\mathbb{I}_{N_{s}}\right)^{-1}+\gamma\nu\left(\mathbf{W}_{t}+\mathbf{W}^{T}_{t}\right)\otimes\mathbf{\Lambda}_{s}\right)\\\
\left(\mathbb{I}_{n_{t}}\otimes\mathbf{U}_{s}\right)^{-1}.$ (5.8)

Note that the second factor in (5.8) is sum of two matrices, the second one
being

$\gamma\nu\left(\mathbf{W}_{t}+\mathbf{W}^{T}_{t}\right)\otimes\mathbf{\Lambda}_{s}$
(5.9)

where $\left(\mathbf{W}_{t}+\mathbf{W}^{T}_{t}\right)$ is a rank 1 matrix that
can be written as
$\left(\mathbf{W}_{t}+\mathbf{W}^{T}_{t}\right)=\boldsymbol{e}_{n_{t}}\boldsymbol{e}_{n_{t}}^{T}$,
with $\boldsymbol{e}_{n_{t}}$ being the last element of the canonical basis of
$\mathbb{R}^{n_{t}}$. Here it is possible to introduce a vector
$\boldsymbol{v}_{n_{t}}\in\mathbb{R}^{n_{t}}$ such that

$\boldsymbol{e}_{n_{t}}\boldsymbol{e}_{n_{t}}^{T}=\widetilde{\mathbf{U}}_{t}^{-T}\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}\widetilde{\mathbf{U}}_{t}^{-1},$

by defining
$\boldsymbol{v}_{n_{t}}:=\widetilde{\mathbf{U}}_{t}^{T}\boldsymbol{e}_{n_{t}}$
that is the last column of $\widetilde{\mathbf{U}}_{t}^{T}$. Thus equation
(5.8) can be rewritten as

$\widehat{\mathbf{P}}=\left(\widetilde{\mathbf{U}}_{t}^{T}\otimes\mathbf{U}_{s}^{T}\right)^{-1}(\mathbf{H}+\mathbf{K})\left(\widetilde{\mathbf{U}}_{t}\otimes\mathbf{U}_{s}\right)^{-1}$
(5.10)

where we have introduced

$\mathbf{H}:=\left(\gamma^{2}\mathbf{\Lambda}_{t}\otimes\mathbb{I}_{N_{s}}+\nu^{2}\mathbb{I}_{t}\otimes\mathbf{\Lambda}_{s}^{2}\right),$

for the full diagonal term, and

$\mathbf{K}:=\gamma\nu\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}\otimes\mathbf{\Lambda}_{s}.$

Let us introduce $\mathbf{S}_{n_{t},N_{s}}$, the space and time shuffle
matrix, such that
$\mathbf{A}\otimes\mathbf{B}=\mathbf{S}_{n_{t},N_{s}}(\mathbf{B}\otimes\mathbf{A})\mathbf{S}_{n_{t},N_{s}}^{T}$
for all $\mathbf{A}\in\mathbb{R}^{n_{t}\times n_{t}}$ and
$\mathbf{B}\in\mathbb{R}^{N_{s}\times N_{s}}$. Under reshuffling with
$\mathbf{S}_{n_{t},N_{s}}$, we have the following block diagonal structure

$\begin{split}\mathbf{H}=\mathbf{S}_{n_{t},N_{s}}\begin{bmatrix}\mathbf{H}_{1}&&\\\
&\ddots&\\\
&&\mathbf{H}_{N_{s}}\end{bmatrix}\mathbf{S}_{n_{t},N_{s}}^{T},\quad\text{and}\quad\mathbf{K}=\mathbf{S}_{n_{t},N_{s}}\begin{bmatrix}\mathbf{K}_{1}&&\\\
&\ddots&\\\
&&\mathbf{K}_{N_{s}}\end{bmatrix}\mathbf{S}_{n_{t},N_{s}}^{T},\end{split}$

where, for $i=1,\dots,N_{s}$, the matrices $\mathbf{H}_{i}$ are diagonal
defined as
$\mathbf{H}_{i}:={\gamma^{2}}\mathbf{\Lambda}_{t}+{\nu^{2}}[\mathbf{\Lambda}_{s}^{2}]_{i,i}\otimes\mathbb{I}_{n_{t}}$
, while
$\mathbf{K}_{i}=\gamma\nu[\mathbf{\Lambda}_{s}]_{i,i}\otimes\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}$.

In order to invert $\mathbf{H}+\mathbf{K}$, it is now sufficient to reshuffle
the data, and invert the following independent $N_{s}$ problems of size
$n_{t}\times n_{t}$:

$\left(\mathbf{H}_{i}+\mathbf{K}_{i}\right)\mathbf{x}_{i}=\mathbf{y}_{i}\quad\text{
for }i=1,\dots,N_{s}.$

Notice that, each $\mathbf{K}_{i}$ is a rank 1 perturbation in the above
systems, therefore the Sherman-Morrison formula gives

$\displaystyle\mathbf{x}_{i}=\mathbf{H}_{i}^{-1}\mathbf{y}_{i}-\gamma\nu[\Lambda_{s}]_{i,i}\dfrac{\mathbf{H}_{i}^{-1}\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}\mathbf{H}_{i}^{-1}\mathbf{y}_{i}}{1+\gamma\nu[\Lambda_{s}]_{i,i}\boldsymbol{v}_{n_{t}}^{T}\mathbf{H}_{i}^{-1}\boldsymbol{v}_{n_{t}}}\quad\text{
for }i=1,\dots,N_{s}.$ (5.11)

Summarizing, the solution of (5.6) can be computed by the following algorithm.

Algorithm 2 Sherman-Morrison FD

1:Compute the factorizations (4.4) and (5.4).

2:Compute
$\widetilde{\mathbf{r}}=(\widetilde{\mathbf{U}}_{t}^{T}\otimes\mathbf{U}_{s}^{T})\mathbf{r}$.

3:Reshuffle ${\mathbf{q}}=\mathbf{S}_{n_{t},N_{s}}^{T}\widetilde{\mathbf{r}}.$

4:Compute
$[\widetilde{\mathbf{q}}]_{i}=\left(\mathbf{H}_{i}+\mathbf{K}_{i}\right)^{-1}\mathbf{q}_{i}\quad\text{
for }i=1,\dots,N_{s}.$

5:Reshuffle
$\widetilde{\mathbf{s}}=\mathbf{S}_{n_{t},N_{s}}\widetilde{\mathbf{q}}$

6:Compute $\mathbf{s}=(\widetilde{\mathbf{U}}_{t}^{\otimes}\mathbf{U}_{s})\
\widetilde{\mathbf{s}}.$

##  6 Least square formulation

In this section we present the Least Squares space-time formulation introduced
in [21], which is an alternative well posed space-time variational
formulation, w.r.t. the one presented in Section 3. Let us recall the model
problem (3.1) here. We seek for a solution $u$ such that

$\left\\{\begin{array}[]{rcllrcl}\gamma\partial_{t}u-\nabla\cdot(\nu\nabla
u)&=&f&\mbox{in }&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!(0,T),\\\\[1.0pt]
u&=&0&\mbox{on
}&\partial\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\![0,T],\\\\[1.0pt]
u&=&0&\mbox{in
}&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!\\{0\\},\end{array}\right.$ (6.1)

this time assuming $f\in L^{2}(\Omega\times(0,T))$, while $\gamma>0$ is the
heat capacity constant and $\nu>0$ is the thermal conductivity constant. We
define the space

$\mathcal{V}:=\\{v\in[(H^{1}_{0}(\Omega)\cap H^{2}(\Omega))\otimes
L^{2}(0,T)]\cap[L^{2}(\Omega)\otimes H^{1}(0,T)]\text{ s.t. }v=0\text{ on
}\Omega\times\\{0\\}\\},$

endowed with the norm

$\|v\|_{\mathcal{V}}^{2}:=\int_{0}^{T}\|\Delta
v(\cdot,t)\|^{2}_{L^{2}(\Omega)}\mathrm{dt}+\int_{0}^{T}\|\partial_{t}v(\cdot,t)\|^{2}_{L^{2}(\Omega)}\mathrm{dt}.$
(6.2)

The minimum regularity of the spline spaces that we assume is the following.

######  Assumption 4.

We assume that $p_{s}\geq 2$, $p_{t}\geq 1$ and that
$\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{h_{s}}\subset
C^{1}(\widehat{\Omega})$ and $\widehat{\mathcal{S}}^{{p}_{t}}_{h_{t}}\subset
C^{0}\left((0,1)\right)$.

Under Assumptions 3 and 4, $\mathcal{V}$ is a Hilbert space and the
$\|\cdot\|_{\mathcal{V}}$-norm is equivalent to

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}:=\|v\|^{2}_{H^{2}(\Omega)\otimes
L^{2}(0,T)}+\|v\|^{2}_{L^{2}(\Omega)\otimes H^{1}(0,1)}.$ (6.3)

The Least Squares space-time variational formulation for system (3.1) reads:
find $u\in\mathcal{V}$ such that

$u=\operatorname*{arg\,min}_{v\in\mathcal{V}}\dfrac{1}{2}\|\gamma\partial_{t}v-\nu\Delta
v-f\|^{2}_{L^{2}(\Omega\times(0,T))}.$ (6.4)

Its Euler-Lagrange equation is

$\mathcal{B}(u,v)=\mathcal{F}(v),\quad\forall v\in\mathcal{V},$ (6.5)

where the bilinear form $\mathcal{B}(\cdot,\cdot)$ and the linear form
$\mathcal{F}(\cdot)$ are defined as

$\begin{split}\mathcal{B}(v,w)&:=\int_{0}^{T}\int_{\Omega}(\gamma^{2}\partial_{t}v\partial_{t}w+\nu^{2}\Delta
v\Delta w-\gamma\nu\partial_{t}v\Delta w-\gamma\nu\Delta
v\partial_{t}w)\mathrm{d}\Omega\,\mathrm{dt},\\\
\mathcal{F}(w)&:=\int_{0}^{T}\int_{\Omega}f(\gamma\partial_{t}w-\nu\Delta
w)\mathrm{d}\Omega\,\mathrm{dt}.\end{split}$ (6.6)

Notice that $\mathcal{B}$ is a $\mathcal{V}$-elliptic continuous bilinear form
and $\mathcal{F}$ is continuous in $\mathcal{V}$. Therefore, the well-
posedness of the variational formulation above is a classical result that
follows from Lax-Milgram Theorem.

The previous setting can be generalized to non-homogeneous initial and
boundary conditions. For example, suppose that in (3.1) we have the initial
condition $u=u_{0}\in\Omega\times\\{0\\}$, with $u_{0}\in H^{1}_{0}(\Omega)$,
we lift $u_{0}$ to $\widetilde{u}_{0}\in(H^{1}_{0}(\Omega)\cap
H^{2}(\Omega))\otimes L^{2}(0,T)\cap L^{2}(\Omega)\otimes H^{1}(0,T)$. Then
$\widetilde{u}=u-\widetilde{u}_{0}\in\mathcal{V}$ is the solution of

$\left\\{\begin{array}[]{rcllrcl}\gamma\partial_{t}\widetilde{u}-\nabla\cdot(\nu\nabla\widetilde{u})&=&\widetilde{f}&\mbox{in
}&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!(0,T),\\\\[1.0pt]
\widetilde{u}&=&0&\mbox{on
}&\partial\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\![0,T],\\\\[1.0pt]
\widetilde{u}&=&0&\mbox{in
}&\Omega\\!\\!\\!\\!&\times&\\!\\!\\!\\!\\{0\\},\end{array}\right.$ (6.7)

where
$\widetilde{f}:=f-\gamma\partial_{t}\widetilde{u}_{0}+\nabla\cdot(\nu\nabla\widetilde{u}_{0})$.
For a detailed description of the variational formulation of problems
(6.1)-(6.7) and their well-posedness see, for example, [32, 45].

###  6.1 Space-time Galerkin method

Denote by $\mathcal{V}_{h}:=\mathcal{X}_{h}$ endowed with the
$\|\cdot\|_{\mathcal{V}}$-norm. Thanks to Assumption 4, it holds

$\mathcal{V}_{h}\subset(H^{1}_{0}(\Omega)\cap H^{2}(\Omega))\otimes
H^{1}(0,T)\subset\mathcal{V}.$ (6.8)

Therefore, we consider a Galerkin method for (6.5), that is, the least square
approximation of the system (6.1): find $u_{h}\in\mathcal{X}_{h}$ such that

$u_{h}=\operatorname*{arg\,min}_{v_{h}\in\mathcal{V}_{h}}\dfrac{1}{2}\|\gamma\partial_{t}v_{h}-\nu\Delta
v_{h}-f\|^{2}_{L^{2}(\Omega\times(0,T))}.$ (6.9)

Its Euler-Lagrange equation is

$\mathcal{B}(u_{h},v_{v})=\mathcal{F}(v_{h}),\quad\forall
v_{h}\in\mathcal{V}_{h},$ (6.10)

Well-posedness and quasi optimality follow from standard arguments.

######  Theorem 3.

The minimization problem (6.9) and the variational problem (6.10) are
equivalent and they admit a unique solution $u_{h}\in\mathcal{V}_{h}$. It also
holds:

$\|u-u_{h}\|_{\mathcal{V}}\leq\sqrt{2}\inf_{v_{h}\in\mathcal{V}_{h}}\|u-v_{h}\|_{\mathcal{V}}$
(6.11)

###### Proof.

The proof of the equivalence and the existence and uniqueness of the solution
follow by using Lax-Milgram Theorem, while the proof of (6.11) is a
consequence of the Ceà Lemma and the symmetry of the bilinear form
$\mathcal{B}$. ∎

We have then the following a-priori estimate for h-refinement.

######  Theorem 4.

Let $q_{s}$ and $q_{t}$ be two positive integers such that $q_{s}\geq 2$ and
$q_{t}\geq 1$. If $u\in\mathcal{V}\cap(H^{q_{s}}(\Omega)\otimes
H^{1}(0,T))\cap(H^{2}(\Omega)\otimes H^{q_{t}}(0,T))$ is the solution of (6.1)
and $u_{h}\in\mathcal{V}_{h}$ is the solution of (6.10), then

$\|u-u_{h}\|_{\mathcal{V}}\leq
C(h_{s}^{k_{s}-2}\|u\|_{H^{k_{s}}(\Omega)\otimes
H^{1}(0,T)}+h_{t}^{k_{t}-1}\|u\|_{H^{2}(\Omega)\otimes H^{k_{t}}(0,T)}),$
(6.12)

where $k_{s}:=\min\\{q_{s},p_{s}+1\\}$, $k_{t}:=\min\\{q_{t},p_{t}+1\\}$,$C$
is a constant that depends only on $p_{s}$,$p_{t}$,$\alpha$ and the
parameterization $\boldsymbol{G}$

The result follows from the anisotropic approximation estimates that are
developed in [33]. An overview of the proof is given in [21].

###  6.2 Discrete system

Before introducing the discrete system, we rewrite the bilinear form
$\mathcal{B}(\cdot,\cdot)$ in an equivalent way, through the following Lemma.

######  Lemma 1.

The bilinear form $\mathcal{B}(\cdot,\cdot)$ can be written as

$\displaystyle\mathcal{B}(v_{h},w_{h})=\gamma^{2}\int_{0}^{T}\int_{\Omega}\partial_{t}v_{h}\partial_{t}w_{h}\mathrm{d}\Omega\,\mathrm{dt}++\nu^{2}\int_{0}^{T}\int_{\Omega}\Delta
v_{h}\Delta w_{h}\mathrm{d}\Omega\,\mathrm{dt}++\gamma\nu\int_{\Omega}\nabla
v_{h}(\boldsymbol{x},T)\cdot\nabla w_{h}(\boldsymbol{x},T)\mathrm{d}\Omega,$
(6.13)

for all $v_{h},w_{h}\in\mathcal{V}_{h}$.

###### Proof.

Let $v_{h},w_{h}\in\mathcal{V}_{h}$. First note that
$\partial_{t}v_{h},\partial_{t}w_{h}\in(H^{1}_{0}(\Omega)\cap
H^{2}(\Omega))\otimes L^{2}(0,T)$, from (6.8), and
$\partial_{t}v_{h},\partial_{t}w_{h}=0$ on $\partial\Omega\times[0,T]$. Using
Green’s formula and integrating by parts yields to

$\displaystyle\begin{split}-\int_{0}^{T}\int_{\Omega}(\partial_{t}v_{h}\Delta
w_{h}+\Delta v_{h}\partial_{t}w_{h})\mathrm{d}\Omega\,\mathrm{dt}&=\\\
-\int_{0}^{T}\int_{\partial\Omega}(\partial_{t}v_{h}\nabla w_{h}+\nabla
v_{h}\partial_{t}w_{h})\cdot\boldsymbol{n}\mathrm{d}\Omega\,\mathrm{dt}+\int_{0}^{T}\int_{\Omega}[\nabla(\partial_{t}v_{h})\nabla
w_{h}+\nabla v_{h}\nabla(\partial_{t}w_{h})]\mathrm{d}\Omega\,\mathrm{dt}&=\\\
+\int_{0}^{T}\partial_{t}\left[\int_{\Omega}\nabla v_{h}\cdot\nabla
w_{h}\mathrm{d}\Omega\right]\,\mathrm{dt}&=\\\ +\int_{\Omega}[\nabla
v_{h}(\boldsymbol{x},T)\cdot\nabla w_{h}(\boldsymbol{x},T)-\nabla
v_{h}(\boldsymbol{x},0)\cdot\nabla
w_{h}(\boldsymbol{x},0)]\mathrm{d}\Omega&=\\\ +\int_{\Omega}[\nabla
v_{h}(\boldsymbol{x},T)\cdot\nabla
w_{h}(\boldsymbol{x},T)]\mathrm{d}\Omega&,\end{split}$

where $\boldsymbol{n}\in\mathbb{R}^{d}$ is the external normal unit vector to
$\partial\Omega$. Then (6.13) follows. ∎

######  Remark 2.

Note that the identity (6.13) holds also in the continuous setting (see [21,
Appendix B]).

The linear system associate to (6.10) is

$\mathbf{B}\mathbf{u}=\mathbf{f},$ (6.14)

where
$[\mathbf{B}]_{i,j}:=\mathcal{B}(B_{i,\boldsymbol{p}},B_{j,\boldsymbol{p}})$
and $[\mathbf{f}]_{i}:=\mathcal{F}(B_{i,\boldsymbol{p}})$. The discrete system
matrix $\mathbf{B}$ can be written as sum of Kronecker product matrices

$\mathbf{B}=\gamma^{2}\mathbf{L}_{t}\otimes\mathbf{M}_{s}+\nu^{2}\mathbf{M}_{t}\otimes\mathbf{J}_{s}+\gamma\nu\mathbf{R}_{t}\otimes\mathbf{L}_{s},$
(6.15) where the time matrices involved are, for $i,j=1,\dots,n_{t}$
$\displaystyle[\mathbf{L}_{t}]_{i,j}:=\int_{0}^{T}b^{\prime}_{j,{p}_{t}}(t)\,b^{\prime}_{i,{p}_{t}}(t)\,\mathrm{dt},\quad[\mathbf{M}_{t}]_{i,j}:=\int_{0}^{T}\,b_{i,p_{t}}(t)\,b_{j,p_{t}}(t)\,\mathrm{dt},\quad[\mathbf{R}_{t}]_{i,j}:=b_{i,p_{t}}(T)\,b_{j,p_{t}}(T),$
(6.16a) and the space matrices are $i,j=1,\dots,N_{s}$

$\begin{split}[{\mathbf{J}_{s}}]_{i,j}:=\int_{\Omega}\Delta
B_{i,p_{s}}(\boldsymbol{x})\Delta&B_{j,p_{s}}(\boldsymbol{x})\
\mathrm{d}\Omega,\quad[\mathbf{M}_{s}]_{i,j}:=\int_{\Omega}B_{i,p_{s}}(\boldsymbol{x})\
B_{j,p_{s}}(\boldsymbol{x})\ \mathrm{d}\Omega,\\\
[{\mathbf{L}_{s}}]_{i,j}&:=\int_{\Omega}\nabla
B_{i,p_{s}}(\boldsymbol{x})\cdot\nabla B_{j,p_{s}}(\boldsymbol{x})\
\mathrm{d}\Omega.\end{split}$ (6.16b)

##  7 Quadratic preconditioner

The matrix $\mathbf{B}$ in (6.14) is symmetric and positive definite. Thus, we
design and analyze a suitable symmetric positive definite preconditioner to be
used for a preconditioned Conjugate Gradients method. Recall
$\widehat{\mathcal{X}}_{h}$ is the spline space defined in Section 2.2,
satisfying the regularity condition of Assumption 4.

The definition of the preconditioner is associated with the bilinear form
$\widehat{\mathcal{B}}:\widehat{\mathcal{X}}_{h}\times\widehat{\mathcal{X}}_{h}\to\mathbb{R}$
defined as

$\widehat{\mathcal{B}}(\widehat{v}_{h},\widehat{w}_{h}):=\gamma^{2}\int_{0}^{1}\int_{\widehat{\Omega}}\partial_{\tau}\widehat{v}_{h}\partial_{\tau}\widehat{w}_{h}\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau+\nu^{2}\sum_{l=1}^{d}\int_{0}^{1}\int_{\widehat{\Omega}}\dfrac{\partial^{2}\widehat{v}_{h}}{\partial\eta_{l}^{2}}\dfrac{\partial^{2}\widehat{w}_{h}}{\partial\eta_{l}^{2}}\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau$
(7.1)

and with the corresponding norm

$\|\widehat{v}_{h}\|^{2}_{\widehat{\mathcal{B}}}:=\widehat{\mathcal{B}}(\widehat{v}_{h},\widehat{v}_{h}).$
(7.2)

Thus, the definition of the preconditioner is

$[\widehat{\mathbf{B}}]_{i,j}:=\widehat{\mathcal{B}}(\widehat{B}_{i,\boldsymbol{p}},\widehat{B}_{j,\boldsymbol{p}}),\quad
i,j=1,\dots,N_{dof},$ (7.3)

and has the following structure

$\widehat{\mathbf{B}}=\gamma^{2}\widehat{\mathbf{L}}_{t}\otimes\widehat{\mathbf{M}}_{s}+\nu^{2}\widehat{\mathbf{M}}_{t}\otimes\widetilde{\mathbf{J}}_{s},$
(7.4)

where for $i,j=1,\dots,n_{t}$

$[\widehat{\mathbf{L}}_{t}]_{i,j}:=\int_{0}^{1}\widehat{b}^{\prime}_{i,p_{t}}(\tau)\widehat{b}^{\prime}_{j,p_{t}}(\tau)\mathrm{d}\tau,\quad[\widehat{\mathbf{M}}_{t}]_{i,j}:=\int_{0}^{1}\widehat{b}_{i,p_{t}}(\tau)\widehat{b}_{j,p_{t}}(\tau)\mathrm{d}\tau,$

while for $i,j=1,\dots,N_{s}$

$[\widetilde{\mathbf{J}}_{s}]_{i,j}:=\sum_{l=1}^{d}\int_{0}^{1}\dfrac{\partial^{2}\widehat{B}^{\prime}_{i,\boldsymbol{p}_{s}}(\boldsymbol{\eta})}{\partial\eta_{l}^{2}}\dfrac{\partial^{2}\widehat{B}^{\prime}_{j,\boldsymbol{p}_{s}}(\boldsymbol{\eta})}{\partial\eta_{l}^{2}}\mathrm{d}\widehat{\Omega},\quad[\widehat{\mathbf{M}}_{s}]_{i,j}:=\int_{0}^{1}\widehat{B}_{i,\boldsymbol{p}_{s}}(\boldsymbol{\eta})\widehat{B}_{j,\boldsymbol{p}_{s}}(\boldsymbol{\eta})\mathrm{d}\widehat{\Omega}-$

Notice that $\widehat{\mathbf{L}}_{t},\widehat{\mathbf{M}}_{t}$ and
$\widehat{\mathbf{M}}_{s}$ correspond to ${\mathbf{L}}_{t},{\mathbf{M}}_{t}$
and $\mathbf{M}_{s}$, respectively, where the integration is performed on the
parametric domain $\widehat{\Omega}$. Moreover, the matrices
$\widetilde{\mathbf{J}}_{s}$ and $\widehat{\mathbf{M}}_{s}$ can be further
factorized as sum of Kronecker products as

$[\widetilde{\mathbf{J}}_{s}]_{i,j}=\sum_{l=1}^{d}\widehat{\mathbf{M}}_{d}\otimes\dots\otimes\widehat{\mathbf{M}}_{l-1}\otimes\widehat{\mathbf{J}}_{l}\otimes\widehat{\mathbf{M}}_{l-1}\otimes\dots\otimes\widehat{\mathbf{M}}_{1},\quad\widehat{\mathbf{M}}_{s}=\widehat{\mathbf{M}}_{d}\otimes\dots\otimes\widehat{\mathbf{M}}_{1},$

where for $l=1,\dots,d$ and for $i,j=1,\dots,n_{s,l}$

$[\widehat{\mathbf{J}}_{l}]_{i,j}:=\int_{0}^{1}\widehat{b}^{\prime\prime}_{i,p_{s}}(\eta_{l})\widehat{b}^{\prime\prime}_{j,p_{s}}(\eta_{l})\mathrm{d}\eta_{l},\quad[\widehat{\mathbf{M}}_{l}]_{i,j}:=\int_{0}^{1}\widehat{b}_{i,p_{s}}(\eta_{l})\widehat{b}_{j,p_{s}}(\eta-l)\mathrm{d}\eta_{l}.$

The efficient application of the proposed preconditioner, that is, the
solution of a liner system with matrix $\widehat{\mathbf{B}}$, should exploit
the structure highlighted above. Again a possible approach is Fast
Diagonalization method.

Finally, the following spectral stability of the preconditioned matrix
$\widehat{\mathbf{B}}^{-1}\mathbf{B}$ holds true.

######  Theorem 5.

Under Assumptions 13 and 4, it holds

$\theta\leq\lambda_{min}(\widehat{\mathbf{B}}^{-1}\mathbf{B}),\quad\lambda_{max}(\widehat{\mathbf{B}}^{-1}\mathbf{B})\leq\Theta,$

where $\theta$ and $\Theta$ are positive constants that do not depend on
$h_{s},h_{t},p_{s},p_{t}$ but dependent on $\boldsymbol{G}$.

For the proof we refer to [21].

###  7.1 Stable factorization of
${({\widehat{\mathbf{J}}_{l}},\widehat{\mathbf{M}}_{l})}$ for ${l=1,\dots,d}$

The spatial matrices ${\widehat{\mathbf{J}}_{l}}$ and
$\widehat{\mathbf{M}}_{l}$ are symmetric and positive definite for
$l=1,\dots,d$. Thus, the pencils
$({\widehat{\mathbf{J}}_{l}},\widehat{\mathbf{M}}_{l})$ for $l=1,\dots,d$
admit the generalized eigendecomposition

${\
\widehat{\mathbf{J}}_{l}}\mathbf{V}_{l}=\widehat{\mathbf{M}}_{l}\mathbf{V}_{l}\mathbf{\Lambda}_{l},$

where the matrices $\mathbf{V}_{l}$ contain in each column the
$\widehat{\mathbf{M}}_{l}$-orthonormal generalized eigenvectors and
$\mathbf{\Lambda}_{l}$ are diagonal matrices whose entries contain the
generalized eigenvalues. Therefore we have for $l=1,\dots,d$ the
factorizations

$\mathbf{V}^{T}_{l}{\widehat{\mathbf{J}}_{l}}\mathbf{V}_{l}=\mathbf{\Lambda}_{l}\quad\text{
and
}\quad\mathbf{V}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{V}_{l}=\mathbb{I}_{n_{s,l}},$
(7.5)

where $\mathbb{I}_{n_{s,l}}$ denotes the identity matrix of dimension
$n_{s,l}\times n_{s,l}$. Figure 5 shows the shape of the generalized
eigenvectors in $\mathbf{V}_{l}$, with associated eigenvalue in $\Lambda_{l}$,
for a fixed univariate direction $l=1,\dots,d$ discretized with degree
$p_{s}=3$ B-Splines and uniform partition with $n_{el}=32$. The stability of
the decomposition (7.5) is expressed by the condition number of the
eigenvector matrix. In particular
$\mathbf{V}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{V}_{l}=\mathbb{I}_{n_{s,l}}$
implies that
$\kappa_{2}(\mathbf{V}_{l})=\sqrt{\kappa_{2}(\widehat{\mathbf{M}}_{l})}$ and
it does not depend on the mesh-size, but it depends on the polynomial degree.
Indeed, we report in Table 5 the behavior of $\kappa_{2}(\mathbf{V}_{l})$ for
different values of spline degree $p_{s}$ and for different uniform
discretizations with number of elements denoted by $n_{el}$. We observe that
$\kappa_{2}(\mathbf{V}_{l})$ exhibits a dependence only on $p_{s}$, but stays
moderately low for all low polynomial degrees that are in the range of
interest.

Figure 5: Generalized eigenvectors in space with associated eigenvalues for $p_{s}=3$ and $n_{el}=32$ elements. The real part is expressed in blue, while the imaginary part (null) is in red. $n_{el}$ | $p_{s}=2$ | $p_{s}=3$ | $p_{s}=4$ | $p_{s}=5$ | $p_{s}=6$ | $p_{s}=7$ | $p_{s}=8$
---|---|---|---|---|---|---|---
0032 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0064 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0128 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0256 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0512 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
1024 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
Table 5: $\kappa_{2}(\mathbf{V}_{l})$ for different polynomial degrees $p_{s}$
and number of elements $n_{el}$.

###  7.2 Stable factorization in time

The time matrices $\widehat{\mathbf{L}}_{t}$ and $\widehat{\mathbf{M}}_{t}$
are symmetric positive definite, therefore they admit a generalized
eigendecomposition of the kind

$\widehat{\mathbf{L}}_{t}\mathbf{V}_{t}=\mathbf{\Lambda}_{t}\widehat{\mathbf{M}}_{t}\mathbf{V}_{t},$
(7.6)

where the matrix $\mathbf{V}_{t}$ contains in each column the
${\widehat{\mathbf{M}}}_{t}$-orthonormal generalized eigenvectors and
$\mathbf{\Lambda}_{t}$ is the diagonal matrix whose entries contain the
generalized eigenvalues. We have the following factorization

$\mathbf{V}_{t}^{T}\widehat{\mathbf{L}}_{t}\mathbf{V}_{t}=\mathbf{\Lambda}_{t}\quad\text{
and
}\quad\mathbf{V}_{t}^{T}\widehat{\mathbf{M}}_{t}\mathbf{V}_{t}=\mathbb{I}_{n_{t}}.$
(7.7)

Figure 6 shows the shape of the generalized eigenvectors in $\mathbf{V}_{t}$,
with associated eigenvalue in $\mathbf{\Lambda}_{t}$, for a fixed univariate
direction $l=1,\dots,d$ discretized with degree $p_{s}=3$ B-Splines and
uniform partition. The stability of the decomposition (7.6) is again expressed
by the condition number of the eigenvector matrix, and since
$\mathbf{V}^{T}_{t}\widehat{\mathbf{M}}_{t}\mathbf{V}_{t}=\mathbb{I}_{n_{t}}$,
it holds
$\kappa_{2}(\mathbf{V}_{t})=\sqrt{\kappa_{2}(\widehat{\mathbf{M}}_{t})}$.

We investigate the stability of the diagonalization (7.6) by computing the
condition number $\kappa_{2}(\mathbf{V}_{t})$ for different values of spline
degree $p_{t}$ and for various uniform discretizations with number of elements
$n_{el}$. The results, reported in Table 6, show that the condition numbers
$\kappa_{2}({\mathbf{V}_{t}})$ are uniformly bounded w.r.t. the mesh
refinement, they grow with respect to the polynomial degree but they are
moderately small for all the degrees of interest. We conclude that the
generalized diagonalization (7.6) for the time pencil
$(\widehat{\mathbf{L}}_{t},\widehat{\mathbf{M}}_{t})$ is stable.

Figure 6: Real part (blue) and imaginary part (null) of columns of $\mathbf{V}_{t}$ with associated diagonal entry in $\mathbf{\Lambda}_{t}$. Discretization with $p_{t}=3$ and $n_{el}=32$. $n_{el}$ | $p_{t}=2$ | $p_{t}=3$ | $p_{t}=4$ | $p_{t}=5$ | $p_{t}=6$ | $p_{t}=7$ | $p_{t}=8$
---|---|---|---|---|---|---|---
0032 | $3.2\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0064 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0128 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0256 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
0512 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
1024 | $3.3\cdot 10^{0}$ | $5.2\cdot 10^{0}$ | $8.3\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.2\cdot 10^{1}$ | $3.6\cdot 10^{1}$ | $5.9\cdot 10^{1}$
Table 6: $\kappa_{2}(\mathbf{V}_{t})$ for different degrees $p_{t}$ and number
of elements $n_{el}$.

###  7.3 Preconditioner application

The application of the preconditioner involves the solution of the linear
system

$\widehat{\mathbf{B}}^{-1}\mathbf{s}=\mathbf{r},$ (7.8)

where $\widehat{\mathbf{B}}$ has the structure (7.4). We are able to
efficiently solve system (7.8) by Fast Diagonalization (FD) method. The
starting point, is the setup of the preconditioner, that is

  * 1.

for the pencils $({\widehat{\mathbf{J}}_{l}},\widehat{\mathbf{M}}_{l})$, with
$l=1,\dots,d$ we have the factorization (7.5).

  * 2.

for the pencil $({\widehat{\mathbf{L}}_{t}},\widehat{\mathbf{M}}_{t})$ we have
the factorization (7.7).

Then, by defining
$\mathbf{V}_{s}:=\mathbf{V}_{d}\otimes\dots\otimes\mathbf{V}_{1}$ and
$\mathbf{\Lambda}_{s}:=\sum_{l=1}^{d}\mathbb{I}_{n_{s,d}}\otimes\dots\otimes\mathbb{I}_{n_{s,l+1}}\otimes\mathbf{\Lambda}_{l}\otimes\mathbb{I}_{n_{s,l-1}}\otimes\dots\otimes\mathbb{I}_{n_{s,1}}$,
we have for the matrix $\widehat{\mathbf{B}}$ the factorization

$\widehat{\mathbf{B}}=\left(\mathbf{V}_{t}^{T}\otimes\mathbf{V}_{s}^{T}\right)^{-1}\left({\gamma}^{2}\mathbf{\Lambda}_{t}\otimes\mathbb{I}_{N_{s}}+{\nu}^{2}{\mathbb{I}}_{n_{t}}\otimes{\mathbf{\Lambda}}_{s}\right)\left(\mathbf{V}_{t}\otimes\mathbf{V}_{s}\right)^{-1}.$
(7.9)

Notice that, the second factor in (7.9) is diagonal. Therefore, the solution
of (7.8) can be computed by the following algorithm.

Algorithm 3 Fast Diagonalization

1:Compute the factorizations (7.5) and (7.7).

2:Compute
$\widetilde{\mathbf{q}}=(\mathbf{V}_{t}^{T}\otimes\mathbf{V}_{s}^{T})\mathbf{r}$.

3:Compute
$\widetilde{\mathbf{s}}=\left({\gamma}^{2}\mathbf{\Lambda}_{t}\otimes\mathbb{I}_{N_{s}}+{\nu}^{2}{\mathbb{I}}_{n_{t}}\otimes{\mathbf{\Lambda}}_{s}\right)^{-1}\widetilde{\mathbf{q}}.$

4:Compute $\mathbf{s}=(\mathbf{V}_{t}\otimes\mathbf{V}_{s})\
\widetilde{\mathbf{s}}.$

###  7.4 Preconditioner robustness: partial inclusion of the geometry

The spectral estimate in Theorem 5 show the dependence on $\boldsymbol{G}$,
that is, the geometry parametrization affects the performance of our
preconditioner (7.3), as it is confirmed by the numerical tests in Section 10.
In this section, we present a strategy to partially incorporate
$\boldsymbol{G}$ in the preconditioner, without increasing its computational
cost. The same idea has been used in [36] for the stokes problem.

Let us split the bilinear form $\mathcal{B}(\cdot,\cdot)$ as

$\mathcal{B}(v_{h},w_{h})=\mathcal{B}_{t}(v_{h},w_{h})+\mathcal{B}_{s}(v_{h},w_{h})-\mathcal{R}(v_{h},w_{h}),\quad\forall
v_{h},w_{h}\in\mathcal{V}_{h},$

where

$\begin{split}\mathcal{B}_{t}(v_{h},w_{h}):=\gamma^{2}\int_{0}^{T}\int_{\Omega}\partial_{t}v_{h}\partial_{t}w_{h}&\mathrm{d}\Omega\,\mathrm{dt},\quad\mathcal{B}_{s}(v_{h},w_{h}):=\nu^{2}\int_{0}^{T}\int_{\Omega}\Delta
v_{h}\Delta w_{h}\mathrm{d}\Omega\,\mathrm{dt},\\\
\mathcal{R}(v_{h},w_{h}):=\gamma\nu\int_{0}^{T}\int_{\Omega}&\partial_{t}v_{h}\Delta
w_{h}+\Delta v_{h}\partial_{t}w_{h}\mathrm{d}\Omega\,\mathrm{dt}.\end{split}$

Using that $v_{h}=\widehat{v}_{h}\circ\boldsymbol{G}^{-1}$,
$w_{h}=\widehat{w}_{h}\circ\boldsymbol{G}^{-1}$ and for $i=1,\dots,d,$

$\dfrac{\partial^{2}v_{h}}{\partial
x_{i}^{2}}=\sum_{j,k=1}^{d}\dfrac{\partial^{2}\widehat{v}_{h}\circ\boldsymbol{G}^{-1}}{\partial\eta_{j}\partial\eta_{k}}[J_{\boldsymbol{F}}^{-1}]_{k,i}[J_{\boldsymbol{F}}^{-1}]_{j,i}+\sum_{j=1}^{d}\dfrac{\partial\widehat{v}_{h}\circ\boldsymbol{G}^{-1}}{\partial\eta_{j}}\dfrac{\partial[J_{\boldsymbol{F}}^{-1}]_{j,i}}{\partial\eta_{i}},$

we can write $\mathcal{B}_{t}$ and $\mathcal{B}_{s}$ as

$\mathcal{B}_{t}(v_{h},w_{h})=\int_{0}^{1}\int_{\widehat{\Omega}}c_{d+1}\partial_{\tau}\widehat{v}_{h}\partial_{\tau}\widehat{w}_{h}\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau,\quad\mathcal{B}_{s}(v_{h},w_{h})=\mathcal{B}_{s,1}(\widehat{v}_{h},\widehat{w}_{h})+\mathcal{B}_{s,2}(\widehat{v}_{h},\widehat{w}_{h}),$

with

$\begin{split}\mathcal{B}_{s,1}(\widehat{v}_{h},\widehat{w}_{h})&:=\sum_{l=1}^{d}\int_{0}^{1}\int_{\widehat{\Omega}}c_{l}\dfrac{\partial^{2}\widehat{v}_{h}}{\partial\eta_{l}^{2}}\dfrac{\partial^{2}\widehat{w}_{h}}{\partial\eta_{l}^{2}}\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau,\\\
\mathcal{B}_{s,2}(\widehat{v}_{h},\widehat{w}_{h})&:=\sum_{\begin{subarray}{c}r,s=1\\\
r\neq s\end{subarray}}^{d}\sum_{\begin{subarray}{c}j,k=1\\\ j\neq
k\end{subarray}}^{d}\int_{0}^{1}\int_{\widehat{\Omega}}g_{rsjk}^{1}\dfrac{\partial^{2}\widehat{v}_{h}}{\partial\eta_{k}\partial\eta_{j}}\dfrac{\partial^{2}\widehat{w}_{h}}{\partial\eta_{r}\partial\eta_{s}}\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau+\\\
&+\sum_{j,k=1}^{d}\int_{0}^{1}\int_{\widehat{\Omega}}g^{2}_{kj}\dfrac{\partial\widehat{v}_{h}}{\partial\eta_{k}}\dfrac{\partial\widehat{w}_{h}}{\partial\eta_{j}}\mathrm{d}\widehat{\Omega}\mathrm{d}\tau+\\\
&+\sum_{r=1}^{d}\sum_{j,k=1}^{d}\int_{0}^{1}\int_{\widehat{\Omega}}g_{rjk}^{3}\dfrac{\partial^{2}\widehat{v}_{h}}{\partial\eta_{k}\partial\eta_{j}}\dfrac{\partial\widehat{w}_{h}}{\partial\eta_{r}}+\dfrac{\partial\widehat{v}_{h}}{\partial\eta_{r}}\dfrac{\partial^{2}\widehat{w}_{h}}{\partial\eta_{k}\partial\eta_{j}}\mathrm{d}\widehat{\Omega}\,\mathrm{d}\tau,\end{split}$

and where we have defined

$c_{l}:=\nu^{2}\left(\|[J_{\boldsymbol{F}}^{-1}]_{\cdot,l}\|_{2}\right)^{4}|\text{det}(J_{\boldsymbol{F}})|T,\quad\text{for
}l=1,\dots,d,\quad c_{d+1}:=\gamma^{2}|\text{det}(J_{\boldsymbol{F}})|T^{-1},$

while $g^{1}_{rsjk},g^{2}_{j,k},g^{3}_{rjk}$ are functions that depend on the
parameterization $\boldsymbol{G}$. The preconditioner will be based on an
approximation of $\mathcal{B}_{t},\mathcal{B}_{s,1}$ only. In particular we
approximate $c_{l}$ for $l=1,\dots,d+1$ as

$\displaystyle
c_{l}(\boldsymbol{\eta},\tau)\approx{\mu_{1}(\eta_{1})\dots\mu_{l-1}(\eta_{l-1})\omega_{l}(\eta_{l})\mu_{l+1}(\eta_{l+1})\dots\mu_{d}(\eta_{d})\mu_{d+1}(\tau)},$
(7.10a) $\displaystyle
c_{d+1}(\boldsymbol{\eta},\tau)\approx{\mu_{1}(\eta_{1})\dots\mu_{d}(\eta_{d})\omega_{d+1}(\tau).}$
(7.10b)

The functions $c_{l}$ in (7.10) are first interpolated by constants in each
element and then the construction of the univariate factors $\mu_{l}$ and
$\omega_{l}$ is performed by the separation of variables algorithm detailed in
[21, Appendix C]. The resulting computational cost is proportional to the
number of elements, which for smooth splines i roughly equal to $N_{dof}$,
independent of the degrees $p_{s}$ and $p_{t}$, and therefore negligible in
the whole iterative strategy.

This first step leads to a matrix of this form

$\widetilde{\mathbf{B}}=\widehat{\mathbf{L}}_{t}^{\boldsymbol{G}}\otimes\widehat{\mathbf{M}}_{s}^{\boldsymbol{G}}+\widehat{\mathbf{M}}_{t}^{\boldsymbol{G}}\otimes\widetilde{\mathbf{J}}_{s}^{\boldsymbol{G}}$

where, with the notation of the basis functions in Section 2.2, for
$i,j=1,\dots,n_{t}$

$[\widehat{\mathbf{L}}_{t}]_{i,j}^{\boldsymbol{G}}:=\int_{0}^{1}\omega_{d+1}(\tau)\widehat{b}^{\prime}_{i,p_{t}}(\tau)\widehat{b}^{\prime}_{j,p_{t}}(\tau)\mathrm{d}\tau,\quad[\widehat{\mathbf{M}}_{t}]_{i,j}^{\boldsymbol{G}}:=\int_{0}^{1}\mu_{d+1}(\tau)\widehat{b}_{i,p_{t}}(\tau)\widehat{b}_{j,p_{t}}(\tau)\mathrm{d}\tau,$

while

$\widetilde{\mathbf{J}}_{s}^{\boldsymbol{G}}=\sum_{l=1}^{d}\widehat{\mathbf{M}}_{d}^{\boldsymbol{G}}\otimes\dots\otimes\widehat{\mathbf{M}}_{l+1}^{\boldsymbol{G}}\otimes\widehat{\mathbf{J}}_{l}^{\boldsymbol{G}}\otimes\widehat{\mathbf{M}}_{l-1}^{\boldsymbol{G}}\otimes\dots\otimes\widehat{\mathbf{M}}_{1}^{\boldsymbol{G}},\quad\widehat{\mathbf{M}}_{s}^{\boldsymbol{G}}=\widehat{\mathbf{M}}_{d}^{\boldsymbol{G}}\otimes\dots\otimes\widehat{\mathbf{M}}_{1}^{\boldsymbol{G}},$

and for $l=1,\dots,d$ and for $i,j=1,\dots,n_{s,l}$

$[\widehat{\mathbf{J}}_{l}]_{i,j}^{\boldsymbol{G}}:=\int_{0}^{1}\omega_{l}(\eta_{l})\widehat{b}^{\prime\prime}_{i,p_{s}}(\eta_{l})\widehat{b}^{\prime\prime}_{j,p_{s}}(\eta_{l})\mathrm{d}\eta_{l},\quad[\widehat{\mathbf{M}}_{l}]_{i,j}^{\boldsymbol{G}}:=\int_{0}^{1}\mu_{l}(\eta_{l})\widehat{b}_{i,p_{s}}(\eta_{l})\widehat{b}_{j,p_{s}}(\eta_{l})\mathrm{d}\eta_{l}.$

The matrix $\widetilde{\mathbf{B}}$ maintains the Kronecker structure of (7.3)
and Algorithm 3 can still be used to compute its application.

Finally, we apply a diagonal scaling and define the preconditioner as

$\widehat{\mathbf{B}}^{\boldsymbol{G}}:=\mathbf{D}^{1/2}\widetilde{\mathbf{B}}\mathbf{D}^{1/2}$
(7.11)

where $\mathbf{D}$ is the diagonal matrix whose entries are
$[\mathbf{D}]_{i,i}:=\dfrac{[\mathbf{B}]_{i,i}}{[\widehat{\mathbf{B}}^{\boldsymbol{G}}]_{i,i}}$
for $i=1,\dots,N_{dof}$.

######  Remark 3.

For the model problem considered, the approximation of the geometry
parametrization in the time direction is trivial. Notice that the coefficients
in (7.10) do not depend on $\tau$. Indeed, in our case it holds

$\mathbf{L}_{t}=\dfrac{1}{T}\widehat{\mathbf{L}}_{t},\quad\mathbf{M}_{t}={T}\widehat{\mathbf{M}}_{t},$

and hence we could set explicitly
$\widehat{\mathbf{L}}_{t}^{\boldsymbol{G}}=\mathbf{L}_{t}$ and
$\widehat{\mathbf{M}}_{t}^{\boldsymbol{G}}=\mathbf{M}_{t}$, which is exact.
However, we want to present the more general approximating strategy above
which could be used also when the spatial geometry or equation’s coefficients
depend on time.

##  8 Continuous Least Square preconditioner

In this section, we propose a second preconditioner for problem (6.14), which
exploits the Kronecker structure of (6.15). The definition of the
preconditioner is

$\widehat{\mathbf{Q}}:=\gamma^{2}\mathbf{L}_{t}\otimes\widehat{\mathbf{M}}_{s}+\nu^{2}\mathbf{M}_{t}\otimes\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}^{-1}_{s}\widehat{\mathbf{L}}_{s}+\gamma\nu\mathbf{R}_{t}\otimes\widehat{\mathbf{L}}_{s},$
(8.1)

where $\widehat{\mathbf{L}}_{s}$ and $\widehat{\mathbf{M}}_{s}$ are defined in
(4.2), while $\mathbf{L}_{t},\mathbf{M_{t}}$ and $\mathbf{R}_{t}$ are defined
in (6.16a). We recall that, by Remark 3,
$\mathbf{L}_{t}=\dfrac{1}{T}\widehat{\mathbf{L}}_{t}$ and
$\mathbf{M}_{t}=T\widehat{\mathbf{M}}_{t}$. Therefore, the preconditioner
$\widehat{\mathbf{Q}}$ can be written as

$\widehat{\mathbf{Q}}=\dfrac{\gamma^{2}}{T}\widehat{\mathbf{L}}_{t}\otimes\widehat{\mathbf{M}}_{s}+\nu^{2}T\widehat{\mathbf{M}}_{t}\otimes\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}^{-1}_{s}\widehat{\mathbf{L}}_{s}+\gamma\nu\mathbf{R}_{t}\otimes\widehat{\mathbf{L}}_{s}.$
(8.2)

Notice that,instead of the matrix $\widetilde{\mathbf{J}}_{s}$ of the previous
section, we consider the matrix
$\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}^{-1}_{s}\widehat{\mathbf{L}}_{s}$
in the spacial factor of the middle term. This choice is motivated by the
spectral equivalence between the two matrices, see [46, Proposition 4.1]. We
investigate numerically the stability of such spectral equivalence, and report
in Figure 7 the eigenvalues of
$(\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}^{-1}_{s}\widehat{\mathbf{L}}_{s})^{-1}\mathbf{J}_{s}$,
which are clustered close to $1$, for different mesh sizes $h$ and different
polynomial degrees $p$. In conclusion the spectral equivalence is stable under
$h$-refinement and $p$-refinement.

Figure 7: Eigenvalues of
$(\widehat{\mathbf{L}}_{s}^{T}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s})^{-1}\widetilde{\mathbf{J}}_{s}$
for different degrees $p$ and number of elements $n_{el}$.

Finally, the efficient application of the preconditioner exploits the
Kronecker structure in (8.1), and again the implementation involves Fast
Diagonalization method, together with Sherman-Morrison formula.

###  8.1 Stable factorization in space and time

The space pencils $(\widehat{\mathbf{L}}_{l},\widehat{\mathbf{M}}_{l})$ for
$l=1,\dots,d$ admit the stable factorizations described in Section 4.1, that
is

$\mathbf{U}^{T}_{l}\widehat{\mathbf{L}}_{l}\mathbf{U}_{l}=\mathbf{\Lambda}_{l}\quad\text{
and
}\quad\mathbf{U}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}=\mathbb{I}_{n_{s,l}},$

where $\mathbb{I}_{n_{s,l}}$ denotes the identity matrix of dimension
$n_{s,l}\times n_{s,l}$.

Whereas, the time pencil $(\widehat{\mathbf{L}}_{t},\widehat{\mathbf{M}}_{t})$
admit the stable factorization described in Section 7.2, that is

$\mathbf{V}^{T}_{t}\widehat{\mathbf{L}}_{t}\mathbf{V}_{t}=\mathbf{\Lambda}_{t}\quad\text{
and
}\quad\mathbf{V}^{T}_{t}\widehat{\mathbf{M}}_{t}\mathbf{U}_{t}=\mathbb{I}_{n_{t}},$

where $\mathbb{I}_{n_{t}}$ denotes the identity matrix of dimension
$n_{t}\times n_{t}$.

###  8.2 Preconditioner application

The application of the preconditioner involves the solution of the linear
system

$\widehat{\mathbf{Q}}\mathbf{s}=\mathbf{r},$ (8.3)

where $\widehat{\mathbf{Q}}$ has the structure (8.2). We are able to
efficiently solve system (8.3) by the FD method and the Sherman-Morrison
formula. The starting points, that are involved in the setup of the
preconditioner, are the following ones:

  * 1.

for the pencils $({\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l})$ for
$l=1,\dots,d$ we have the factorizations (4.4);

  * 2.

for the pencil $(\widehat{\mathbf{L}}_{t},\widehat{\mathbf{M}}_{t})$ we have
the factorization (7.7).

Then, define $\mathbf{U}_{s}:=\mathbf{U}_{d}\otimes\dots\otimes\mathbf{U}_{1}$
and
$\mathbf{\Lambda}_{s}:=\sum_{l=1}^{d}\mathbb{I}_{n_{s,d}}\otimes\dots\otimes\mathbb{I}_{n_{s,l+1}}\otimes\mathbf{\Lambda}_{l}\otimes\mathbb{I}_{n_{s,l-1}}\otimes\dots\otimes\mathbb{I}_{n_{s,1}}$.
Notice that
$\widehat{\mathbf{M}}_{s}^{-1}={\mathbf{U}}_{s}{\mathbf{U}}_{s}^{T}$,
therefore the matrix
$\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}$
admits the stable factorization (5.7), that is

$\mathbf{U}^{T}_{s}\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}\mathbf{U}_{s}=\mathbf{\Lambda}_{s}^{2}.$

The preconditioner $\widehat{\mathbf{Q}}$ admits the following factorization

$\widehat{\mathbf{Q}}=\left(\mathbb{I}_{n_{t}}\otimes\mathbf{U}_{s}^{T}\right)^{-1}\\\
\left(\left({\mathbf{V}}_{t}^{T}\otimes\mathbb{I}_{N_{s}}\right)^{-1}\left(\dfrac{\gamma^{2}}{T}\mathbf{\Lambda}_{t}\otimes\mathbb{I}_{N_{s}}+\nu^{2}T\mathbb{I}_{t}\otimes\mathbf{\Lambda}_{s}^{2}\right)\left({\mathbf{V}}_{t}\otimes\mathbb{I}_{N_{s}}\right)^{-1}+\gamma\nu\mathbf{R}_{t}\otimes\mathbf{\Lambda}_{s}\right)\\\
\left(\mathbb{I}_{n_{t}}\otimes\mathbf{U}_{s}\right)^{-1}.$ (8.4)

Note that the second factor in (8.4) is sum of two matrices, the second one
being

$\gamma\nu\mathbf{R}_{t}\otimes\mathbf{\Lambda}_{s}$ (8.5)

where $\mathbf{R}_{t}$ is a rank 1 matrix that can be written as
$\mathbf{R}_{t}=\boldsymbol{e}_{n_{t}}\boldsymbol{e}_{n_{t}}^{T}$, with
$\boldsymbol{e}_{n_{t}}$ being the last element of the canonical basis of
$\mathbb{R}^{n_{t}}$. Here it is possible to introduce a vector
$\boldsymbol{v}_{n_{t}}\in\mathbb{R}^{n_{t}}$ such that

$\boldsymbol{e}_{n_{t}}\boldsymbol{e}_{n_{t}}^{T}={\mathbf{V}}_{t}^{-T}\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}{\mathbf{V}}_{t}^{-1},$

by defining
$\boldsymbol{v}_{n_{t}}:={\mathbf{V}}_{t}^{T}\boldsymbol{e}_{n_{t}}$ that is
the last column of ${\mathbf{V}}_{t}^{T}$. Thus equation (8.4) can be
rewritten as

$\widehat{\mathbf{Q}}=\left({\mathbf{V}}_{t}^{T}\otimes\mathbf{U}_{s}^{T}\right)^{-1}(\mathbf{H}+\mathbf{K})\left({\mathbf{V}}_{t}\otimes\mathbf{U}_{s}\right)^{-1}$
(8.6)

where we have introduced

$\mathbf{H}:=\left(\dfrac{\gamma^{2}}{T}\mathbf{\Lambda}_{t}\otimes\mathbb{I}_{N_{s}}+\nu^{2}T\mathbb{I}_{t}\otimes\mathbf{\Lambda}_{s}^{2}\right)$

for the full diagonal term, and

$\mathbf{K}:=\gamma\nu\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}\otimes\mathbf{\Lambda}_{s}.$
(8.7)

Introducing $\mathbf{S}_{n_{t},N_{s}}$, the space and time shuffle matrix,
such that
$\mathbf{A}\otimes\mathbf{B}=\mathbf{S}_{n_{t},N_{s}}(\mathbf{B}\otimes\mathbf{A})\mathbf{S}_{n_{t},N_{s}}^{T}$
for all $\mathbf{A}\in\mathbb{R}^{n_{t}\times n_{t}}$ and
$\mathbf{B}\in\mathbb{R}^{N_{s}\times N_{s}}$. Under reshuffling with
$\mathbf{S}_{n_{t},N_{s}}$, we have the following block diagonal structure

$\begin{split}\mathbf{H}=\mathbf{S}_{n_{t},N_{s}}\begin{bmatrix}\mathbf{H}_{1}&&\\\
&\ddots&\\\
&&\mathbf{H}_{N_{s}}\end{bmatrix}\mathbf{S}_{n_{t},N_{s}}^{T},\quad\text{and}\quad\mathbf{K}=\mathbf{S}_{n_{t},N_{s}}\begin{bmatrix}\mathbf{K}_{1}&&\\\
&\ddots&\\\
&&\mathbf{K}_{N_{s}}\end{bmatrix}\mathbf{S}_{n_{t},N_{s}}^{T},\end{split}$

where, for $i=1,\dots,N_{s}$, the matrices $\mathbf{H}_{i}$ are diagonal
matrices defined as
$\mathbf{H}_{i}:=\dfrac{\gamma^{2}}{T}\mathbf{\Lambda}_{t}+{\nu^{2}T}[\mathbf{\Lambda}_{s}^{2}]_{i,i}\otimes\mathbb{I}_{n_{t}}$
, while
$\mathbf{K}_{i}=\gamma\nu[\mathbf{\Lambda}_{s}]_{i,i}\otimes\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}$.

In order to invert $\mathbf{H}+\mathbf{K}$, it is now sufficient to reshuffle
the data, and invert the following independent $N_{s}$ problems of size
$n_{t}\times n_{t}$:

$\left(\mathbf{H}_{i}+\mathbf{K}_{i}\right)\mathbf{x}_{i}=\mathbf{y}_{i}\quad\text{
for }i=1,\dots,N_{s}.$

Notice that, each $\mathbf{K}_{i}$ is a rank 1 perturbation in the above
systems, therefore the Sherman-Morrison formula gives

$\displaystyle\mathbf{x}_{i}=\mathbf{H}_{i}^{-1}\mathbf{y}_{i}-\gamma\nu[\Lambda_{s}]_{i,i}\dfrac{\mathbf{H}_{i}^{-1}\boldsymbol{v}_{n_{t}}\boldsymbol{v}_{n_{t}}^{T}\mathbf{H}_{i}^{-1}\mathbf{y}_{i}}{1+\gamma\nu[\Lambda_{s}]_{i,i}\boldsymbol{v}_{n_{t}}^{T}\mathbf{H}_{i}^{-1}\boldsymbol{e}_{n_{t}}}\quad\text{
for }i=1,\dots,N_{s}.$ (8.8)

Summarizing, the solution of (8.3) can be computed by Algorithm 2, with the
proper notation for the involved matrices.

##  9 Computational cost and memory requirement

In this section we discuss the computational costs and memory requirements in
the implementation of Algorithms 1,2 and 3. First, notice that the matrix
$\mathbf{A}$ in (3.6) is neither positive definite nor symmetric and we choose
GMRES as linear solver for the system (3.5). Whereas, the matrix
$\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{A}$ in (5.2) and $\mathbf{B}$ in (6.15)
are symmetric positive definite, therefore we choose Conjugate Gradients (CG)
as linear solver for solving (5.1) and (6.14). Clearly, the computational cost
of each iteration of the CG solver, or GMRES solver, depends on both the
preconditioner setup and application cost, and on the residual computation
cost.

We assume for simplicity that, for each univariate direction $l=1,\dots,d$,
the space matrices have dimension $n_{s}\times n_{s}$, while the time matrices
involved in the preconditioners, have dimension $n_{t}\times n_{t}$. Thus the
total number of degrees-of-freedom is $N_{dof}=N_{s}n_{t}=n_{s}^{d}n_{t}$.

We remark that the setup of the preconditioners has to be performed only once,
since the matrices involved do not change during the iterative procedure.

###  9.1 Setup and application cost of Algorithm 1

The setup of $\widehat{\mathbf{A}}$ and
$\widehat{\mathbf{A}}^{\boldsymbol{G}}$ includes the operations performed in
Step 1 of Algorithm 1, i.e. $d$ spatial eigendecompositions, that have a total
cost of $O(dn_{s}^{3})$ FLOPs, and the factorization of the time matrices. The
computational cost of the latter, that is the sum of the cost of the
eigendecomposition (4.12) and of the cost to compute the solution $\mathbf{v}$
of the linear system (4.14), yields a cost of $O(n_{t}^{3})$ FLOPs. Then, the
total cost of the spatial and time factorizations is $O(dn_{s}^{3}+n_{t}^{3})$
FLOPs. Note that, if $n_{t}=O(n_{s})$, this cost is optimal for $d=2$ and
negligible for $d=3$. The setup cost of
$\widehat{\mathbf{A}}^{\boldsymbol{G}}$ includes also the the construction of
the diagonal matrix $\mathbf{D}$, that has a negligible cost, and the
computation of the $2d$ approximations $\varphi_{1},\dots,\varphi_{d}$ and
$\Phi_{1},\dots,\Phi_{d}$ in (4.20), whose cost is negligible too, as
mentioned in Section 4.4.

The application of the preconditioner is performed by Steps 2-4 of Algorithm
1. Exploiting (2.5), Step 2 and Step 4 costs
$4(dn_{s}^{d+1}n_{t}+n_{t}^{2}n_{s}^{d})=4N_{dof}(dn_{s}+n_{t})$ FLOPs. The
use of the block LU decomposition (4.19) makes the cost for Step 3 equal to
$O(N_{dof})$ FLOPs.

In conclusion, the total cost of Algorithm 1 is
$4N_{dof}(dn_{s}+n_{t})+O(N_{dof})$ FLOPs. The non-optimal dominant cost of
Step 2 and Step 4 is determined by the dense matrix-matrix products. However,
these operations are usually implemented on modern computers in a very
efficient way. For this reason, in our numerical tests, the overall serial
computational time grows almost as $O(N_{dof})$, as reported in [22, Figure
5].

###  9.2 Setup and application cost of Algorithm 2

We apply Algorithm 2 both for the preconditioner $\widehat{\mathbf{P}}$ of
Section 5 and for the preconditioner $\widehat{\mathbf{Q}}$ of Section 8. We
discuss in detail the first strategy, since the second is analogous.

First notice that, for the discrete least square formulation of Section 5, we
have to consider also the setup of the right hand side
$\mathbf{g}=\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{f}$, which requires to invert
the mass matrix $\mathbf{M}=\mathbf{M}_{t}\otimes\mathbf{M}_{s}$. Here we
follow the preconditioned iterative technique proposed in [41], whose
computational cost of each iteration for inverting the mass is $O(pN_{dof})$
FLOPs, assuming $p=p_{s}\approx p_{t}$. We recall no inversion of the mass is
need for the strategy presented in 8.

The setup of $\widehat{\mathbf{P}}$ includes the operations performed in Step
1 of Algorithm 2, i.e. $d$ spatial eigendecompositions, that have a total cost
of $O(dn_{s}^{3})$ FLOPs, and one time eigendecomposition, that have a total
cost of $O(n_{t}^{3})$ FLOPs. Then, the total cost of the spatial and time
factorizations is $O(dn_{s}^{3}+n_{t}^{3})$ FLOPs. Again, if $n_{t}=O(n_{s})$,
this cost is optimal for $d=2$ and negligible for $d=3$.

The application of the preconditioner is performed by Steps 2-6 of Algorithm
2. Exploiting again (2.5), Step 2 and Step 6 costs $4N_{dof}(dn_{s}+n_{t})$
FLOPs. The cost of the reshuffling in Step 3 and Step 5 is negligible. Step 4
requires to solve $N_{s}$ independent problems of size $n_{t}\times n_{t}$.
The Sherman-Morrison formula together with the time decomposition, allow to
compute each solution by first inverting the diagonal matrices
$\mathbf{H}_{i}$ for $i=1,\dots,N_{s}$, and then computing the correction term
in (5.11). Both of these steps cost $O(n_{t}N_{s})=O(N_{dof})$ FLOPs. In
conclusion, the total cost of Algorithm 2 is
$4N_{dof}O(dn_{s}+n_{t})+O(N_{dof})$.

###  9.3 Setup and application cost of Algorithm 3

The setup of the preconditioners $\widehat{\mathbf{B}}$ and
$\widehat{\mathbf{B}}^{\boldsymbol{G}}$ include again the eigendecomposition
of the pencils, that is Step 1 of Algorithm 3. Therefore, the cost of the
factorizations is $O(dn_{s}^{3}+n_{t}^{3})$ FLOPs. We remark this cost is
optimal for $d=2$ and negligible for $d=3$, provided that $n_{t}\approx
n_{s}$. The setup cost of $\widehat{\mathbf{B}}^{\boldsymbol{G}}$ includes
also the the construction of the diagonal matrix $\mathbf{D}$, that has a
negligible cost, and the computation of the $2(d+1)$ univariate approximations
$\mu_{1},\dots,\mu_{d+1}$ and $\omega_{1},\dots,\omega_{d+1}$, that are used
to incorporate the information of the geometry into the preconditioner. As
mentioned above, this cost is negligible.

The application of the preconditioner is performed by Steps 2-4 of Algorithm
(3). Exploiting (2.5), Step 2 and Step 4 costs $4N_{dof}(dn_{s}+n_{t})$ FLOPs,
while Step 3 has an optimal cost, as it requires $O(N_{dof})$ FLOPs. Thus, the
total cost of Algorithm 3 is $4N_{dof}(dn_{s}+n_{t})+O(N_{dof})$ FLOPs. We
remark that the non-optimal dominant cost is given by the dense matrix-matrix
products of Step 2 and Step 4, which, however, are usually implemented on
modern computers in a high efficient way, as they are BLAS level 3 operations.

###  9.4 Computational cost of the residuals

The other dominant computational cost in a CG or GMRES iteration is the cost
of the residual computation. In Algorithm 1, this involves the multiplication
of the matrix $\mathbf{A}$ with a vector. This multiplication is done by
exploiting the special structure (3.6), that allows a matrix-free approach and
the use of formula (2.5). As mentioned above, the computational cost of a
single matrix-vector product is $O(N_{dof}p^{d})$ FLOPs, if we assume
$p=p_{s}\approx p_{t}$.

In Algorithm 2, the multiplication by
$\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{A}$ is done by exploiting the kronecker
structure in (3.6) together with the above mentioned iterative technique for
the inversion of the mass matrix. Notice that, this allows a matrix free
approach and uses the formula (2.5). In particular we do not need to compute
and to store the whole matrix $\mathbf{A}^{T}\mathbf{M}^{-1}\mathbf{A}$, but
only the time and spatial factors of the matrix $\mathbf{A}$ are enough. The
time matrices $\mathbf{M}_{t}$ and $\mathbf{W}_{t}$ are banded with a band of
width $2p_{t}+1$ and the spatial matrices $\mathbf{K}_{s}$ and
$\mathbf{M}_{s}$ have roughly $N_{s}(2p_{s}+1)^{d}$ nonzero entries. Assuming
$p=p_{s}\approx p_{t}$, the computational cost of a single matrix-vector
product is given by:

  * 1.

the multiplication by $\mathbf{A}$, whose cost in FLOPs is

$\displaystyle
O((2p_{s}+1)^{d}N_{s}n_{t}+(2p_{t}+1)n_{t}N_{s}]=O(p^{d}N_{dof})$

  * 2.

the inversion of the mass $\mathbf{M}$, which costs $O(pN_{dof})$ FLOPs;

  * 3.

the multiplication by $\mathbf{A}^{T}$, that costs again $O(p^{d}N_{dof})$
FLOPs.

The overall cost is thus $O(p^{d}N_{dof})$ FLOPs. Notice that the number of
iterations required in the CG solver for inverting the mass $\mathbf{M}_{s}$
over 3D objects, with a tolerance of $10^{-8}$, may be around 5, see [41,
Tables 1]. This factor, together with the possible higher number of iterations
required by the iterative strategy preconditioned by $\widehat{\mathbf{P}}$,
may cause lack of performance for this preconditioning strategy w.r.t the
previous one.

In Algorithm 3, the residual computation consists in the multiplication
between ${\mathbf{B}}$ and a vector. This multiplication can be computed by
exploiting the special structure (6.15) and the formula (2.5). Again, we do
not need to compute and store the whole matrix $\mathbf{B}$ but only its
factors
$\mathbf{L}_{t},\mathbf{M}_{t},\mathbf{R}_{t},\mathbf{J}_{s},\mathbf{L}_{s}$and
$\mathbf{M}_{s}$. With this matrix-free approach, noting again that the time
matrices $\mathbf{L}_{t},\mathbf{M}_{t},\mathbf{R}_{t}$ are banded matrices
with a band of width $2p_{t}+1$ and the spatial matrices
$\mathbf{J}_{s},\mathbf{L}_{s},\mathbf{M}_{s}$ have a number of non-zeros per
row approximately equal to $(2p_{s}+1)^{d}$, the computational cost of a
single matrix-vector product is $O(N_{dof}p^{d})$ FLOPs, if we assume
$p=p_{s}\approx p_{t}$.

The numerical experiments reported in [22, Table 5] and [21, Table 3] show
that the dominant cost in the iterative solver is represented by the residual
computation. This is a typical behavior of the FD-based preconditioning
strategies, see [22, 21, 36, 26].

###  9.5 Memory requirements

We now investigate the memory consumption of the preconditioning strategies
proposed, giving the details for the preconditioner $\widehat{\mathbf{A}}$,
since the other cases are analogous. For the preconditioner, we have to store
the eigenvector spatial matrices, $\mathbf{U}_{1},\dots,\mathbf{U}_{d}$, the
time matrix $\mathbf{U}_{t}$ and the block-arrowhead matrix (4.18). The memory
required is roughly

$n_{t}^{2}+dn_{s}^{2}+2N_{dof}.$

This extends analogously to the preconditioners proposed in Sections 5, 7 and
8.

For the system matrix $\mathbf{A}$, we have to store the time factors
$\mathbf{M}_{t}$ and $\mathbf{W}_{t}$ and the spatial factors $\mathbf{M}_{s}$
and $\mathbf{K}_{s}$. Thus the memory required is roughly

$2(2p_{t}+1)n_{t}+2(2p_{s}+1)^{d}N_{s}\approx
4p_{t}n_{t}+2^{d+1}p_{s}^{d}N_{s}.$

The discrete least squares approach of Section 5, further requires the
application of $\mathbf{M}^{-1}$. Using the approach of [41], requires a
further spatial eigendecomposition of the mass matrix $\mathbf{M}$, thus the
memory required is again $n_{t}^{2}+dn_{s}^{2}+2N_{dof}$, plus a diagonal
scaling $\mathbf{D}$ whose memory consumption is $N_{dof}$.

These numbers show that memory-wise our space-time strategies are very
appealing when compared to other approaches, even when space and time
variables are discretized separately, e.g., with finite differences in time or
other time-stepping schemes. For example if we assume $d=3$, $p_{t}\approx
p_{s}=p$ and $n_{t}^{2}\leq Cp^{3}N_{s}$, then the total memory consumption is
$O(p^{3}N_{s}+N_{dof})$, that is equal to the sum of the memory needed to
store the Galerkin matrices associated to spatial variables and the memory
needed to store the solution of the problem.

##  10 Numerical Results

In this section we present the numerical experiments that assess the
performance of the preconditioners. As regards the orders of convergence of
the discretizations of Section 3 and Section 6, we refer respectively to [22]
and [21].

We consider only sequential executions and we force the use of a single
computational thread in a Intel Core i7-5820K processor, running at 3.30 GHz
and with 64 GB of RAM.

The tests are performed with Matlab R2023a and GeoPDEs toolbox [47]. We use
the eig Matlab function to compute the generalized eigendecompositions present
in Step 1 of Algorithms 1,2 and 3 , while Tensorlab toolbox [48] is employed
to perform the multiplications with Kronecker matrices. The solution of the
linear system (4.14) is performed by Matlab direct solver (backslash operator
“$\backslash$”). The linear system (3.5) is solved by GMRES (Matlab routine
gmres), while the linear systems 5.1 and 6.14 are solved with CG (Matlab
routine pcg). The tolerance in the iterative solvers is set equal to $10^{-8}$
and the null vector is the initial guess in all tests. We remark that GMRES
computes and stores a full orthonormal basis for the Krylov space, and this
might be unfeasible if the number of iterations is too large. This issue could
be addressed by switching to a different solver for nonsymmetric systems, like
e.g. BiCGStab, or using the restarted version of GMRES.

We use the same mesh-size in space and in time $h_{s}=h_{t}=:h$, and use
splines of maximal continuity and same degree in space and in time
$p_{t}=p_{s}=:p$. For the sake of simplicity, we also consider uniform knot
vectors, and denote the number of elements in each parametric direction by
$n_{el}:=\frac{1}{h}$.

To assess the performance of our preconditioning strategies, we set $T=1$ and
we focus on two 3D spatial domains $\Omega\subset\mathbb{R}^{3}$, represented
in Figure LABEL:heat_fig:cube and Figure LABEL:heat_fig:rev-quarter: the cube
and the rotated quarter of annulus, respectively.

In out tables, the symbol $``\ast"$ denotes that the construction of the
matrix factors of $\mathbf{A}$ (see (3.6)) goes out of memory, while the
symbol $``\ast\ast"$ indicates that the dimension of the Krylov subspace is
too high and there is not enough memory to store all the GMRES iterates. We
remark that in all the tables the total solving time of the iterative
strategies includes also the setup time of the considered preconditioner.

(a)

(b)

Figure 8: Computational domains.

###  10.1 Performance of the preconditioners: cube domain

In the cube domain, we set homogeneous Dirichlet and zero initial boundary
conditions and we fix $f$ such that the exact solution is $u=\sin{(\pi
x)}\sin{(\pi y)}\sin{(\pi z)}\sin{(t)}$. Clearly, in this computational
domain, the preconditioners $\widehat{\mathbf{A}}$ and $\widehat{\mathbf{P}}$
are direct solvers. Moreover, since we are solving on the parametric domain,
we have $\widehat{\mathbf{A}}^{\boldsymbol{G}}=\widehat{\mathbf{A}}$ and
$\widehat{\mathbf{B}}^{\boldsymbol{G}}=\widehat{\mathbf{B}}$.

In view of this observations, we analyze the performance of
$\widehat{\mathbf{B}}$ and $\widehat{\mathbf{Q}}$. As a comparison we also
consider as preconditioner for CG the Incomplete Cholesky with zero fill-in
(IC(0)) factorization of $\mathbf{B}$, that is executed with Matlab routine
ichol. Table 7 report the number of iterations and the total solving time,
that includes the setup time of the preconditioner. The results for
$\widehat{\mathbf{B}}$ (upper table) and IC(0) (lower table) are taken from
[21, Table 1]. The matrix-vector products of CG are computed in a matrix-free
way using its time and spacial factors. Matrix $\mathbf{B}$ is assembled when
we want to use the IC(0) preconditioner. In any case, the assembly times are
never included in the reported times. The number of iterations obtained with
$\widehat{\mathbf{B}}$ and $\widehat{\mathbf{Q}}$ are stable with respect to
the polynomial degree $p:=p_{t}=p_{s}$. The performance of
$\widehat{\mathbf{B}}$ is stable also with respect to the number of elements
$n_{el}$, while $\widehat{\mathbf{Q}}$ seems to converge to a direct solver.
Even in the case when the number of iterations of $\widehat{\mathbf{B}}$ might
be larger than that of IC(0), the overall computational time is significantly
lower, up to two orders of magnitude for the problems considered. This is due
to the higher setup and application cost of the IC(0) preconditioner.

|  $\widehat{\mathbf{B}}$ \+ CG Iterations / Time
---|---
$n_{el}$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 |  09 / 0000.06 | 11 / 0000.07 | 11 / 000.18 | 11 / 000.28
16 | 11 / 0000.27 | 11 / 0000.69 | 12 / 001.80 | 12 / 003.80
32 | 12 / 0005.10 | 12 / 0013.17 | 12 / 027.31 | 12 / 052.95
64 | 13 / 0100.09 | 13 / 0227.93 | 13 / 458.86 | 13 / 924.44
128 | 13 / 2012.94 | 13 / 4235.96 | $\ast$ | $\ast$

|  $\widehat{\mathbf{Q}}$ \+ CG Iterations / Time
---|---
$n_{el}$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 | 3 / 000.03 | 5 / 000.08 | 5 / 000.19 | 4 / 000.26
16 | 3 / 000.13 | 5 / 000.78 | 5 / 002.18 | 3 / 002.61
32 | 2 / 001.79 | 4 / 007.17 | 3 / 012.70 | 2 / 019.60
64 | 2 / 032.36 | 4 / 113.51 | 3 / 186.92 | 2 / 356.48
128 | da fare… | da fare… | $\ast$ | $\ast$

|  IC(0) + CG Iterations / Time
---|---
$n_{el}$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 |  09 / 000.18 |  07 / 0001.69 |  06 / 014.04 |  06 / 0080.39
16 | 22 / 005.01 | 16 / 0045.54 | 12 / 355.99 | 10 / 1913.90
32 | 64 / 157.05 | $\ast$ | $\ast$ | $\ast$

Table 7: Cube domain. Performance of $\widehat{\mathbf{B}}$,
$\widehat{\mathbf{Q}}$ and IC(0).

###  10.2 Performance of the preconditioners: rotated quarter of annulus

The second computational domain $\Omega$ is a quarter of annulus with center
in the origin, internal radius 1 and external radius 2, rotated by $\pi/2$
along the axis $y=-1$. Boundary data and forcing function are set such that
the exact solution is
$u=-(x^{2}+y^{2}-1)(x^{2}+y^{2}-4)xy^{2}\sin{(z)}\sin{(t)}$.

We analyze the performance of $\widehat{\mathbf{A}}$,
$\widehat{\mathbf{A}}^{\boldsymbol{G}}$, $\widehat{\mathbf{P}}$,
$\widehat{\mathbf{B}}$, $\widehat{\mathbf{B}}^{\boldsymbol{G}}$ and
$\widehat{\mathbf{Q}}$. In the GMRES solver, the maximum dimension of the
Krylov subspace is set equal to 100 for both the preconditioners
$\widehat{\mathbf{A}}$ and $\widehat{\mathbf{A}}^{\boldsymbol{G}}$, up to
$n_{el}=64$. We are able to reach convergence and to perform the tests with
$\widehat{\mathbf{A}}^{\boldsymbol{G}}$, $n_{el}=128$ and $p=1,2,3$ by setting
the maximum Krylov subspace dimension equal to 25.

In Table 8 we first report the number of iterations and the total solving time
of GMRES preconditioned with $\widehat{\mathbf{A}}$ (upper table) and
$\widehat{\mathbf{A}}^{\boldsymbol{G}}$ (middle table), taken from [22, Table
4]. The non-trivial geometry clearly affects the performance of
$\widehat{\mathbf{A}}$, but, when we include some information on the
parametrization by using $\widehat{\mathbf{A}}^{\boldsymbol{G}}$, the number
of iterations is more than halved and it is stable w.r.t. $p$ and $n_{el}$.
Moreover, the computational times are one order of magnitude lower for the
highest degrees $p$ and numbers of elements $n_{el}$.

Next, in the lower table of Table 8 we report the results for the
preconditioner $\widehat{\mathbf{P}}$ obtained by solving the discrete least-
squares problem (5.1). Recall in this case the iterative solver is the
preconditioned conjugate gradient method, with tolerance $10^{-8}$ and initial
guess the null vector. In this case the number of iterations with respect to
the preconditioner $\widehat{\mathbf{A}}$ is more than doubled, although they
are stable with respect to the degrees and number of elements, we suggest to
use $\widehat{\mathbf{P}}$ as a last resort.

|  $\widehat{\mathbf{A}}$ Iterations / Time
---|---
$n_{el}$ | $p=1$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 | 34 / 000.20 | 37 / 000.21 | 42 / 0000.42 | 46 / 0000.63 | 50 / 0001.13
16 | 43 / 001.15 | 46 / 001.65 | 50 / 0003.42 | 54 / 0005.80 | 57 / 0011.87
32 | 50 / 022.75 | 53 / 031.10 | 57 / 0054.02 | 61 / 0096.06 | 64 / 0184.84
64 | 57 / 586.73 | 60 / 764.26 | 67 / 1254.81 | 67 / 1858.55 | 71 / 3188.51
128 | $\ast\ast$ | $\ast\ast$ | $\ast\ast$ | $\ast$ | $\ast$

|  $\widehat{\mathbf{A}}^{\boldsymbol{G}}$ Iterations / Time
---|---
$n_{el}$ | $p=1$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 | 11 / 0000.06 | 12 / 0000.09 | 12 / 0000.11 | 13 / 000.18 | 14 / 000.29
16 | 13 / 0000.26 | 14 / 0000.52 | 14 / 0001.18 | 14 / 001.44 | 15 / 003.85
32 | 15 / 0004.73 | 15 / 0006.76 | 15 / 0012.67 | 15 / 021.47 | 16 / 040.54
64 | 16 / 0107.24 | 16 / 0135.74 | 18 / 0249.27 | 16 / 370.31 | 17 / 695.44
128 | 17 / 2623.57 | 17 / 3105.76 | 17 / 5614.10 | $\ast$ | $\ast$

|  $\widehat{\mathbf{P}}$ Iterations / Time
---|---
$n_{el}$ | $p=1$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 |  093 / 0001.90 |  095 / 0003.15 |  098 / 0002.84 | 104 / 0006.28 | 107 / 00019.14
16 | 119 / 0005.60 | 118 / 0016.24 | 117 / 0038.84 | 119 / 00112.74 | 117 / 00232.96
32 | 129 / 0178.08 | 127 / 0285.19 | 125 / 0645.68 | 125 / 01711.45 | 124 / 04038.55
64 | 133 / 1607.33 | 130 / 3834.92 | 129 / 8981.30 | 144 / 19868.87 | 165 / 43930.63
128 | da fare … | da fare … | da fare … | $\ast$ | $\ast$

Table 8: Rotated quarter domain. Performance of $\widehat{\mathbf{A}}$,
$\widehat{\mathbf{A}}^{\boldsymbol{G}}$ and $\widehat{\mathbf{P}}$.

In Table 9, we report the results obtained in [21, Table 2], for
$\widehat{\mathbf{B}}$ (top section of the table) and
$\widehat{\mathbf{B}}^{\boldsymbol{G}}$ (middle section of the table), applied
to the same problem in the least-squares formulation framework. For the
preconditioner $\widehat{\mathbf{B}}$, the numbers of iterations have more
than doubled, while still remaining stable with respect to degrees and numbers
of elements. Finally, in the lower section of Table 9 we present the results
obtained with the preconditioner $\widehat{\mathbf{Q}}$, which is stable
w.r.t. the degrees and numbers of elements, and behaves like
$\widehat{\mathbf{B}}$, suggesting that including the information of the
geometry may boost its performance.

|  $\widehat{\mathbf{B}}$ \+ CG Iterations / Time
---|---
$n_{el}$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 | 107 / 00000.21 | 107 / 00000.48 | 114 / 0001.17 | 123 / 0002.73
16 | 126 / 00002.56 | 128 / 00006.90 | 133 / 0017.04 | 135 / 0035.17
32 | 142 / 00052.77 | 143 / 00132.24 | 148 / 0292.53 | 151 / 0572.84
64 | 153 / 01056.21 | 155 / 02415.23 | 156 / 4956.68 | 159 / 9906.33
128 | 164 / 22106.01 | 166 / 47539.02 | $\ast$ | $\ast$

|  $\widehat{\mathbf{B}}^{\boldsymbol{G}}$ \+ CG Iterations / Time
---|---
$n_{el}$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 | 24 / 0000.09 | 24 / 00000.13 | 26 / 0000.37 | 26 / 0000.60
16 | 35 / 0000.77 | 34 / 00001.96 | 33 / 0004.62 | 33 / 0009.35
32 | 42 / 0017.03 | 41 / 00039.57 | 40 / 0082.35 | 41 / 0161.73
64 | 46 / 0333.20 | 44 / 0716.03 | 49 / 1577.55 | 53 / 3384.08
128 | 48 / 6767.08 | 50 / 14814.09 | $\ast$ | $\ast$

|  $\widehat{\mathbf{Q}}$ \+ CG Iterations / Time
---|---
$n_{el}$ | $p=2$ | $p=3$ | $p=4$ | $p=5$
8 | 101 / 0000.35 | 111 / 0001.27 | 136 / 0003.02 | 172 / 00006.78
16 | 132 / 0003.08 | 140 / 0014.98 | 161 / 0041.14 | 175 / 00092.71
32 | 139 / 0072.11 | 151 / 0190.92 | 169 / 0444.41 | 186 / 00938.47