Robust preconditioner for MPETE. Piersanti, J. J. Lee, K.-A. Mardal, M. E.
Rognes, and T. Thompson

# Parameter robust preconditioning by congruence for multiple-network
poroelasticity ††thanks: Submitted to the editors . The work of J. J. Lee has
been supported by the European Research Council under the European Union’s
Seventh Framework Programme (FP7/2007-2013) ERC grant agreement 339643. The
work of K.-A. Mardal, T. Thompson and M. E. Rognes have been supported by the
Research Council of Norway under the FRINATEK Young Research Talents Programme
through project #250731/F20 (Waterscape). E. Piersanti is a doctoral fellow in
the Simula-UCSD-University of Oslo Research and PhD training (SUURPh) program,
an international collaboration in computational biology and medicine funded by
the Norwegian Ministry of Education and Research.

E. Piersanti Simula Research Laboratory, P. O. Box 134, 1325 Lysaker, Norway
<EMAIL_ADDRESS>J. J. Lee Department of Mathematics, Baylor University,
One Bear Place # 97328, Waco, Texas 78798, USA<EMAIL_ADDRESS>T.
Thompson Simula Research Laboratory, Norway. Currently at: Mathematical
Institute, Oxford University<EMAIL_ADDRESS>K.-A. Mardal
Department of Mathematics, University of Oslo, P. O. Box 1053 Blindern, 0316
Oslo, Norway and Simula Research Laboratory, P. O. Box 134, 1325 Lysaker,
Norway<EMAIL_ADDRESS>M. E. Rognes Simula Research Laboratory, P. O.
Box 134, 1325 Lysaker, Norway<EMAIL_ADDRESS>
###### Abstract

The mechanical behaviour of a poroelastic medium permeated by multiple
interacting fluid networks can be described by a system of time-dependent
partial differential equations known as the multiple-network poroelasticity
(MPET) equations or multi-porosity/multi-permeability systems. These equations
generalize Biot’s equations, which describe the mechanics of the one–network
case. The efficient numerical solution of the MPET equations is challenging,
in part due to the complexity of the system and in part due to the presence of
interacting parameter regimes. In this paper, we present a new strategy for
efficiently and robustly solving the MPET equations numerically. In
particular, we introduce a new approach to formulating finite element methods
and associated preconditioners for the MPET equations. The approach is based
on designing transformations of variables that simultaneously diagonalize (by
congruence) the equations’ key operators and subsequently constructing
parameter-robust block-diagonal preconditioners for the transformed system.
Our methodology is supported by theoretical considerations as well as by
numerical results.

##  1 Introduction

In this paper, we consider the preconditioned iterative solution of finite
element discretizations of the multiple-network poroelasticity (MPET)
equations. These equations traditionally originate in geomechanics where they
are also known under the term multi-porosity/multi-permeability systems [1].
The MPET equations generalize Biot’s equations [2] from the one network to the
multiple network case, and multi-compartment Darcy (MPT) equations [19] from a
porous (but rigid) to a poroelastic medium. Over the last decade, the MPT and
MPET equations have seen a surge of interest in biology and physiology; e.g.
to model perfusion in the heart [19, 14], cancer [25] brain [13] or liver [3],
or to model the interaction between elastic deformation and fluid flow and
transport in the brain [5, 6, 21, 28, 29, 30].

Concretely, the quasi–static MPET equations read as follows [1]: for a given
number of networks $J\in\mathbb{N}$, find the displacement $u$ and the network
pressures $p_{j}$ for $j=1,\dots,J$ such that

(1.1a)
$\displaystyle-\operatorname{div}(2\mu\varepsilon(u)+\lambda\operatorname{div}u\mathbb{I})+\operatorname{\textstyle\sum}_{j}\alpha_{j}\operatorname{\nabla}p_{j}$
$\displaystyle=f,$ (1.1b) $\displaystyle
s_{j}\dot{p}_{j}+\alpha_{j}\operatorname{div}\dot{u}-\operatorname{div}K_{j}\operatorname{\nabla}p_{j}+\operatorname{\textstyle\sum}_{i}\xi_{j\leftarrow
i}(p_{j}-p_{i})$ $\displaystyle=g_{j},$

where $u=u(x,t)$, $p_{j}=p_{j}(x,t)$ for $x\in\Omega\subset{\mathbb{R}}^{d}$
($d=1,2,3$), $t\in(0,T]$, and $\mathbb{I}$ is the $d\times d$ identity matrix.
Physically, the equations (1.1) describe a porous and elastic medium permeated
by a number of fluid networks under the assumptions that the solid matrix can
be modeled as isotropic and linearly elastic with Lamé constants $\mu>0$ and
$\lambda>0$, and the transfer between the networks is regulated by the
corresponding pressure differences with exchange coefficients
$\xi_{j\leftarrow i}\geq 0$. For each network $j$, we define the Biot-Willis
coefficient $\alpha_{j}\in(0,1]$ such that $\sum_{j}\alpha_{j}\leq 1$, the
storage coefficient $s_{j}>0$, and the hydraulic conductivity tensor
$K_{j}=\kappa_{j}/\nu_{j}>0$ with $\kappa_{j}$ and $\nu_{j}$ being the
permeability and fluid viscosity, respectively. Moreover,
$\operatorname{\nabla}$ denotes the column-wise gradient, $\varepsilon$ is the
symmetric gradient, $\operatorname{div}$ denotes the (row-wise) divergence the
superposed dot denotes the time derivative(s), and $I$ denotes the identity
matrix. On the right hand side, $f$ represents body forces and $g_{j}$ sources
(or sinks) in network $j$ for $j=1,\dots,J$.

The MPT equations represents a reduced version of (1.1) that result from
ignoring the elastic contribution of the solid matrix. These equations then
read as follows: for a given number of networks $J\in\mathbb{N}$, find the
network pressures $p_{j}$ for $j=1,\dots,J$ such that

(1.2)
$-\operatorname{div}K_{j}\operatorname{\nabla}p_{j}+\operatorname{\textstyle\sum}_{i=1}^{J}\xi_{j\leftarrow
i}(p_{j}-p_{i})=g_{j},$

where for $i,j=1,\dots,J$, $p_{j}=p_{j}(x)$ for
$x\in\Omega\subset{\mathbb{R}}^{d}$ ($d=1,2,3$), the parameters $K_{j}$ and
$\xi_{j\leftarrow i}$ remain the hydraulic conductivity and exchange
coefficients, respectively, and $g_{j}$ again represents other sources (or
sinks) in each network.

Parameter | Unit | Value | Reference
---|---|---|---
Hydraulic conductivities ($K_{j}$) | mm2 (kPa s)-1 |  |
Brain gray matter |  | $2.0\times 10^{-3}$ | [27]
Brain white matter |  | $2.0\times 10^{-2}$ | [27]
Cardiac arteries |  | $1.0$ | [19]
Cardiac capillaries |  | $2.0$ | [19]
Cardiac veins |  | $10.0$ | [19]
Brain vasculature |  | $3.75\times 10^{1}$ | [30]
Brain fluid exchange |  | $1.57\times 10^{-2}$ | [30]
Exchange coefficients ($\xi_{j\rightarrow i}$) | (kPa s)-1 |  |
Brain capillary-vasculature |  | $1.5\times 10^{-16}$ | [30]
Brain capillary-tissue fluid |  | $2.0\times 10^{-16}$ | [30]
Brain tissue fluid-veins |  | $2.0\times 10^{-10}$ | [30]
Cardiac capillary-arteries |  | $2.0\times 10^{-2}$ | [19]
Cardiac capillary-veins |  | $5.0\times 10^{-2}$ | [19]
Table 1: Sample parameter values for hydraulic conductivities and exchange
coefficients with reference to (1.1) and/or (1.2).

The relative size of the conductivities $K_{j}$ and the exchange coefficients
$\xi_{j\leftarrow i}$ may vary tremendously in applications. Large parameter
variation is certainly present in applied problems of a physiological nature;
a selection of representative parameter values, from research literature, is
given in Table 1.

Here, we see that the hydraulic conductivities span four orders of magnitude
while the exchange coefficients span fourteen orders of magnitude. Hence,
there is a need for preconditioners that are robust with respect to variations
in parameters. Physiological applications, in particular, can benefit from
preconditioners which are robust with respect to $K_{j}$, $\xi_{j\leftarrow
i}$ and $\lambda$ as in (1.1) and (1.2).

With this in mind, parameter-robust numerical approximations and solution
algorithms for (1.1) is currently an active research topic. In the nearly
incompressible case $\lambda\gg 1$, the standard two-field variational
formulation of (1.1) is not robust. To address this challenge, we introduced
and analyzed a mixed finite element method for the MPET equations based on a
total pressure formulation in [17]. We note that the total pressure in case of
one network was presented in [22, 16]. Hong et al. [9] shortly thereafter
extended the three-field formulation in [8] to MPET equations taking the
displacement, the network fluid fluxes and the network pressures as unknowns
and targeting a range of parameter regimes. As an alternative to these fully
coupled approaches a form of splitting schemes has been analyzed by Lee [15].
Regarding the iterative solution and preconditioning of the fully coupled
formulations, a robust preconditioner for Biot’s equations (the case for
$J=1$) was presented by Lee et al. [16]. Hong et al. [9] presented both
theoretical results and numerical examples regarding parameter-robust
preconditioners for the MPET equations with their extended three-field-type
formulation. Hong et al. further developed parameter-robust solver algorithms,
an iterative solver algorithm using the iterative coupling approach (cf. [20])
in [11], and an Uzawa-type algorithm in [10].

In this paper, we present a new approach for preconditioning linear systems of
equations resulting from a conforming finite element discretization of the
total pressure variational formulation of the MPET equations. The key idea, as
introduced for the MPT equations in [23], to design a parameter-dependent
transformation of the pressure variables $p=(p_{1},\dots,p_{J})$ into a set of
transformed variables $\tilde{p}$. The transformation should be such that the
originally coupled exchange operator decouples while the originally decoupled
diffusion operator stays decoupled (diagonal). The design of such a transform
hinges on the concept of _diagonalization by congruence_ and associated matrix
theory. After transformation, we then define a natural block diagonal
preconditioner for the transformed system of equations. This strategy yields a
parameter–robust preconditioner, which we both prove theoretically and
demonstrate numerically.

This manuscript is organized as follows. We introduce notation and review
relevant preconditioning and matrix theory in Section 2. We briefly consider
the reduced case of the MPT equations in Section 3 before turning to the
analysis of the preconditioner for the MPET equations in Section 4. Finally,
we present some conclusions and outlook in Section 5.

##  2 Preliminaries

In Section 2.2 we briefly review preconditioning of parameter-dependent
systems and state a known result regarding simultaneous diagonalization by
congruence. Notation for the remainder of the manuscript is discussed in
Section 2.1.

###  2.1 Notation

In the subsequent manuscript, we use the following notation. Let $\Omega$ be
an open, bounded domain in $\mathbb{R}^{d}$, $d=2,3$, with Lipschitz
polyhedral boundary $\partial\Omega$. We denote by $L^{2}(\Omega)$ the space
of square integrable functions on $\Omega$ with inner product
$\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. We denote by $H^{m}(\Omega)$
the standard Sobolev space with norm $\lVert\cdot\rVert_{H^{m}}$ and semi-norm
$|\cdot|_{H^{m}}$ for $m\geq 1$ and $H^{m}(\Omega;{\mathbb{R}}^{d})$ the
corresponding $d$-vector fields. We use $H^{m}_{0}(\Omega)$ to denote the
subspace of $H^{m}(\Omega)$ with vanishing trace on the boundary of $\Omega$.
Let $\Gamma$ be a subset of $\partial\Omega$ such that
$\partial\Omega\setminus\Gamma$ has a positive $(d-1)$-dimensional Lebesgue
measure. $H_{\Gamma}^{m}(\Omega)$ is the subspace of $H^{m}(\Omega)$ such that
the elements in $H_{\Gamma}^{m}(\Omega)$ have vanishing trace on $\Gamma$.
$H_{\Gamma}^{m}(\Omega;{\mathbb{R}}^{d})$ is the subspace of
$H^{m}(\Omega;{\mathbb{R}}^{d})$ such that every $v_{j}$ in
$(v_{1},\ldots,v_{d})\in H^{m}(\Omega;{\mathbb{R}}^{d})$ is an element in
$H_{\Gamma}^{m}(\Omega)$. Throughout this paper we set $\Gamma$ a fixed subset
of $\partial\Omega$ satisfying the aforementioned assumption.

We introduce the parameter-dependent $L^{2}$-inner product and norm:

$\|p\|_{\beta}^{2}=\langle p,p\rangle_{\beta}=\langle\beta p,p\rangle$

for $\beta\in L^{\infty}(\Omega)$, $\beta(x)>0$, and $p\in L^{2}(\Omega)$ (and
similarly for vector or tensor fields). The notation $\mathbb{I}$ will denote
an identity $d\times d$ matrix while $I_{V}$ will denote the identity operator
on a Hilbert space $V$.

To be self-contained we recall the Kronecker product of matrices. If $A$ in
$\mathbb{R}^{m\times n}$ and $B\in\mathbb{R}^{r\times s}$ are two real-valued
matrices then $A\otimes B$ is the $mr\times ns$ matrix defined by multiplying
each entry of $A$ by the matrix $B$. That is,

(2.5) $\displaystyle A\otimes
B=\left[\begin{array}[]{cccc}a_{11}B&a_{12}B&\dots&a_{1n}B\\\
a_{21}B&a_{22}B&\dots&a_{2n}B\\\ \vdots&\ddots&\ddots&\vdots\\\
a_{m1}B&a_{m2}B&\dots&a_{mn}B\end{array}\right].$

We can consider its natural extension for a matrix $A$ and a linear operator
$B$. More specifically, if $W$ is a Hilbert space, $Q$ is the $n-$fold product
$Q=W\times W\times\dots\times W$, $A$ is an $n\times n$ matrix, and $B$ is a
linear operator on $W$, then $A\otimes B$ is the linear operator on $Q$
defined by (2.5).

Finally, we introduce a notation for uniform proportionality, used throughout
the manuscript, as

$X\lesssim Y.$

That is, $X\lesssim Y$ implies the existence of some real constant $c_{0}>0$
such that $X\leq c_{0}Y$; any relationship between $c_{0}$ and pertinent
mathematical objects, such as the total number of porous media networks
considered, will be specified.

###  2.2 Preconditioning of parameter-dependent systems

In this paper, we consider the preconditioning of discretizations of the
systems (1.1) and (1.2) under large parameter variations. Therefore, we begin
by summarizing core aspects of the theory of parameter–robust preconditioning
as presented in [18]. We will apply this theory for formulations of the MPT
equations (1.2) and MPET equations (1.1) in the subsequent sections.

Let $X$ be a separable, real Hilbert space with inner product
$\langle\cdot,\cdot\rangle_{X}$, norm $\|\cdot\|_{X}$ and dual space $X^{*}$.
Let $\mathcal{A}:X\rightarrow X$ be an invertible, symmetric isomorphism on
$X$ such that $\mathcal{A}\in\mathcal{L}(X,X^{*})$ where
$\mathcal{L}(X,X^{*})$ is the set of bounded linear operators mapping $X$ to
its dual. Given $f\in X^{*}$ consider the problem of finding $x\in X$ such
that

(2.6) $\mathcal{A}x=f.$

The preconditioned problem reads as follows

(2.7) $\mathcal{BA}x=\mathcal{B}f,$

where $\mathcal{B}\in\mathcal{L}(X^{*},X)$ is a symmetric isomorphism defining
the preconditioner. The convergence rate of a Krylov space method for this
problem can be bounded in terms of the condition number $\kappa(\mathcal{BA})$
where

$\kappa(\mathcal{BA})=\lVert\mathcal{BA}\rVert_{\mathcal{L}(X,X)}\lVert(\mathcal{BA})^{-1}\rVert_{\mathcal{L}(X,X)}.$

Here, the operator norm $\lVert\mathcal{A}\rVert_{\mathcal{L}(X,X^{*})}$ is
defined by

(2.8) $\lVert\mathcal{A}\rVert_{\mathcal{L}(X,X^{*})}=\sup_{x\in
X}\frac{\lVert\mathcal{A}x\rVert_{X^{*}}}{\lVert x\rVert_{X}}.$

Now, for a parameter $\varepsilon$ (or more generally a set of parameters
$\varepsilon$) consider the parameter-dependent operator
$\mathcal{A}_{\varepsilon}$ and its preconditioner
$\mathcal{B}_{\varepsilon}$. Assume that we can choose appropriate spaces
$X_{\varepsilon}$ and $X^{*}_{\varepsilon}$ such that the norms

$\|\mathcal{A}_{\varepsilon}\|_{\mathcal{L}(X_{\varepsilon},X^{*}_{\varepsilon})}\text{
and
}\|\mathcal{A}^{-1}_{\varepsilon}\|_{\mathcal{L}(X^{*}_{\varepsilon},X_{\varepsilon})}$

are bounded independently of $\varepsilon$. Similarly, we assume that we can
find a preconditioner $\mathcal{B}_{\varepsilon}$ such that the norms
$\lVert\mathcal{B}_{\varepsilon}\rVert_{\mathcal{L}(X_{\varepsilon},X^{*}_{\varepsilon})}$
and
$\lVert\mathcal{B}^{-1}_{\varepsilon}\rVert_{\mathcal{L}(X^{*}_{\varepsilon},X_{\varepsilon})}$
are bounded independently of $\varepsilon$. Given these assumptions, the
condition number $\kappa(\mathcal{B_{\varepsilon}A_{\varepsilon}})$ will be
bounded independently of $\varepsilon$. We will refer to such a preconditioner
as robust in (or with respect to) $\varepsilon$. We conclude this section with
a change of variables result, recalled from basic matrix analysis [12], that
will prove effective in the sections that follow.

######  Lemma 2.1.

Let $W$ be a real Hilbert space and $Q=W\times W\times\dots\times W$ be the
$n$-fold direct product of $W$ for a fixed $n\in\mathbb{N}$. Let
$A:W\rightarrow W^{\ast}$, and $B:W\rightarrow W^{\ast}$ be linear operators.
Suppose that $K,E\in\mathbb{R}^{n\times n}$ are symmetric matrices, and $K$ is
invertible. Define the operators $S:Q\rightarrow Q^{\ast}$ and $T:Q\rightarrow
Q^{\ast}$ by

$S=K\otimes A,\quad\text{and}\quad T=E\otimes B,$

where $\otimes$ is the Kronecker product. Consider the variational problem:
given $f\in Q^{*}$ find $p=(p_{1},p_{2},\dots,p_{n})^{T}\in Q$ such that

(2.9) $\langle Sp,q\rangle+\langle Tp,q\rangle=\langle
f,q\rangle,\quad\forall\,q\in Q$

where $\langle\cdot,\cdot\rangle$ is the duality pairing of $Q^{*}$ and $Q$.
Then there exists an invertible matrix $P\in\mathbb{R}^{n\times n}$ such that
the above variational problem is equivalent to: find $\tilde{p}\in Q$ such
that

(2.10) $\langle D_{S}\tilde{p},q\rangle+\langle
D_{T}\tilde{p},q\rangle=\langle F,q\rangle,\quad\forall\,q\in Q,$

where $F=(P^{T}\otimes I_{W^{*}})f$ for $I_{W^{*}}$ the identity operator on
$W^{*}$, and $D_{S}=\left(P^{T}KP\right)\otimes A$ and
$D_{T}=\left(P^{T}EP\right)\otimes B$ are block diagonal linear operators from
$Q$ to $Q^{*}$.

######  Proof 2.2.

Apply [12, Theorem 4.5.17a-b p. 287] the hypotheses on the matrices $K$ and
$E$ and properties of the tensor product; see Appendix A for detail.

##  3 Preconditioning the MPT equations via diagonalization

We begin by summarizing our novel approach to variational formulations and
associated preconditioning of the MPT equations. This approach was introduced
in [23]. The core idea is to reformulate the MPT (and MPET) equations using a
change of pressure variables $p$. In particular, we aim to find a
transformation of the variables $p\mapsto\tilde{p}$ such that the transformed
system of pressure equations decouple. Here, we will consider a Hilbert space
$W$ and the $J$-fold product $Q=W\times W\times\dots\times W$. Each pressure
$p_{j}$, for $j=1,2,\dots,J$ satisfies $p_{j}\in W$ and we will write
$p=(p_{1},p_{2},\dots,p_{J})\in Q$. In the sections that follow, we briefly
illustrate the core idea, formulation of the MPT equations and resulting
preconditioner, and refer to [23] for more details. This approach is then
extended to the MPET equations in Section 4.

###  3.1 The MPT equations in operator form

We consider the MPT equations as defined by (1.2). We further impose
homogeneous Dirichlet boundary conditions for all pressures: $p_{j}=0$ on
$\partial\Omega$ for $1\leq j\leq J$. Define
$\xi_{j}=\operatorname{\textstyle\sum}_{i=1}^{J}\xi_{j\leftarrow i}$ for each
$1\leq j\leq J$. Let us define two $J\times J$ matrices:

(3.1) $K=\begin{pmatrix}K_{1}&0&\cdots&0\\\ 0&K_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&K_{J}\\\ \end{pmatrix},\quad
E=\begin{pmatrix}\xi_{1}&-\xi_{1\leftarrow 2}&\cdots&-\xi_{1\leftarrow J}\\\
-\xi_{1\leftarrow 2}&\xi_{2}&\cdots&-\xi_{2\leftarrow J}\\\
\vdots&\vdots&\ddots&\vdots\\\ -\xi_{1\leftarrow J}&-\xi_{2\leftarrow
J}&\cdots&\xi_{J}\\\ \end{pmatrix}.$

The system (1.2) can be expressed in operator form as: given $g\in Q$ find
$p\in Q$ satisfying

(3.2) $\mathcal{A}_{\rm MPT}p=g\quad\text{ where }\quad\mathcal{A}_{\rm
MPT}=-K\otimes\Delta+E\otimes I_{W}.$

In the above, $-K\otimes\Delta:Q\rightarrow Q^{*}$ is the block diagonal
operator such that its $j$-th block is given by the bilinear form $\langle
K_{j}\nabla p_{j},\nabla q_{j}\rangle$ for $p_{j},q_{j}\in Q_{j}=W$, and
$E\otimes I_{W}:Q\rightarrow Q^{*}$ is the block operator such that its
$(i,j)$-block $E_{ij}$ is defined by the bilinear forms

$-\langle\xi_{i\leftarrow j}p_{i},q_{j}\rangle\text{ if
}i\not=j,\qquad\langle\xi_{j}p_{j},p_{j}\rangle\text{ if }i=j.$

We note that $K$ is real, positive definite and diagonal (and thus
invertible), and that $E$ is real, symmetric and (weakly row) diagonally
dominant by definition. In particular, $E$ is symmetric positive semi-definite
because of the identity

(3.3) $wEw^{T}=\sum_{1\leq i,j\leq J}\xi_{i\leftarrow j}(w_{i}-w_{j})^{2},$

for $w=(w_{1},w_{2},\ldots,w_{J})$ with the convention $\xi_{i\leftarrow
i}=0$. A naive block diagonal preconditioner $\mathcal{B}_{\rm MPT}$ can be
constructed by taking the inverse of the diagonal blocks of $\mathcal{A}_{\rm
MPT}$. However, as we demonstrated in [23], the resulting preconditioner is
not robust with respect to variations in the conductivity and exchange
parameters. In fact, the condition numbers increased linearly with the ratio
between the exchange and conductivity coefficients.

###  3.2 Change of variables using diagonalization by congruence

In this section we discuss a new formulation for the MPT equations, which in
turn easily offers a parameter-robust preconditioner. Let
$P\in{\mathbb{R}}^{J\times J}$ be an invertible linear transformation defining
a change of variables and let $\tilde{p}$ and $\tilde{q}$ be the new set of
variables such that

(3.4) $p=\left(P\otimes I_{W}\right)\tilde{p},\quad q=\left(P\otimes
I_{W}\right)\tilde{q},$

with $q=(q_{1},q_{2},\dots,q_{J})$ and similarly for $\tilde{q},\tilde{p}$.
Since $K$ and $E$ are symmetric, we apply Lemma 2.1, with $A=\Delta$ and
$B=I_{W}$, to obtain a matrix, $P$, simultaneously diagonalizing $K$ and $E$
by congruence; that is, the equivalent operators
$\left(P^{T}KP\right)\otimes\Delta$ and $(P^{T}EP)\otimes I_{W}$ are block
diagonal. The resulting formulation (c.f. (2.10)) of the MPT equations reads
as follows: find the transformed pressures

$\tilde{p}=(\tilde{p}_{1},\dots,\tilde{p}_{J})$ such that, for a given $g\in
Q$, we have the equality

(3.5) $\mathcal{\tilde{A}_{\rm
MPT}}\tilde{p}=(-\tilde{K}\otimes\Delta+\tilde{E}\otimes
I_{W})\tilde{p}=\left(P^{T}\otimes I_{W}\right)g,$

where $\tilde{K}=P^{T}KP$ and $\tilde{E}=P^{T}EP$ are diagonal with

(3.6)
$\tilde{K}=\operatorname{diag}(\tilde{K}_{1},\dots,\tilde{K}_{J}),\quad\tilde{E}=\operatorname{diag}(\tilde{\xi}_{1},\dots,\tilde{\xi}_{J}).$

###  3.3 Preconditioning the transformed MPT system

As in [23], we can immediately identify the parameter-dependent norm

$\lVert\tilde{p}\rVert^{2}_{\tilde{\mathcal{B}}_{\rm
MPT}}=\sum_{j=1}^{J}\langle\tilde{K}_{j}\operatorname{\nabla}\tilde{p}_{j},\operatorname{\nabla}\tilde{p}_{j}\rangle+\langle\tilde{\xi}_{j}\tilde{p}_{j},\tilde{p}_{j}\rangle$

and the following preconditioner associated to the above norm for (3.5):

(3.7) $\mathcal{\tilde{B}}_{\rm
MPT}=\begin{pmatrix}(-\tilde{K}_{1}\Delta+\tilde{\xi}_{1}I)^{-1}&0&\cdots&0\\\
0&(-\tilde{K}_{2}\Delta+\tilde{\xi}_{2}I)^{-1}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&(-\tilde{K}_{J}\Delta+\tilde{\xi}_{J}I)^{-1}\\\ \end{pmatrix}.$

Clearly, $\mathcal{\tilde{A}_{\rm MPT}}$ and $\mathcal{\tilde{B}}^{-1}_{\rm
MPT}$ are trivially spectrally equivalent. We refer to [23] for numerical
experiments comparing the standard and transformed formulation and
preconditioners.

###  3.4 Finding the transformation matrix

There are two cases that we will consider; the first case is when the matrix
$C=K^{-1}E$ has $J$ distinct eigenvalues, while the second case will be for
the case where at least one of the eigenvalues is repeated. In the case of
distinct eigenvalues, the number of distinct eigenvalues of $C=K^{-1}E$ will
depend on the material parameter values $K_{j}$ and $\xi_{j\rightarrow i}$ for
$1\leq i,j\leq J$. In the common case where $C$ has $J$ distinct eigenvalues,
the transformation matrix is easily defined as follows. Let
$\lambda_{1},\dots,\lambda_{J}$ be the real eigenvalues of $C$, and let
$v_{1},\dots,v_{J}$ be the corresponding normalized eigenvectors. Then,

(3.8) $P=[v_{1},\dots,v_{J}],$

will diagonalize $K$ and $E$ by congruence. In [23], we presented numerical
examples for the case of $J$ distinct eigenvalues (with $J=2$).

The congruence matrix for the case of repeated eigenvalues is also easily
constructed. For these cases, the transform $P$ can be constructed by repeated
application of block-wise eigenvector matrices, see [12] for the general
procedure. In Example 3.1 below, we present an example on how to obtain the
transformation matrix $P$ in the case where one of the eigenvalues has
algebraic multiplicity $2$ with $J=3$.

######  Example 3.1.

In this example we show how to obtain the transformation matrix $P$ for a
three–network case when one of the eigenvalues of $K^{-1}E$ has algebraic
multiplicity $2$. In this example, due to the presence of the repeated
eigenvector, the construction of $P$ does not follow directly from the use of
normalized eigenvectors and, thus, $P$ is not normalized a priori. We will,
however, normalize $P$ following construction to maintain consistency with the
previous case; in practice, either the normalized or non-normalized version of
$P$ may be used.

(3.9) $K=\begin{pmatrix}1.0&0&0\\\ 0&0.0001&0\\\ 0&0&0.01\\\
\end{pmatrix},\quad E=\begin{pmatrix}1.01&-0.01&-1.0\\\
-0.01&0.0101&-0.0001\\\ -1.0&-0.0001&1.0001\end{pmatrix}.$

By definition

(3.10) $C=K^{-1}E=\begin{pmatrix}1.01&-0.01&-1.0\\\ -100&101&-1.0\\\
-100&-0.01&100.01\end{pmatrix}.$

The eigenvalues $\lambda_{1},\lambda_{2},\lambda_{3}$ and eigenvectors
$[v_{1},v_{2},v_{3}]=P_{1}$ of $C$ are then:

(3.11) $\begin{split}\lambda_{1}&=0,\quad\lambda_{2}=\lambda_{3}=101.01;\\\
P_{1}&=\begin{pmatrix}-0.5773&-0.0071&-0.0091\\\ -0.5773&0.7070&-0.4031\\\
-0.5773&0.7070&0.9150\end{pmatrix}.\end{split}$

In this specific case the eigenvalues $\lambda_{2},\lambda_{3}$ have algebraic
multiplicity $2$ and geometrical multiplicity $1$. If we try to diagonalize
$K$ and $E$ by congruence via $P_{1}$, we obtain

(3.12) $\begin{split}P_{1}^{T}KP_{1}&=\begin{pmatrix}3.3670\times
10^{-1}&0&0\\\ 0&5.1007\times 10^{-3}&6.5069\times 10^{-3}\\\ 0&6.5069\times
10^{-3}&8.4729\times 10^{-3}\end{pmatrix},\\\
P_{1}^{T}EP_{1}&=101.01\begin{pmatrix}0&0&0\\\ 0&5.1007\times
10^{-3}&6.5069\times 10^{-3}\\\ 0&6.5069\times 10^{-3}&8.4729\times
10^{-3}\end{pmatrix}.\end{split}$

In this case, the resulting matrices are block diagonal. The lower right
blocks are multiples of each other. We can diagonalize the lower right blocks
by computing the eigen-decomposition of either of these. The lower right block
of $P_{1}^{T}KP_{1}$ is

(3.13) $\begin{pmatrix}5.1007\times 10^{-3}&6.5069\times 10^{-3}\\\
6.5069\times 10^{-3}&8.4729\times 10^{-3}\end{pmatrix}$

and its eigenpairs are

(3.14) $\begin{split}\lambda_{1}&=6.4967\times
10^{-5},\lambda_{2}=1.3508\times 10^{-2};\\\
P_{2}&=\begin{pmatrix}-0.79083&-0.6120\\\
0.6120&-0.7908\end{pmatrix}.\end{split}$

The final transformation matrix $P$ that diagonalizes $K$ and $E$ by
congruence is then:

(3.15) $P=P_{1}\begin{pmatrix}1&\begin{matrix}0&0\end{matrix}\\\
\begin{matrix}0\\\
0\end{matrix}&P_{2}\end{pmatrix}=\begin{pmatrix}-5.7735\times
10^{-1}&7.1935\times 10^{-5}&1.1575\times 10^{-2}\\\ -5.7735\times
10^{-1}&-8.0594\times 10^{-1}&-1.1391\times 10^{-1}\\\ -5.7735\times
10^{-1}&8.6590\times 10^{-4}&-1.1564\end{pmatrix}.$

Note that despite the columns of $P_{1}$ and $P_{2}$ are normalized with norm
$1,$ the resulting matrix $P^{\prime}$s columns are not normalized. After the
normalization, the matrix $P$ looks as follows:

(3.16) $P=\begin{pmatrix}-5.7735\times 10^{-1}&8.9255\times
10^{-5}&9.9611\times 10^{-2}\\\ -5.7735\times 10^{-1}&-9.9999\times
10^{-1}&-9.8026\times 10^{-2}\\\ -5.7735\times 10^{-1}&1.0743\times
10^{-4}&-9.9513\times 10^{-1}\end{pmatrix}.$

and the diagonalized matrices are as follows

(3.17) $\begin{split}\tilde{K}=P^{T}KP&=\begin{pmatrix}3.3670\times
10^{-1}&0&0\\\ 0&1.0001\times 10^{-4}&0\\\ 0&0&1.0003\times
10^{-2}\end{pmatrix},\\\ \tilde{E}=P^{T}EP&=\begin{pmatrix}0&0&0\\\
0&1.0102\times 10^{-2}&0\\\ 0&0&1.0104\end{pmatrix}.\end{split}$

##  4 Preconditioning the MPET equations via diagonalization

In this section, we present a change of variables for the total pressure
formulation of the time-discrete MPET equations and propose and analyze a
preconditioning strategy for the resulting variational formulation. The change
of MPET variables is guided by the change of MPT variables presented in the
previous section. As in Sec. 3, $Q$ will be defined as the $J$-fold product of
a Hilbert space $W$.

###  4.1 Total pressure formulation of the MPET equations and its time-
discrete form

The total pressure formulation of Biot’s equations [16] and more generally the
MPET equations [17] is a robust mixed variational formulation targeting the
nearly incompressible case and incompressible limit ($\lambda\gg 1$). The
total pressure, which we will see satisfies $p_{0}\in L^{2}(\Omega)$, is
defined by

(4.1) $p_{0}=\lambda\operatorname{div}u-\alpha\cdot p,$

where 111Note that we start counting at $1$ in the definition of $p$ here and
throughout, in contrast to e.g. in [23].
$\alpha=(\alpha_{1},\dots,\alpha_{J})\in\mathbb{R}^{J}$,
$p=(p_{1},\dots,p_{J})\in Q$ and $\alpha\cdot
p=\operatorname{\textstyle\sum}_{i=1}^{J}\alpha_{i}p_{i}\in W$. The total
pressure formulation of (1.1) then reads as follows: for $t\in(0,T]$, find the
displacement vector field $u$ and the pressure scalar fields $p_{0}$ and
$p_{j}$ for $j=1,\dots,J$ such that

(4.2a)
$\displaystyle-\operatorname{div}\left(2\mu\varepsilon(u)+p_{0}\mathbb{I}\right)=f,$
(4.2b)
$\displaystyle\operatorname{div}u-\lambda^{-1}p_{0}-\lambda^{-1}\alpha\cdot
p=0,$ (4.2c)
$\displaystyle\lambda^{-1}\dot{p}_{0}+s_{j}\dot{p}_{j}-\operatorname{div}(K_{j}\nabla
p_{j})+\alpha_{j}\lambda^{-1}\alpha\cdot\dot{p}+\operatorname{\textstyle\sum}_{i=1}^{J}\xi_{j\leftarrow
i}(p_{j}-p_{i})=g_{j},$

for $j=1,\dots,J$.

We consider an implicit Euler discretization in time of the total pressure
formulation of the time-dependent MPET equations (4.2) and examine the
resulting stationary problem at each time step. The resulting time-discrete
version of (4.2) with time step $\tau>0$ reads as follows: find the
displacement $u$ and the pressures $p_{j}$ for $0\leq j\leq J$ such that

(4.3a)
$\displaystyle-\operatorname{div}\left(2\mu\varepsilon(u)+p_{0}\mathbb{I}\right)=f,$
(4.3b)
$\displaystyle\operatorname{div}u-\lambda^{-1}p_{0}-\lambda^{-1}\alpha\cdot
p=0,$ (4.3c) $\displaystyle-
s_{j}p_{j}-\alpha_{j}\lambda^{-1}p_{0}-\alpha_{j}\lambda^{-1}\alpha\cdot
p+\tau\operatorname{div}(K_{j}\nabla
p_{j})-\tau\operatorname{\textstyle\sum}_{i=1}^{J}\xi_{j\leftarrow
i}(p_{j}-p_{i})=g_{j},$

for $1\leq j\leq J$ where the new right hand sides $g_{j}$ for $j=1,\dots,J$
have been negated and contain also terms from the previous time-step. Again,
we impose homogeneous Dirichlet boundary conditions for all network pressures:
$p_{j}=0$ on $\partial\Omega$ for $1\leq j\leq J$.

Let $V=H^{1}_{\Gamma}(\Omega;{\mathbb{R}}^{d})$, $Q_{0}=L^{2}(\Omega)$ and
$Q_{j}=W=H^{1}_{0}(\Omega)$ for $1\leq j\leq J$ and
$\Omega\subset{\mathbb{R}}^{d}$. Let $Q=Q_{1}\times\dots\times Q_{J}$. As in
Section 3, we write $p=(p_{1},\dots,p_{J})$, $q=(q_{1},\ldots,q_{J})$, and
$g=(g_{1},\ldots,g_{J})$. Multiplying by test functions, and integrating
second-order derivatives by parts, we obtain the following variational
formulation of (4.3): find $u\in V$ and $p_{i}\in Q_{i}$ for $i=0,\dots,J$
such that

(4.4a) $\displaystyle a(u,v)+b(v,p_{0})$ $\displaystyle=\langle f,v\rangle$
$\displaystyle\quad\forall\,v\in V,$ (4.4b) $\displaystyle
b(u,q_{0})-c_{1}(p_{0},q_{0})-c_{2}(q_{0},p)$ $\displaystyle=0$
$\displaystyle\quad\forall\,q_{0}\in Q_{0},$ (4.4c) $\displaystyle-
c_{2}(p_{0},q)-c_{3}(p,q)$ $\displaystyle=\langle g,q\rangle$
$\displaystyle\quad\forall\,q\in Q.$

The bilinear forms $a:V\times V\rightarrow{\mathbb{R}}$ and $b:V\times
Q_{0}\rightarrow{\mathbb{R}}$ are defined as:

(4.5) $\displaystyle a(u,v)$ $\displaystyle=\langle
2\mu\varepsilon(u),\varepsilon(v)\rangle,$ $\displaystyle b(v,q_{0})$
$\displaystyle=\langle\operatorname{div}v,q_{0}\rangle,$

while $c_{1}:Q_{0}\times Q_{0}\rightarrow{\mathbb{R}}$, $c_{2}:Q_{0}\times
Q\rightarrow{\mathbb{R}}$, and $c_{3}:Q\times Q\rightarrow{\mathbb{R}}$ are
defined as:

(4.6) $\displaystyle c_{1}(p_{0},q_{0})$
$\displaystyle=\langle\lambda^{-1}p_{0},q_{0}\rangle,$ (4.7) $\displaystyle
c_{2}(p_{0},q)$ $\displaystyle=\langle\lambda^{-1}\alpha\cdot q,p_{0}\rangle,$
(4.8) $\displaystyle c_{3}(p,q)$ $\displaystyle=\tau\sum_{j=1}^{J}\langle
K_{j}\operatorname{\nabla}p_{j},\operatorname{\nabla}q_{j}\rangle+\sum_{j=1}^{J}\langle
s_{j}p_{j},q_{j}\rangle,$
$\displaystyle\quad+\tau\sum_{j=1}^{J}\sum_{i=1}^{J}\langle\xi_{j\leftarrow
i}(p_{j}-p_{i}),q_{j}\rangle+\langle\lambda^{-1}\alpha\cdot p,\alpha\cdot
q\rangle.$

For future reference we define $c:(Q_{0}\times Q)\times(Q_{0}\times
Q)\rightarrow{\mathbb{R}}$ via

(4.9) $\displaystyle c((p_{0},p),(q_{0},q))$
$\displaystyle=c_{1}(p_{0},q_{0})+c_{2}(p_{0},q)+c_{2}(q_{0},p)+c_{3}(p,q).$

###  4.2 MPET as a parameter-dependent saddle point system

The purpose of this subsection is to explain difficulties to construct
parameter-robust block preconditioners for the system (4.4). The difficulties
lead us to develop the change of variables idea in the MPT problem in
consideration of its extension to MPET problems.

Note that the system (4.3) or equivalently (4.4) can be viewed as a saddle
point problem with a stabilization term (given by the bilinear form $c$).
Thus, the equations fit well into Brezzi saddle point theory [4]. However,
various material parameters in different ranges are involved in the system, so
constructing parameter-robust preconditioners for this system is not a
straightforward application of the Brezzi theory. Let us recall the parameter
ranges we are concerned with in this paper. The existing literature cover the
parameters

$0\leq s_{j}\lesssim 1,\qquad 0<K_{j}\ll 1,\qquad
1\lesssim\mu\lesssim\lambda<+\infty,\qquad(1\leq j\leq J).$

In addition to these we are also interested in developing preconditioners
which are robust for the ratios of the exchange coefficients
$\xi_{i\rightarrow j}$’s.

We first consider construction of preconditioners utilizing the saddle point
problem structure. To reveal the saddle point problem structure of (4.4) let
us look at the operator form of (4.4), which is

(4.10) $\mathcal{A}_{\rm MPET}\begin{pmatrix}u\\\ p_{0}\\\
p\end{pmatrix}=\begin{pmatrix}-2\operatorname{div}(\mu\varepsilon)&-\operatorname{\nabla}&\textbf{0}\\\
\operatorname{div}&-C_{1}&-{C}_{2}^{*}\\\ \textbf{0}&-{C}_{2}&-{C}_{3}\\\
\end{pmatrix}\begin{pmatrix}u\\\ p_{0}\\\ p\end{pmatrix}=\begin{pmatrix}f\\\
0\\\ g\\\ \end{pmatrix}$

where $C_{1}:Q_{0}\rightarrow Q_{0}^{*}$, $C_{2}:Q_{0}\rightarrow Q^{*}$,
${C}_{3}:Q\rightarrow Q^{*}$ are the operators associated to the bilinear
forms $c_{1}$, $c_{2}$, $c_{3}$ in (4.6), (4.7), (4.8). Here $C_{2}^{*}$ is
the adjoint operator of $C_{2}$. We can rewrite $\mathcal{A}_{\rm MPET}$ of
(4.10) in the standard saddle point form

$\mathcal{A}_{\rm MPET}=\begin{pmatrix}A&B_{0}^{*}\\\ B_{0}&-C\end{pmatrix}$

by considering the product space grouping $V\times(Q_{0}\times Q)$ and
identifying

(4.11) $A=-2\operatorname{div}(\mu\varepsilon),\qquad
B_{0}=(\operatorname{div},\textbf{0})^{T},\qquad
C=\begin{pmatrix}C_{1}&C_{2}^{*}\\\ C_{2}&C_{3}\\\ \end{pmatrix}.$

One of natural approaches to construct block preconditioners for this system
is to use the block diagonal operator

$\begin{pmatrix}A^{-1}&0\\\ 0&(C+B_{0}A^{-1}B_{0}^{*})^{-1}\end{pmatrix}$

or its approximation. However, the operator $(C+B_{0}A^{-1}B_{0}^{*})^{-1}$ is
not easy to implement efficiently in practice. Moreover, the analysis for
spectral equivalence of this type of preconditioners is related to a non-
trivial generalized eigenvalue problem. More precisely, the spectral
equivalence is equivalent to uniform upper and lower bounds of the generalized
eigenvalues, so it requires a deep analysis of the non-trivial generalized
eigenvalue problem utilizing block matrix structures. In this paper we
consider a general MPET model with general $J$ and general (constant) exchange
coefficients, so the number of blocks in block matrices is not restricted.
This makes an analysis of the generalized eigenvalue problem even more
challenging, so we will not pursue this approach further in this paper.

Another natural choice of block preconditioners for this system is a direct
extension of the preconditioner in [16]. In other words, we use the block
diagonal operator of the form

(4.12) $\mathcal{B}_{\rm MPET}=\begin{pmatrix}(-\mu\Delta)^{-1}&0&0\\\
0&I^{-1}&0\\\ 0&0&D^{-1}\\\ \end{pmatrix},\;$

where $I:Q_{0}\rightarrow Q_{0}^{*}$ is the operator defined by the bilinear
form $\langle p_{0},q_{0}\rangle$ for $p_{0},q_{0}\in Q_{0}$, $D:Q\rightarrow
Q^{*}$ is the block diagonal operator such that its $j$-th diagonal block
($1\leq j\leq J$) is defined by selecting the $j^{\text{th}}$ diagonal entry
of the operator $C_{3}$ associated to the bilinear form (4.8); that is

$\tau\langle
K_{j}\operatorname{\nabla}p_{j},\operatorname{\nabla}q_{j}\rangle+\langle
s_{j}p_{j},q_{j}\rangle+\tau\langle\xi_{j}p_{j},q_{j}\rangle+\langle\lambda^{-1}\alpha_{j}p_{j},\alpha_{j}q_{j}\rangle,\qquad
p,q\in Q.$

However, this preconditioner is not robust with respect to the material
parameters, particularly for the hydraulic conductivity and the exchange
coefficients. We illustrate numerical experiment results in Example 4.1.

######  Example 4.1.

Let $\Omega=[0,1]^{2}\subset{\mathbb{R}}^{2}$, and consider a structured
triangulation $\mathcal{T}_{h}$ of $\Omega$ constructed by dividing $\Omega$
into $N\times N$ squares and then subdividing each square by a fixed diagonal.
Let $J=2$. Consider a finite element discretization of (4.4) using the lowest
order Taylor–Hood-type elements i.e. continuous piecewise quadratics for each
displacement component, and continuous piecewise linear for all pressures
[17]. Let $\tau=1.0$, $\mu=1.0$, $s_{j}=1.0$, $\alpha_{j}=0.5$ and
$K_{1}=1.0$, and consider ranges of values for $\lambda,\xi_{1\leftarrow 2}$
and $K_{2}$. We consider the case for $s_{1}=s_{2}=1.0$ and $s_{1}=s_{2}=0.0.$
Starting from an initial random guess, we consider a MinRes solver of the
resulting linear system of equations with an algebraic multigrid (Hypre AMG)
preconditioner of the form (4.12). The convergence criterion used was

$(\mathcal{B}r_{k},r_{k})/(\mathcal{B}r_{0},r_{0})\leq 10^{-6}$

where $r_{k}$ is the residual of the $k$-th iteration. The resulting number of
Krylov iterations are shown in Table 2 for $\xi_{1\rightarrow 2}=10^{6}$ and
ranges of $K_{2}$ and $\lambda$. We observe that the number of iterations is
moderate ($\approx 30$) for $K_{2}$ of comparable magnitude ($10^{6}$) to
$\xi_{1\leftarrow 2}$. The number of iterations increase with decreasing
$K_{2}$: up to $\approx 1000$ for $K_{2}=1$. For large $K_{2}$, the number of
iterations seems independent of the mesh resolution $N$. In contrast, for
smaller $K_{2}$ (relative to $\xi_{1\rightarrow 2}$), the number of iterations
also increase with the mesh resolution. We note that the iteration counts do
not vary substantially with $\lambda$.

$K_{2}$ | $\lambda$ | $N=16$ | $32$ | $64$ | $128$
---|---|---|---|---|---
$10^{0}$ | $10^{0}$ | $738$ | $1271$ | $1756$ | $1938$
$10^{2}$ | $1024$ | $1505$ | $1679$ | $1631$
$10^{4}$ | $1028$ | $1506$ | $1666$ | $1628$
$10^{6}$ | $1004$ | $1499$ | $1677$ | $1633$
$10^{2}$ | $10^{0}$ | $396$ | $424$ | $406$ | $353$
$10^{2}$ | $337$ | $368$ | $351$ | $333$
$10^{4}$ | $364$ | $352$ | $348$ | $332$
$10^{6}$ | $345$ | $357$ | $361$ | $328$
$10^{4}$ | $10^{0}$ | $65$ | $65$ | $62$ | $60$
$10^{2}$ | $64$ | $60$ | $56$ | $55$
$10^{4}$ | $62$ | $60$ | $57$ | $55$
$10^{6}$ | $63$ | $61$ | $58$ | $55$
$10^{6}$ | $10^{0}$ | $30$ | $30$ | $30$ | $28$
$10^{2}$ | $34$ | $31$ | $29$ | $29$
$10^{4}$ | $32$ | $31$ | $31$ | $29$
$10^{6}$ | $33$ | $31$ | $31$ | $29$
Table 2: Number of MinRes iterations (c.f. Example 4.1): (4.4) as discretized
with Taylor-Hood type elements and an algebraic multigrid preconditioner of
the form (4.12). Of note is the fact that the number of iterations grow for
$K_{2}$ decreasing relative to $\xi_{2\rightarrow 1}=10^{6}$, and for
increasing $N$.

###  4.3 Transformed MPET equations via change of variables

In this subsection, we present MPET equations which are transformed via change
of variables for construction of block preconditioners. As in the MPT problem
we will find an invertible linear map $P\in{\mathbb{R}}^{J\times J}$ that
provides a fortuitous co-diagonalization; we will then consider the change of
variables

$\displaystyle p=\left(P\otimes I_{W}\right)\tilde{p},$

which will lead to a (partial) diagonalization, in the spirit of Lemma 2.1,
for the transformed MPET system in the new unknowns $(u,p_{0},\tilde{p})$. For
the discussions below let us give remarks on the block operators defined by
$c_{1}$, $c_{2}$, $c_{3}$. Specifically, regarding $\alpha$ as a column
vector,

(4.13) ${C}_{2}=\left(\lambda^{-1}\alpha\right)\otimes
I_{W},\quad{C}_{3}=-\tau K\otimes\Delta+\left(S+\tau E+L\right)\otimes I_{W},$

with $K$ and $E$ as in (3.1), $L$ is the matrix
$L_{ij}=\lambda^{-1}\alpha_{i}\alpha_{j}$, $S$ is the diagonal matrix such
that its $j$-th entry is $s_{j}$, and $I_{W}$ is the identity (functional) on
$W$; we recall that $Q$ is the $J$-fold cartesian product of $W$. We will
first describe the transformed MPET equations for general coordinate
transformation, $P$. From the form of the transformed equations, we will
extract the conditions for $P$ that yield a system that is suitable for the
construction of parameter-robust block preconditioners.

Suppose we have an, fixed but otherwise arbitrary, invertible coordinate
transformation matrix $P\in\mathbb{R}^{J\times J}$. Applying this
transformation of variables to the semi-discretized total pressure variational
formulation of the MPET equations (4.4), we obtain the following variational
formulation: find the displacement $u\in V$, the total pressure $p_{0}\in
Q_{0}$ and the transformed pressures
$\tilde{p}=(\tilde{p}_{1},\dots,\tilde{p}_{J})\in Q$ such that

(4.14a) $\displaystyle a(u,v)+b(v,p_{0})$ $\displaystyle=\langle f,v\rangle$
$\displaystyle\quad\forall\,v\in V,$ (4.14b) $\displaystyle
b(u,q_{0})-c_{1}(p_{0},q_{0})-\tilde{c}_{2}(q_{0},\tilde{p})$
$\displaystyle=0$ $\displaystyle\quad\forall\,q_{0}\in Q_{0},$ (4.14c)
$\displaystyle-\tilde{c}_{2}(p_{0},\tilde{q})-\tilde{c}_{3}(\tilde{p},\tilde{q})$
$\displaystyle=\langle g,\left(P\otimes I_{W}\right)\tilde{q}\rangle$
$\displaystyle\quad\forall\,\tilde{q}\in Q$

where

(4.15) $\tilde{c}_{2}(q_{0},\tilde{q})\equiv c_{2}(q_{0},\left(P\otimes
I_{W}\right)\tilde{q}),\qquad\tilde{c}_{3}(\tilde{p},\tilde{q})\equiv
c_{3}(\left(P\otimes I_{W}\right)\tilde{p},\left(P\otimes
I_{W}\right)\tilde{q}).$

We define $\mathcal{\tilde{A}}_{\rm MPET}:V\times Q_{0}\times
Q\rightarrow(V\times Q_{0}\times Q)^{*}$ as the operator corresponding to the
bilinear form (4.14). The operator form of the transformed system (4.14) then
reads as:

(4.16) $\mathcal{\tilde{A}}_{\rm MPET}\begin{pmatrix}u\\\ p_{0}\\\
\tilde{p}\end{pmatrix}=\begin{pmatrix}f\\\ 0\\\
\tilde{g}\end{pmatrix},\quad\mathcal{\tilde{A}}_{\rm
MPET}=\begin{pmatrix}A&B^{T}&\textbf{0}\\\ B&-C_{1}&-\tilde{C_{2}}^{*}\\\
\textbf{0}&-\tilde{C_{2}}&-\tilde{C_{3}}\end{pmatrix},$

where $A=-2\operatorname{div}(\mu\varepsilon)$, $B=\operatorname{div}$ as
before, and $\tilde{g}=\left(P^{T}\otimes I_{W}\right)g$. By inserting (4.13)
and reordering, we note that

$\displaystyle\tilde{C_{2}}$
$\displaystyle=P^{T}C_{2}=(\lambda^{-1}P^{T}\alpha)\otimes I_{W},$
$\displaystyle\tilde{C_{3}}$
$\displaystyle=P^{T}C_{3}P=-\tau(P^{T}KP)\otimes\Delta+\left(P^{T}SP+P^{T}\left(\tau
E+L\right)P\right)\otimes I_{W}.$

For simplicity we will write

(4.17) $\tilde{\alpha}=\left(P^{T}\otimes I_{W}\right)\alpha.$

We now look to apply Lemma 2.1 with the choice of operators $S=K\otimes\Delta$
and $T=\left(\tau E+L\right)\otimes I_{W}$. The matrices $K$ and $\tau E+L$
satisfy the required conditions and, thus, there exists (c.f. Appendix A) an
invertible transformation $P$ simultaneously diagonalizing $K$ and $\tau E+L$
by congruence. That is, we have matrices $\tilde{K}$ and $\tilde{\Gamma}$
given by the formulas

(4.18) $\displaystyle\tilde{K}$
$\displaystyle=P^{T}KP=\operatorname{diag}(\tilde{K}_{1},\dots,\tilde{K}_{J}),$
(4.19) $\displaystyle\tilde{\Gamma}$ $\displaystyle=P^{T}\left(\tau
E+L\right)P=\operatorname{diag}(\tilde{\gamma}_{1},\dots,\tilde{\gamma}_{J}).$

We point out that the storage coefficients $\\{s_{j}\\}_{j=1}^{J}$ are not
involved in this simultaneous diagonalization process. This is critically
important in order to achieve a preconditioner that is parameter-robust, even
in the presence of vanishing storage coefficients. For future reference we
briefly note that

(4.20) $\tilde{\gamma}_{j}\geq\tilde{\alpha}_{j}^{2}/\lambda,$

follows from the definition of $\tilde{\Gamma}$ in (4.19) since $E$ is
positive semi-definite and therefore $\tilde{\gamma}_{j}$ is greater than or
equal to the $j$-th diagonal entry of the matrix
$\left(P^{T}LP\right)_{ij}=\lambda^{-1}\tilde{\alpha}_{i}\tilde{\alpha}_{j}$.

We also remark that the following identity holds for $\tilde{c}_{3}$:

(4.21) $\displaystyle\tilde{c}_{3}(\tilde{p},\tilde{q})$
$\displaystyle=\tau\sum_{j=1}^{J}\langle\tilde{K}_{j}\operatorname{\nabla}\tilde{p}_{j},\operatorname{\nabla}\tilde{q}_{j}\rangle+\sum_{j=1}^{J}\langle
s_{j}\left(\left(P\otimes
I_{W}\right)\tilde{p})\right)_{j},\left(\left(P\otimes
I_{W}\right)\tilde{q}\right)_{j}\rangle,$
$\displaystyle\quad+\tau\sum_{j=1}^{J}\sum_{i=1}^{J}\langle\xi_{j\leftarrow
i}(\left(\left(P\otimes I_{W}\right)\tilde{p}\right)_{j}-\left(\left(P\otimes
I_{W}\right)\tilde{p}\right)_{i}),\left(\left(P\otimes
I_{W}\right)\tilde{q}\right)_{j}\rangle$
$\displaystyle\quad+\langle\lambda^{-1}\tilde{\alpha}\cdot\tilde{p},\tilde{\alpha}\cdot\tilde{q}\rangle$

where $\left(\left(P\otimes I_{W}\right)\tilde{p}\right)_{j}$ is the $j$-th
component of $\left(P\otimes I_{W}\right)\tilde{p}$.

###  4.4 Preconditioning of the transformed MPET system

In this subsection we show that a parameter-robust preconditioner can be
constructed using an appropriate parameter-dependent norm.

We first define a parameter-dependent norm

(4.22)
$\|(u,p_{0},\tilde{p})\|_{\tilde{\mathcal{B}}}^{2}=\|\varepsilon(u)\|_{2\mu}^{2}+\|p_{0}\|_{(2\mu)^{-1}}^{2}+\sum_{j=1}^{J}\|\operatorname{\nabla}\tilde{p}_{j}\|_{\tau\tilde{K}_{j}}^{2}+\sum_{j=1}^{J}\|\tilde{p}_{j}\|_{\tilde{\gamma}_{j}}^{2}$

and consider the associated block preconditioner of the form

(4.23) $\mathcal{\tilde{B}}_{\rm
MPET}=\begin{pmatrix}(-2\mu\Delta)^{-1}&0&0&\cdots&0\\\ 0&2\mu
I^{-1}&0&\cdots&0\\\
0&0&(-\tau\tilde{K}_{1}\Delta+\tilde{\gamma}_{1}I)^{-1}&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&(-\tau\tilde{K}_{J}\Delta+\tilde{\gamma}_{J}I)^{-1}\end{pmatrix}.$

######  Lemma 4.2 (Continuity).

Let $\tilde{\mathcal{A}}_{\rm MPET}$ be defined by (4.16), and consider the
norm defined by (4.22). We assume that $2\mu\leq M_{0}\lambda$ for some
$M_{0}>0$ and $s_{j}\lesssim\tilde{\gamma}_{j}$ for $1\leq j\leq J$. Then
there exists a constant $C>0$, dependent on $M_{0}$, the constants in
$s_{j}\lesssim\tilde{\gamma}_{j}$, the matrix $P$, and the number of networks
$J$ but independent of any other problem parameters, such that

(4.24) $\langle\tilde{\mathcal{A}}_{\rm
MPET}(u,p_{0},\tilde{p}),(v,q_{0},\tilde{q})\rangle\leq
C\|(u,p_{0},\tilde{p})\|_{\tilde{\mathcal{B}}}\|(v,q_{0},\tilde{q})\|_{\tilde{\mathcal{B}}},$

for all $(u,p_{0},\tilde{p}),(v,q_{0},\tilde{q})\in V\times Q_{0}\times Q$.

######  Proof 4.3.

By redistributing the material parameter weights and the Cauchy-Schwarz
inequality, we obtain the preliminary upper bound

$\begin{split}\langle\tilde{\mathcal{A}}_{\rm
MPET}(u,p_{0},\tilde{p}),(v,q_{0},\tilde{q})\rangle\leq
Z_{1}+Z_{2}+Z_{3}=:Z,\end{split}$

where

$\displaystyle Z_{1}$
$\displaystyle=\|\varepsilon(u)\|_{2\mu}\|\varepsilon(v)\|_{2\mu}+\|p_{0}\|_{(2\mu)^{-1}}\|\operatorname{div}v\|_{2\mu}+\|q_{0}\|_{(2\mu)^{-1}}\|\operatorname{div}u\|_{2\mu},$
$\displaystyle Z_{2}$
$\displaystyle=\|p_{0}\|_{\lambda^{-1}}\|q_{0}\|_{\lambda^{-1}}+\|\tilde{\alpha}\cdot\tilde{q}\|_{\lambda^{-1}}\|p_{0}\|_{\lambda^{-1}}+\|\tilde{\alpha}\cdot\tilde{p}\|_{\lambda^{-1}}\|q_{0}\|_{\lambda^{-1}},$
$\displaystyle Z_{3}$
$\displaystyle=\operatorname{\textstyle\sum}_{j=1}^{J}\left(\|\operatorname{\nabla}\tilde{p}_{j}\|_{\tau\tilde{K}_{j}}\|\operatorname{\nabla}\tilde{q}_{j}\|_{\tau\tilde{K}_{j}}+\|\tilde{p}_{j}\|_{\tilde{\gamma}_{j}}\|\tilde{q}_{j}\|_{\tilde{\gamma}_{j}}+\langle
s_{j}(P\tilde{p})_{j},(P\tilde{q})_{j}\rangle\right).$

Since $\|\operatorname{div}u\|\leq\|\varepsilon(u)\|$ and by $2\mu\leq
M_{0}\lambda$ and the assumptions on $s_{j}$ and $\tilde{\gamma}_{j}$, it
follows that

$\begin{split}Z\lesssim&\left(\|\varepsilon(u)\|_{2\mu}+\|p_{0}\|_{(2\mu)^{-1}}+\|\tilde{\alpha}\cdot\tilde{p}\|_{\lambda^{-1}}+\operatorname{\textstyle\sum}_{j=1}^{J}\left(\|\operatorname{\nabla}\tilde{p}_{j}\|_{\tau\tilde{K}_{j}}+\|\tilde{p}_{j}\|_{\tilde{\gamma}_{j}}\right)\right)\\\
&\times\left(\|\varepsilon(v)\|_{2\mu}+\|q_{0}\|_{(2\mu)^{-1}}+\|\tilde{\alpha}\cdot\tilde{q}\|_{\lambda^{-1}}+\operatorname{\textstyle\sum}_{j=1}^{J}\left(\|\operatorname{\nabla}\tilde{q}_{j}\|_{\tau\tilde{K}_{j}}+\|\tilde{q}_{j}\|_{\tilde{\gamma}_{j}}\right)\right).\end{split}$

By the triangle inequality and (4.20) we obtain

(4.25)
$\|\tilde{\alpha}\cdot\tilde{p}\|_{\lambda^{-1}}\leq\sum_{j=1}^{J}\|\tilde{p}_{j}\|_{\tilde{\alpha}_{j}^{2}/\lambda}\leq\sum_{j=1}^{J}\|\tilde{p}_{j}\|_{\tilde{\gamma}_{j}}$

and it completes the proof.

######  Lemma 4.4 (Inf-sup condition).

Let $\tilde{\mathcal{A}}_{\rm MPET}$, $\tilde{\mathcal{B}}_{\rm MPET}$ and all
assumptions be as in Lemma 4.2. Then, there exists a constant $C>0$, dependent
on $M_{0}$, the constants in $s_{j}\lesssim\tilde{\gamma}_{j}$, and the number
of networks $J$ but independent of other material parameters, such that

(4.26)
$\adjustlimits{\inf}_{(u,p_{0},p)}{\sup}_{(v,q_{0},q)}\frac{\langle\tilde{\mathcal{A}}_{\rm
MPET}(u,p_{0},p),(v,q_{0},q)\rangle}{\|(u,p_{0},p)\|_{\tilde{\mathcal{B}}}\|(v,q_{0},q)\|_{\tilde{\mathcal{B}}}}\geq
C,$

where the $\inf$ and $\sup$ are taken over the non-vanishing elements in
$V\times Q_{0}\times Q$.

######  Proof 4.5.

Consider any $(u,p_{0},\tilde{p})\in V\times Q_{0}\times Q$, and choose
$\tilde{q}=-\tilde{p}$, and $q_{0}=-p_{0}$. Let $w\in V$ satisfy

(4.27)
$\langle\operatorname{div}w,p_{0}\rangle=\|p_{0}\|_{(2\mu)^{-1}}^{2},\quad\|\varepsilon(w)\|_{2\mu}\leq
C_{0}\|p_{0}\|_{(2\mu)^{-1}}.$

for a $C_{0}>0$ depending on the domain $\Omega$ via Korn’s inequality, and
next choose $v=u+2\delta w$ for $\delta>0$ to be further specified. We note
that, with this choice of $v,q_{0}$, and $q$,

$\|(v,q_{0},\tilde{q})\|_{\tilde{\mathcal{B}}}\lesssim\|(u,p_{0},\tilde{p})\|_{\tilde{\mathcal{B}}},$

with inequality constant depending only on the domain $\Omega$ and the choice
of $\delta$ since

$\|\varepsilon(v)\|_{2\mu}\leq\|\varepsilon(u)\|_{2\mu}+2\delta
C_{0}\|p_{0}\|_{(2\mu)^{-1}}\lesssim\|(u,p_{0},\tilde{p})\|_{\tilde{\mathcal{B}}}.$

Therefore, it suffices to show that

(4.28) $\langle\tilde{\mathcal{A}}_{\rm
MPET}(u,p_{0},\tilde{p}),(v,q_{0},\tilde{q})\rangle\gtrsim\|(u,p_{0},\tilde{p})\|_{\tilde{\mathcal{B}}}^{2}.$

Using the definition of $\tilde{\mathcal{A}}_{\rm MPET}$ together with (4.27),
we find that

(4.29) $\langle\tilde{\mathcal{A}}_{\rm
MPET}(u,p_{0},\tilde{p}),(v,q_{0},\tilde{q})\rangle=\|\varepsilon(u)\|^{2}_{2\mu}+2\delta\langle\varepsilon(u),\varepsilon(w)\rangle_{2\mu}+2\delta\|p_{0}\|^{2}_{(2\mu)^{-1}}\\\
+c_{1}(p_{0},p_{0})+2\tilde{c}_{2}(p_{0},\tilde{p})+\tilde{c}_{3}(\tilde{p},\tilde{p}).$

Note that $c_{1}(p_{0},p_{0})=\|p_{0}\|^{2}_{\lambda^{-1}}$,
$\tilde{c}_{2}(p_{0},\tilde{p})=\langle\tilde{\alpha}\cdot\tilde{p},p_{0}\rangle_{\lambda^{-1}}$,
and

$\displaystyle\tilde{c}_{3}(\tilde{p},\tilde{p})=\operatorname{\textstyle\sum}_{j=1}^{J}(\|\operatorname{\nabla}\tilde{p}_{j}\|^{2}_{\tau\tilde{K}_{j}}+\|\tilde{p}_{j}\|^{2}_{\tilde{\gamma}_{j}}+\langle
s_{j}(P\tilde{p})_{j},(P\tilde{p})_{j}\rangle)$

from the definitions of $\tilde{c}_{2}$, $\tilde{c}_{3}$, and the congruent
diagonalization. Moreover, we have

(4.30)
$\displaystyle\operatorname{\textstyle\sum}_{j=1}^{J}\|\tilde{p}_{j}\|^{2}_{\tilde{\gamma}_{j}}\geq\|\tilde{\alpha}\cdot\tilde{p}\|^{2}_{\lambda^{-1}}$

by observing the identity (4.21). Then we obtain

$\displaystyle\delta\|p_{0}\|_{(2\mu)^{-1}}^{2}+c_{1}(p_{0},p_{0})+2\tilde{c}_{2}(p_{0},\tilde{p})+\tilde{c}_{3}(\tilde{p},\tilde{p})$
$\displaystyle\geq\delta\|p_{0}\|_{(2\mu)^{-1}}^{2}+\|p_{0}\|_{\lambda^{-1}}^{2}+2\langle\tilde{\alpha}\cdot\tilde{p},p_{0}\rangle_{\lambda^{-1}}+\operatorname{\textstyle\sum}_{j=1}^{J}(\|\tilde{p}_{j}\|_{\tilde{\gamma}_{j}}^{2}+\|\operatorname{\nabla}\tilde{p}_{j}\|^{2}_{\tau\tilde{K}_{j}})$
$\displaystyle\geq\left({\delta}/{M_{0}}+1\right)\|p_{0}\|_{\lambda^{-1}}^{2}+2\langle\tilde{\alpha}\cdot\tilde{p},p_{0}\rangle_{\lambda^{-1}}+\operatorname{\textstyle\sum}_{j=1}^{J}(\|\tilde{p}_{j}\|_{\tilde{\gamma}_{j}}^{2}+\|\operatorname{\nabla}\tilde{p}_{j}\|^{2}_{\tau\tilde{K}_{j}})$
$\displaystyle=\|\sqrt{{\delta}/{M_{0}}+1}p_{0}+\tilde{\alpha}\cdot\tilde{p}\|_{\lambda^{-1}}^{2}-\left(\frac{\delta}{M_{0}}+1\right)^{-1}\|\tilde{\alpha}\cdot\tilde{p}\|_{\lambda^{-1}}^{2}$
$\displaystyle\quad+\operatorname{\textstyle\sum}_{j=1}^{J}(\|\tilde{p}_{j}\|^{2}_{\tilde{\gamma}_{j}}+\|\operatorname{\nabla}\tilde{p}_{j}\|^{2}_{\tau\tilde{K}_{j}})$
$\displaystyle\geq\operatorname{\textstyle\sum}_{j=1}^{J}\left(\frac{\delta}{M_{0}+\delta}\|\tilde{p}_{j}\|^{2}_{\tilde{\gamma}_{j}}+\|\operatorname{\nabla}\tilde{p}_{j}\|^{2}_{\tau\tilde{K}_{j}}\right)$

where the second inequality in the above follows from the assumption $2\mu\leq
M_{0}\lambda$ and the third inequality follows from (4.30). Thus,

$\displaystyle\langle\tilde{\mathcal{A}}_{\rm
MPET}(u,p_{0},\tilde{p}),(v,q_{0},\tilde{q})\rangle$ (4.31)
$\displaystyle\quad\geq\|\varepsilon(u)\|^{2}_{2\mu}+\delta\langle\varepsilon(u),\varepsilon(w)\rangle_{2\mu}+\delta\|p_{0}\|^{2}_{(2\mu)^{-1}}+\operatorname{\textstyle\sum}_{j=1}^{J}\|\operatorname{\nabla}\tilde{p}_{j}\|^{2}_{\tau\tilde{K}_{j}}$
$\displaystyle\qquad+\frac{\delta}{M_{0}+\delta}\operatorname{\textstyle\sum}_{j=1}^{J}\|\tilde{p}_{j}\|^{2}_{\tilde{\gamma}_{j}}.$

On the other hand, the Cauchy-Schwarz inequality, the definition of $w$, and
Young’s inequality give that

(4.32) $\delta|\langle\varepsilon(u),\varepsilon(w)\rangle_{2\mu}|\leq\delta
C_{0}\|\varepsilon(u)\|_{2\mu}\|p_{0}\|_{(2\mu)^{-1}}\leq\frac{1}{2}\|\varepsilon(u)\|^{2}_{2\mu}+\frac{1}{2}\delta^{2}C_{0}^{2}\|p_{0}\|^{2}_{(2\mu)^{-1}}.$

Inserting the negation of (4.32) as a lower bound in (4.5), we thus obtain
that

$\langle\tilde{\mathcal{A}}_{\rm
MPET}(u,p_{0},\tilde{p}),(v,q_{0},\tilde{q})\rangle\\\
\geq\frac{1}{2}\|\varepsilon(u)\|^{2}_{2\mu}+\delta(1-\frac{1}{2}\delta
C_{0}^{2})\|p_{0}\|^{2}_{(2\mu)^{-1}}+\operatorname{\textstyle\sum}_{j=1}^{J}\left(\|\operatorname{\nabla}\tilde{p}_{j}\|^{2}_{\tau\tilde{K}_{j}}+\frac{\delta}{M_{0}+\delta}\|\tilde{p}_{j}\|^{2}_{\tilde{\gamma}_{j}}\right).$

By choosing $\delta$, in particular e.g. by letting $\delta<2/C_{0}^{2}$, the
estimate (4.28) follows.

######  Remark 4.6.

In Lemma 4.2 and Lemma 4.4 we assumed $s_{j}\lesssim\tilde{\gamma}_{j}$ for
$1\leq j\leq J$ and it covers the cases that $s_{j}$’s are degenerate. This
assumption can be removed if $\lambda^{-1}\lesssim s_{j}$, $1\leq j\leq J$,
hold with constants of scale 1. For parameter-robust preconditioners we use
$P$ which gives a different simultaneous diagonalization. More precisely, we
consider $P$ satisfying (4.18), and $P^{T}(S+\tau
E+L)P=\operatorname{diag}(\tilde{\gamma}_{1},\dots,\tilde{\gamma}_{J})$
instead of (4.19). The norm (4.22) with these new $\tilde{K}_{j}$’s and
$\tilde{\gamma}_{j}$’s, will be used to obtain parameter-robust
preconditioners. Since the modification of proofs is straightforward and most
steps are almost same, we omit the detailed proofs.

For concreteness, we here illustrate the form of the MPET equations and of the
proposed preconditioner in a specific example.

######  Example 4.7.

We consider the simple case of two networks with $K_{1}=K_{2}=1.0$,
$s_{1}=s_{2}=1.0$, $\alpha_{1}=\alpha_{2}=0.5$, $\lambda=1.0$,
$\xi_{1\rightarrow 2}=0.0$, and $\tau=1.0$. The transformation matrix in this
case is

(4.33) $P=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\\ 1&1\end{pmatrix}.$

We remark that $P$ is not normalized. The associated transformed MPET operator
(expanded), cf. (4.16) and associated definitions, is then

(4.34) $\tilde{\mathcal{A}}_{\rm
MPET}=\begin{pmatrix}-2\mu\operatorname{div}\varepsilon&-\operatorname{\nabla}&0&0\\\
\operatorname{div}&-\lambda^{-1}&-(\sqrt{2}\lambda)^{-1}&0\\\
0&-(\sqrt{2}\lambda)^{-1}&-\Delta+\frac{3}{2}&0\\\ 0&0&0&-\Delta+1\\\
\end{pmatrix},$

and the proposed preconditioner will be in the following form:

(4.35) $\mathcal{\tilde{B}}_{\rm
MPET}=\begin{pmatrix}(-2\mu\operatorname{div}\varepsilon)^{-1}&0&0&0\\\
0&(2\mu)^{-1}&0&0\\\ 0&0&(-\Delta+1)^{-1}&0\\\
0&0&0&(-\Delta+1)^{-1}\end{pmatrix}.$

The objective of this example was to illustrate the layout of the operators in
a simple case. The results for more general numerical examples will be
presented later.

###  4.5 Numerical performance

######  Example 4.8.

In this example we demonstrate the robustness of the block diagonal
preconditioner (4.23) for a mixed finite element discretization of the
transformed total pressure MPET equations (4.16). We consider the same test
case, discretization and solver set-up as described in Example 4.1; the new
preconditioner is the only modification. Parameter ranges are as follows:
$K_{2}\in[10^{-6},10^{6}]$, $\xi_{1\leftarrow 2}\in[10^{-6},10^{6}]$ and
$\lambda\in[1,10^{6}]$.

The resulting number of iterations are shown in Figure 1 for
$K_{2}\in[10^{-6},1]$ and $\xi_{1\leftarrow 2}\in[1,10^{6}]$; omitted values
demonstrated similar behaviour. Each of the subplots in Figure 1 represent a
fixed choice of $K_{2}$ and $\xi_{1\leftarrow 2}$. In each subplot four curves
are shown; these curves show the number of MinRes iterations corresponding to
different values of $\lambda$, indicated by their respective symbols, at
discretization levels $N=16$, $32$, $64$ and $128$. The stopping criterion was

$(\tilde{\mathcal{B}}r_{k},r_{k})/(\tilde{\mathcal{B}}r_{0},r_{0})\leq
10^{-6}$

where $r_{k}$ is the residual of the $k$-th iteration. We observe that the
number of iterations is moderate in general. Moreover, the number of
iterations does not grow for smaller $K_{2}$’s relative to larger
$\xi_{1\leftarrow 2}$ or larger $N$ – in contrast to what was observed for
Example 4.1.

Figure 1: Number of MinRes iterations: (4.4) discretized with Taylor-Hood type
elements and algebraic multigrid, for $s_{1}=s_{2}=1.0$. $K_{2}$ varies along
the horizontal axis while the vertical axis shows variations in
$\xi_{1\leftarrow 2}$ for $K_{2}$ fixed. Each subplot contains four piecewise
linear curves; each curve is decorated by a symbol indicating a corresponding
value of $\lambda$ and corresponds to results for discretizations $N=16,32,64$
and $128$. Figure 2: Number of MinRes iterations: (4.4) discretized with
Taylor-Hood type elements and algebraic multigrid for $s_{1}=s_{2}=0$. $K_{2}$
varies along the horizontal axis while the vertical axis shows variations in
$\xi_{1\leftarrow 2}$ for $K_{2}$ fixed. Each subplot contains four piecewise
linear curves; each curve is decorated by a symbol indicating a corresponding
value of $\lambda$ and corresponds to results for discretizations $N=16,32,64$
and $128$.

######  Example 4.9.

In this final example we present a modified version of a 3D footing problem
[7, 24, 26]; we demonstrate the problem for two fluid networks ($J=2$) and use
the standard unit cube, $\Omega=[0,1]^{3}\subset{\mathbb{R}}^{3}$, as
computational domain. At the base of the domain, homogeneous Dirichlet
conditions for the displacement and for both fluid pressures are imposed.At
the top-most surface of the domain, i.e. $z=1$, a load of $0.1N/m^{2}$ is
applied on the square $[0.25,0.75]\times[0.25,0.75]$, and a no flow condition
is applied to the fluid pressures. For all remaining boundary sides of the
domain, the zero stress condition is applied alongside a homogeneous Dirichlet
condition for the fluid pressures. In the numerical experiments we vary the
exchange coefficient $\xi_{1\leftarrow 2},$ and the mesh size, the other
physical parameters are reported in Table 3.

(a) fluid pressure $p_{1}$

(b) fluid pressure $p_{2}$

(c) displacement magnitude $u$

Figure 3: A 3D Footing Problem, solution for $1/h=32,$ $\xi_{1\leftarrow 2}=1.0,$ at $t=0.5$ Property | Symbol | Value | Units
---|---|---|---
Young’s modulus | $E$ | $3\times 10^{4}$ | Pa
Poisson ratio | $\nu$ | 0.45 | [-]
Hydraulic conductivities | $K_{1},K_{2}$ | $10^{-6}$ |  m2(Pa s)-1
Storage coefficients | $s_{1},s_{2}$ | $0.0$ |  Pa-1
Biot coefficient | $\alpha_{1},\alpha_{2}$ | $0.5$ | [-]
Table 3: Parameters used in the numerical simulations

In table 4, we report the number of MinRes iterations for each time-step (from
$0.1$ to $0.5$), varying the mesh size and exchange parameters. The initial
guess for the solution is set to zero. Similarly to what observed in Example
4.8, the number of iterations is moderate also for this 3D case. In Fig. 3 the
solution for $t=0.5$ is shown.

h | $\xi_{1\leftarrow 2}$ | Number of iterations
---|---|---
$t=0.1$ | $t=0.2$ | $t=0.3$ | $t=0.4$ | $t=0.5$
$1/8$ | $1.0\times 10^{-6}$ | $87$ | $97$ | $97$ | $97$ | $97$
$1.0$ | $89$ | $102$ | $102$ | $102$ | $102$
$1/16$ | $1.0\times 10^{-6}$ | $90$ | $102$ | $102$ | $102$ | $102$
$1.0$ | $93$ | $108$ | $109$ | $107$ | $109$
$1/32$ | $1.0\times 10^{-6}$ | $95$ | $107$ | $107$ | $107$ | $107$
$1.0$ | $98$ | $112$ | $112$ | $114$ | $111$
Table 4: MinRes iterations for the footing problem (c.f. Example 4.9).

##  5 Conclusions

In this paper, we have presented a new strategy for decoupling the total-
pressure variational formulation of the multiple-network poroelasticity
equations. The decoupling strategy is based on a transformation via a change
of variables, allowing for simultaneous diagonalization by congruence of the
equation operators. In particular, the transformed equations are readily
amenable for block–diagonal preconditioning. Moreover, we have proposed a
block-diagonal preconditioner for the transformed system and shown
theoretically that the preconditioner and the equation operator are norm
equivalent, independently of the material parameters, under reasonable
parameter assumptions. The theoretical results are supported by numerical
examples. Combined, these results allow the efficient iterative solution of
the multiple-network poroelasticity equations, even in the case of nearly
incompressible materials.

We note that our strategy is based on spatially constant material parameters.
The applicability of this approach for spatially varying parameters has not
yet been considered.

##  Appendix A Proof of Lemma 2.1

We first recall a basic [12] definition and result for posterity.

######  Definition A.1.

A matrix $C\in\mathbb{C}^{n\times n}$ is diagonalizable if there exists an
invertible transformation, $P$, such that $P^{-1}CP$ is diagonal. The matrix
$C$ is called diagonalizable by congruence if there exists $P$, not
necessarily invertible, such that $P^{T}CP$ is diagonal.

######  Theorem 0 (4.5.17a-b p. 287, [12]).

Suppose $A$ and $B\in\mathbb{C}^{n\times n}$ are symmetric and that $A$ is
invertible. Then $A$ and $B$ are diagonalizable by congruence if and only if
$C=A^{-1}B$ is diagonalizable.

######  Proof A.3 (Proof of Lemma 2.1).

Assume $K$ and $E$ satisfy the hypotheses of Lemma 2.1; we first show that
$C=K^{-1}E\in\mathbb{R}^{n\times n}$ is diagonalizable. We note that $C$
satisfies

$K^{1/2}CK^{-1/2}=K^{-1/2}EK^{-1/2}.$

The right-hand side, above, is symmetric due to the symmetry of $K$ and $E$.
Thus $C$ is similar to a real, symmetric matrix and is therefore
diagonalizable. From [12, 4.5.17a-b] there exists an invertible matrix
$P\in\mathbb{R}^{n\times n}$ such that

$P^{T}KP=\tilde{D}_{K}\quad\text{ and }\quad P^{T}EP=\tilde{D}_{E},$

where $\tilde{D}_{K}$, $\tilde{D}_{E}\in\mathbb{R}^{n\times n}$ are diagonal
matrices.

Recalling $Q=W\times\cdots\times W$, define the change of variables
$\tilde{p}=(P^{-1}\otimes I_{W})p$ for $p\in Q$ and substitute into (2.9) to
get

$\langle S(P\otimes I_{W})\tilde{p},q\rangle+\langle T(P\otimes
I_{W})\tilde{p},q\rangle=\langle f,q\rangle,\quad\forall\,q\in Q.$

Writing $q=(P\otimes I_{W})(P^{-1}\otimes I_{W})q$ and noting that the adjoint
operator of $P\otimes I_{W}$ is $P^{T}\otimes I_{W^{*}}$, we have

$\langle S(P\otimes I_{W})\tilde{p},q\rangle=\langle(P^{T}\otimes
I_{W^{*}})S(P\otimes I_{W})\tilde{p},(P^{-1}\otimes I_{W})q\rangle.$

Since $S=K\otimes A$ we can obtain $(P^{T}\otimes I_{W^{*}})S(P\otimes
I_{W})=\tilde{D}_{K}\otimes A$ (by Hadamard product). By a similar argument

$\langle TP\tilde{p},q\rangle=\langle(\tilde{D}_{E}\otimes
B)\tilde{p},(P^{-1}\otimes I_{W})q\rangle.$

Then $D_{S}:=\tilde{D}_{K}\otimes A$ and $D_{T}:=\tilde{D}_{E}\otimes B$ are
block diagonal operators from $Q$ to $Q^{*}$. Finally $\langle
f,q\rangle=\langle(P^{T}\otimes I_{W^{*}})f,(P^{-1}\otimes I_{W})q\rangle$ and
the variational problem (2.10) follows because $q\in Q$ is arbitrary.

######  Remark A.4.

The construction of the matrix $P$, yielding both $P^{T}AP=D_{1}$ and
$P^{T}BP=D_{2}$, is straightforward for the case when $C=A^{-1}B$ has $n$
distinct eigenvalues. In this case $C$ has $n$ linearly independent
eigenvectors; if $\\{v_{1},\dots,v_{n}\\}$ denote these eigenvectors then
$P=[v_{1},\dots,v_{n}]$ is the matrix whose $j$-th column is $v_{j}$. When the
eigenvalues of $C$ are not distinct: $P$ can be realized as the product of
block-wise eigenvector matrices. The general procedure for this case is
discussed in [12]; an example has been discussed in Section 3.

## References

  * [1] M. Bai, D. Elsworth, and J.-C. Roegiers, Multiporosity/multipermeability approach to the simulation of naturally fractured reservoirs, Water Resources Research, 29 (1993), pp. 1621–1633.
  * [2] M. A. Biot, General theory of three-dimensional consolidation, Journal of Applied Physics, 12 (1941), pp. 155–164.
  * [3] J. Brašnová, V. Lukeš, and E. Rohan, Identification of multi-compartment Darcy flow model material parameters, in 20th International Conference Applied Mechanics 2018: April 9-11, 2018, Myslovice, Czech republic: conference proceedings, 2019, pp. 9–13.
  * [4] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers, Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique, 8 (1974), pp. 129–151.
  * [5] D. Chou, J. C. Vardakis, L. Guo, B. J. Tully, and Y. Ventikos, A fully dynamic multi-compartmental poroelastic system: Application to aqueductal stenosis, Journal of Biomechanics, 49 (2016), pp. 2306–2312.
  * [6] A. Eisenträger and I. Sobey, Multi-fluid poroelastic modelling of CSF flow through the brain, in Poromechanics V: Proceedings of the Fifth Biot Conference on Poromechanics, 2013, pp. 2148–2157.
  * [7] F. Gaspar, J. Gracia, F. Lisbona, and C. Oosterlee, Distributive smoothers in multigrid for problems with dominating grad-div operators, Numerical Linear Algebra with Applications, 15 (2008), pp. 661–683, https://doi.org/10.1002/nla.587.
  * [8] Q. Hong and J. Kraus, Parameter-robust stability of classical three-field formulation of Biot’s consolidation model, Electron. Trans. Numer. Anal., 48 (2018), pp. 202–226, https://doi.org/10.1553/etna_vol48s202, https://doi.org/10.1553/etna_vol48s202.
  * [9] Q. Hong, J. Kraus, M. Lymbery, and F. Philo, Conservative discretizations and parameter-robust preconditioners for Biot and multiple-network flux-based poroelasticity models, Numerical Linear Algebra with Applications, 26 (2019), p. e2242.
  * [10] Q. Hong, J. Kraus, M. Lymbery, and F. Philo, Parameter-robust Uzawa-type iterative methods for double saddle point problems arising in Biot’s consolidation and multiple-network poroelasticity models, arXiv preprint arXiv:1910.05883, (2019).
  * [11] Q. Hong, J. Kraus, M. Lymbery, and M. F. Wheeler, Parameter-robust convergence analysis of fixed-stress split iterative method for multiple-permeability poroelasticity systems, Multiscale Model. Simul., 18 (2020), pp. 916–941, https://doi.org/10.1137/19M1253988, https://doi.org/10.1137/19M1253988.
  * [12] R. A. Horn and C. R. Johnson, Matrix analysis, Cambridge university press, 2nd ed., 1990.
  * [13] T. Józsa, W. El-Bouri, R. Padmos, S. Payne, and A. Hoekstra, A cerebral circulation model for in silico clinical trials of ischaemic stroke.  In the proceedings of CompBioMed Conference 2019.
  * [14] J. Lee, A. Cookson, R. Chabiniok, S. Rivolo, E. Hyde, M. Sinclair, C. Michler, T. Sochi, and N. Smith, Multiscale modelling of cardiac perfusion, in Modeling the heart and the circulatory system, Springer, 2015, pp. 51–96.
  * [15] J. J. Lee, Unconditionally stable second order convergent partitioned methods for multiple-network poroelasticity, arXiv preprint arXiv:1901.06078, (2019).
  * [16] J. J. Lee, K.-A. Mardal, and R. Winther, Parameter-robust discretization and preconditioning of Biot’s consolidation model, SIAM Journal on Scientific Computing, 39 (2017), pp. A1–A24.
  * [17] J. J. Lee, E. Piersanti, K.-A. Mardal, and M. E. Rognes, A mixed finite element method for nearly incompressible multiple-network poroelasticity, SIAM Journal on Scientific Computing, 41 (2019), pp. A722–A747.
  * [18] K.-A. Mardal and R. Winther, Preconditioning discretizations of systems of partial differential equations, Numerical Linear Algebra with Applications, 18 (2011), pp. 1–40.
  * [19] C. Michler, A. Cookson, R. Chabiniok, E. Hyde, J. Lee, M. Sinclair, T. Sochi, A. Goyal, G. Vigueras, D. Nordsletten, et al., A computationally efficient framework for the simulation of cardiac perfusion using a multi-compartment Darcy porous-media flow model, International Journal for Numerical Methods in Biomedical Engineering, 29 (2013), pp. 217–232.
  * [20] A. Mikelić and M. F. Wheeler, Convergence of iterative coupling for coupled flow and geomechanics, Comput. Geosci., 17 (2013), pp. 455–461, https://doi.org/10.1007/s10596-012-9318-y, https://doi.org/10.1007/s10596-012-9318-y.
  * [21] M. J. M. Mokhtarudin and S. J. Payne, The study of the function of AQP4 in cerebral ischaemia–reperfusion injury using poroelastic theory, International Journal for Numerical Methods in Biomedical Engineering, 33 (2017), p. e02784.
  * [22] R. Oyarzúa and R. Ruiz-Baier, Locking-free finite element methods for poroelasticity, SIAM Journal on Numerical Analysis, 54 (2016), pp. 2951–2973.
  * [23] E. Piersanti, M. E. Rognes, and K.-A. Mardal, Parameter robust preconditioning for multi–compartmental Darcy equations, in , 2019.  (in review for EnuMath 2019).
  * [24] C. Rodrigo, X. Hu, P. Ohm, J. Adler, F. Gaspar, and L. Zikatanov, New stabilized discretizations for poroelasticity and the stokes’ equations, Comput. Methods Appl. Mech. Engrg., 341 (2018), pp. 467–484, https://doi.org/10.1016/j.cma.2018.07.003.
  * [25] R. Shipley, P. Sweeney, S. Chapman, and T. Roose, A four-compartment multiscale model of fluid and drug distribution in vascular tumours, International Journal for Numerical Methods in Biomedical Engineering, (2020), p. e3315.
  * [26] E. Storvik, J. Both, K. Kuman, J. Nordbotten, and F. Radu, On the optimization of the fixed-stress splitting for biot’s equations, Int. J. Numer. Methods. Eng., 120 (2019), pp. 179–194, https://doi.org/10.1002/nme.6130.
  * [27] K. H. Støverud, M. Alnæs, H. P. Langtangen, V. Haughton, and K.-A. Mardal, Poro-elastic modeling of syringomyelia–a systematic study of the effects of pia mater, central canal, median fissure, white and gray matter on pressure wave propagation and fluid movement within the cervical spinal cord, Computer Methods in Biomechanics and Biomedical Engineering, 19 (2016), pp. 686–698.
  * [28] B. Tully and Y. Ventikos, Modelling normal pressure hydrocephalus as a ‘two-hit’disease using multiple-network poroelastic theory, in ASME 2010 Summer Bioengineering Conference, American Society of Mechanical Engineers Digital Collection, 2013, pp. 877–878.
  * [29] B. J. Tully and Y. Ventikos, Cerebral water transport using multiple-network poroelastic theory: application to normal pressure hydrocephalus, Journal of Fluid Mechanics, 667 (2011), pp. 188–215.
  * [30] J. C. Vardakis, D. Chou, B. J. Tully, C. C. Hung, T. H. Lee, P.-H. Tsui, and Y. Ventikos, Investigating cerebral oedema using poroelasticity, Medical Engineering & Physics, 38 (2016), pp. 48–57.