# Efficient reference configuration formulation in fully nonlinear poroelastic
media

Nicolás A. Barnafi  Argyrios Petras  Luca Gerardo-Giorda  Center for
Mathematical Modelling (CMM), Universidad de Chile, Beauchef 851, Santiago,
Chile<EMAIL_ADDRESS>Johann Radon Institute for Computational and
Applied Mathematics (RICAM), Austrian Academy of Sciences, Linz, Austria
Institute for Mathematical Methods in Medicine and Data-Based Modelling,
Johannes Kepler University, Linz, Austria

###### Abstract

Typical pipelines for model geometry generation in computational biomedicine
stem from images, which are usually considered to be at rest, despite the
object being in mechanical equilibrium under several forces. We refer to the
stress-free geometry computation as the reference configuration problem, and
in this work we extend such a formulation to the theory of fully nonlinear
poroelastic media. The main steps are (i) writing the equations in terms of
the reference porosity and (ii) defining a time dependent problem whose steady
state solution is the reference porosity. This problem can be computationally
challenging as it can require several hundreds of iterations to converge, so
we propose the use of Anderson acceleration to speed up this procedure. Our
evidence shows that this strategy can reduce the number of iterations up to
80%. In addition, we note that a primal formulation of the nonlinear mass
conservation equations is not consistent due to the presence of second order
derivatives of the displacement, which we alleviate through adequate mixed
formulations. All claims are validated through numerical simulations in both
idealized and realistic scenarios.

###### keywords:

Nonlinear poroelasticity , Anderson acceleration , Reference configuration ,
Inverse design , Robust discretizations , Mixed formulations

##  1 Introduction

Porous media, understood as objects that exhibit a solid phase and an
interconnected network through which a fluid flows, are ubiquitous in nature
and can be found in materials like wood, rubber, and stones, in animals such
as sponges, and in organs like lungs and hearts. See [10] for further examples
and applications. These objects may deform, but they preserve their porous
structure under the action of forces. These forces can be either external or
internal from within the pores, and in fact internal pressures in organs play
a significant role in their function. For instance, when a baby takes their
first breath, it generates a pressure gradient that pushes the fluid into the
interstitial tissue, allowing it to be removed by pulmonary and lymphatic
circulations [21]. The porous structure of biological tissues is essential not
only for the normal function of organs through blood perfusion, but also for
drug delivery and pharmacological treatment of various diseases [2].

Biomedical models typically consider biological tissues as elastic continua,
disregarding their complex vascular network that provides blood, which is
essential for their function. However, the extreme complexity of vascular
networks has motivated the use of a porous media modeling approach in which
only the average macroscopic behavior of flow in such networks is described
[19]. Most models use multiple compartments to account for the blood vessels
(arteries, veins, and capillaries), the solid part, and even the interstitial
space of tissues undergoing large deformations, which applies mainly to the
lungs [25] and the heart [13, 12, 11, 7].

Biological tissues and organs are subject to external forces (gravity, rib
cage) and internal forces, which can be either given by air pressure (lungs)
or fluid pressure (vessels, lymphatic system, extracellular space), among
others. These forces are already present when a patient’s geometry is computed
from medical imaging, so it becomes fundamental to find a reference
configuration whose deformed state is the one computed. This problem, known as
the inverse design problem [17] or prestress problem [18], has been
successfully applied to the case of a beating (purely elastic) heart [28, 23,
24, 9].

To the best of our knowledge, porous media models using a prestress
configuration have only been developed in [25] under strong assumptions that
remove the mass conservation law from the system. This work has been motivated
by this missing reference configuration model that accommodates to arbitrary
nonlinear porous media undergoing large deformations. For this, we propose a
novel poroelastic media equation that yields the reference porosity, together
with the reference configuration, as an unknown initial condition.

The resulting model presents several difficulties beyond being extremely
nonlinear. The main two are: (i) it can take many time iterations to yield a
reference configuration, and (ii) it presents second order derivatives of the
displacement–an $\vectorsym H^{1}$ field–which we refer to as the _primal
inconsistency_ error. To address these limitations, we propose for (i) to use
a fixed-point acceleration technique to reduce the number of iterations, and
for (ii) to formulate mixed models that do not suffer from primal
inconsistency. The resulting models are accurate and present a framework that
is also compatible with their solution in extreme scale computing. We validate
our claims through several numerical tests.

This work is structured as follows: in Section 2 we present both the forward
equation and the proposed novel reference configuration formulation. In
Section 3 we present all details pertaining the discretization of the obtained
model in order to have a robust and efficient solution strategy. In this
section we also present how the computation of a stationary state can be
interpreted as a fixed-point problem, and how this can be accelerated. In
Section 4, we present our mixed models that avoid the primal consistency
error, and comment on their theoretical properties. We present several
numerical tests in Section 5, and conclude with a discussion of our results in
Section 6.

To avoid confusing notation, throughout this work we will refer to the direct
problem as the _forward_ problem. Instead, we will refer to the proposed
system as the reference configuration model. Also, we will denote scalar,
vector, and tensor quantities with regular, bold, and slanted fonts such that
$L^{2}$, $\vectorsym L^{2}$, and $\tensorsym L^{2}$ are the corresponding
square integrable spaces.

##  2 Model formulation

In this section we show the formulation of the forward and reference
configuration models of nonlinear poroelasticity. We provide all the balance
laws related to nonlinear poroelasticity in both Lagrangian and Eulerian
coordinates, as both configurations are relevant for devising the forward and
reference configuration problems respectively. For further explanation on the
model, see [15].

Consider a reference domain $\Omega_{0}$ and a displacement field $\vectorsym
d:\Omega_{0}\to\mathbb{R}^{3}$ such that a reference (Lagrangian) point
$\vectorsym X$ is warped into a spatial (current, Eulerian) one $\vectorsym x$
through $\vectorsym x(\vectorsym X)=\vectorsym X+\vectorsym d(\vectorsym X)$.
We denote the deformed configuration as $\Omega=\vectorsym x(\Omega_{0})$, and
further define the strain tensor and its determinant as $\tensorsym
F\coloneqq\tensorsym I+\operatorname{\nabla}_{\vectorsym X}\vectorsym d$ and
$J\coloneqq\det\tensorsym F$ respectively. As done with the definition of
$\tensorsym F$, whenever relevant we will clarify the coordinates with respect
to which we are deriving as $\operatorname{\nabla}_{\vectorsym X}$ and
$\operatorname{\nabla}_{\vectorsym x}$. The reference and current outwards
normal vector will be denoted with $\vectorsym N$ and $\vectorsym n$
respectively.

The geometry $\Omega_{0}$ is assumed to be a porous medium, where fluid flows
through a complex interconnected network structure. Being a continuum theory,
both solid and fluid phases coexist spatially. The percentage of the fluid
phase in $\Omega$ is given by the Eulerian porosity $\phi$, so that naturally
the solid porosity is given by $\phi_{s}=1-\phi$. The Lagrangian porosity
$\Phi$ is defined as $\Phi=J\phi$, such that

$\int_{\omega_{0}}\Phi\,dX=\int_{\omega}\phi\,dx,$

for all $\omega_{0}\subset\Omega_{0}$, $\omega=\vectorsym x(\omega_{0})$. Note
that all developments and conclusions given in this work can be easily
extended to the case of more than one fluid phase, known as _multi-
compartment_ of _multi-phase_ fluids.

###  2.1 Forward problem

In this section, we recapitulate the balance laws for the forward model, given
by a quasi-static poroelastic solid. These are: (i) conservation of linear
momentum, (ii) mass conservation, and (iii) solid incompressibility (see [15]
for their derivations). We highlight that incompressibility is not strictly
necessary for this formulation, but we nevertheless decided to include it.
This is motivated by our target applications and also by the fact that it can
be easily approximated using quasi-incompressibility or even removed, thus
rendering the methodology more general.

#### Equilibrium equation.

Consider a given Helmholtz potential $\Psi$ such that the Piola-Kirchhoff
stress tensor is given by $\tensorsym
P\coloneqq\frac{\partial\Psi}{\partial\tensorsym F}$. Then, the linear
momentum conservation is given by

$-\operatorname{\text{div}}\tensorsym P(\tensorsym F)=\vectorsym
g\qquad\text{in $\Omega_{0}$},$ (1)

for any given bulk force $\vectorsym g$, and it represents the balance of
forces. Despite the popularity of this formulation, typically forces are known
in the deformed configuration (denoted by $\vectorsym g_{e}$) and not in the
reference configuration. Note that performing the pull-back of the force
$\vectorsym g_{e}$ to the reference configuration yields the relationship
$\vectorsym g=J\vectorsym g_{e}$, which stems from the change of variables
formula $dx=JdX$. Boundary conditions can be of Dirichlet $\vectorsym
d=\bar{\vectorsym{d}}$, Neumann $\tensorsym P\vectorsym N=\bar{\vectorsym t}$,
and Robin $\tensorsym P\vectorsym N=-\alpha\vectorsym d$ type in non-
overlapping portions of the boundary.

#### Mass conservation.

For the same potential $\Psi$ considered before, we define the pressure $p$ to
be given by $p\coloneqq\frac{\partial\Psi}{\partial\Phi}$. This fact is
deduced from the thermodynamics of poroelastic media [15]. The mass
conservation of the fluid phase is given by

$\frac{d\Phi}{dt}+\operatorname{\text{div}}\vectorsym
U=1/{\rho_{f}}\Theta\qquad\text{in $\Omega_{0}$},$ (2)

where $\Theta$ is a given source/sink term, and $\vectorsym U$ is the relative
velocity of the fluid with respect to the solid phase. If the source/sink term
is given in the Eulerian configuration, then $\Theta=J\theta$ for a given
$\theta$. We adopt a typical choice used in porous media to model the relative
fluid velocity $\vectorsym U$, using Darcy’s law. This law states that the
fluid follows the negative pressure gradient, and is given by

$\vectorsym U\coloneqq-\tensorsym K\operatorname{\nabla}p,$

where the tensor $\tensorsym K$ is obtained as the pullback of a given
Eulerian permeability tensor $\tensorsym k$, i.e. $\tensorsym K=J\tensorsym
F^{-1}\tensorsym k\tensorsym F^{-T}$. Boundary conditions can be of Dirichlet
$\Phi=\bar{\Phi}$, Neumann $\tensorsym K\operatorname{\nabla}p\cdot\vectorsym
N=0$, and Robin type $\tensorsym K\operatorname{\nabla}p\cdot\vectorsym
N=-\alpha\Phi$. Typically, the porous media is considered to be isolated,
meaning that homogeneous Neumann conditions are considered throughout the
boundary. The initial condition is given by a known Lagrangian porosity such
that

$\Phi(0)=\Phi_{0}.$

In case the known initial porosity value is given in Eulerian form, then we
consider an initial Jacobian $J_{0}=1$ such that

$\Phi_{0}=\phi_{0}.$

The lack of an initial Jacobian in our formulation stems from the choice of
using a quasi-static elasticity. Despite this limitation, we maintain this
choice as we found this approach more numerically robust in our preliminary
tests due to oscillations arising from the inertia.

We make two remarks on the chosen law of mass conservation. The first one is
that we are implicitly assuming that the fluid is incompressible as
$\rho_{f}=\rho_{f0}$. This can be justified by adequate thermodynamic
modeling, and is common in the porous media community [14] if the pore content
is given by an incompressible fluid. The second one is that we are only
imposing mass conservation of the fluid phase and not that of the solid phase.
This is the case because mass conservation of the solid phase establishes a
relationship between the solid density in spatial configuration and
deformation, which we do not require in virtue of the quasi-static
approximation and the lack of external body forces. Even in the case of
inertial deformation, the computation of an Eulerian solid density is
typically circumvented by solving the problem in the reference configuration.
This would of course not be true for the Eulerian one.

#### Solid incompressibility.

This condition states that the solid phase is incompressible, and is given in
an arbitrary subdomain $\omega_{0}\subset\Omega_{0}$ by [15]:

$\int_{\omega_{0}}\left(1-\phi_{0}\right)\,dX=\int_{\omega}\left(1-\phi\right)\,dx,$
(3)

for some initial porosity $\phi_{0}$ with $\omega=\vectorsym x(\omega_{0})$.
Equation (3) can be rewritten through a localization argument as

$\Phi_{s}=\phi_{s,0}\qquad\text{in $\Omega_{0}$},$ (4)

where $\Phi_{s}$ and $\phi_{s,0}$ are the porosities of the solid phase in the
Lagrangian and Eulerian form respectively. This yields $J-\Phi=1-\phi_{0}$ and
can be included in the potential as

$\tilde{\Psi}\coloneqq\Psi+\lambda\left(J-\Phi-\left(1-\phi_{0}\right)\right),$

where $\lambda$ is the Lagrange multiplier associated to the constraint. As a
result, the modified first Piola–Kirchhoff stress tensor is given by

$\widetilde{\tensorsym P}=\frac{\partial\tilde{\Psi}}{\partial\tensorsym
F}=\tensorsym P+\lambda J\tensorsym F^{-T},$

and the pressure becomes

$\widetilde{p}=\frac{\partial\tilde{\Psi}}{\partial\Phi}=p-\lambda.$

This fundamental hypothesis establishes the feedback of the fluid phase on the
solid one. To obtain such a feedback without using incompressibility, an
additional energy term must be considered in $\Psi$ that is driven by the
solid porosity and can be interpreted as the solid phase quasi-
incompressibility. See [11] for more details on this topic.

#### Complete forward model.

Putting everything together and assuming homogeneous boundary conditions
(either Dirichlet or Neumann), we obtain that the weak formulation of (1),
(2), and (3) is given by: Find a displacement $\vectorsym d$ in $V^{d}$, a
Lagrangian porosity $\Phi$ in $V^{\Phi}$, and a Lagrange multiplier $\lambda$
in $V^{\lambda}$ such that

$\displaystyle\int_{\Omega_{0}}\tensorsym P(\tensorsym
F,\Phi):\operatorname{\nabla}\vectorsym d^{*}\,dX+\int_{\Omega_{0}}\lambda
J\tensorsym F^{-T}:\operatorname{\nabla}\vectorsym d^{*}\,dX$
$\displaystyle=\int_{\Omega_{0}}\vectorsym g\cdot\vectorsym d^{*}\,dX$
$\displaystyle\forall\vectorsym d^{*}\in V^{d},$ (5)
$\displaystyle\int_{\Omega_{0}}\frac{d\Phi}{dt}\Phi^{*}\,dX+\int_{\Omega_{0}}\tensorsym
K(\tensorsym F)\operatorname{\nabla}\widetilde{p}(\tensorsym
F,\Phi)\cdot\operatorname{\nabla}\Phi^{*}\,dX$
$\displaystyle=\int_{\Omega_{0}}\rho_{f}^{-1}\Theta\Phi^{*}\,dX$
$\displaystyle\forall\Phi^{*}\in V^{\Phi},$
$\displaystyle\int_{\Omega_{0}}\left(J-\Phi\right)\lambda^{*}\,dX$
$\displaystyle=\int_{\Omega_{0}}\left(1-\phi_{0}\right)\lambda^{*}\,dX$
$\displaystyle\forall\lambda^{*}\in V^{\lambda},$

where the first Piola-Kirchhoff tensor $\tensorsym P$ and the pressure $p$ are
computed from the Helmhotz potential $\Psi$ as $\tensorsym
P=\frac{\partial\Psi}{\partial\tensorsym F}$ and
$p=\frac{\partial\Psi}{\partial\Phi}$.

#### Note: Hidden second order derivatives.

We note that problem (2) has second order order derivatives coming from the
pressure gradient, which is given by

$\operatorname{\nabla}p(\tensorsym F,\Phi)=\frac{\partial
p}{\partial\tensorsym F}:\operatorname{\nabla}\tensorsym F+\frac{\partial
p}{\partial\Phi}\operatorname{\nabla}\Phi.$

The conflicting term is $\operatorname{\nabla}\tensorsym F$ because $H^{1}$
conforming finite element schemes do not have a well defined second
derivative, and in any case the displacement is typically an $\vectorsym
H^{1}$ function, not $\vectorsym H^{2}$. This difficulty, which we refer to as
_primal inconsistency_ , will be addressed in Section 4 by formulating mixed
models through physically significant auxiliary variables.

###  2.2 Eulerian field equations

To define the Eulerian formulation, we follow the inverse displacement model
from [17]. Consider the inverse displacement $\hat{\vectorsym d}$ such that
$\vectorsym X(\vectorsym x)=\vectorsym x+\hat{\vectorsym d}(\vectorsym x)$,
which yields a strain tensor $\tensorsym f=\tensorsym
I+\operatorname{\nabla}_{\vectorsym x}\hat{\vectorsym d}$ and Jacobian
$j=\det\tensorsym f$ that are related to the forward quantities through
$\tensorsym f=\tensorsym F^{-1}$ and $j=1/J$. Interestingly, the elasticity
problem in Eulerian configuration yields a problem that completely determines
the inverse displacement $\hat{\vectorsym d}$. In what follows, we pose the
current configuration equations of linear momentum (1), mass conservation (2),
and solid incompressibility (3). We highlight that in the current
configuration, the natural porosity to consider as a variable is the Eulerian
one, i.e. $\phi$.

#### Equilibrium equation.

Changing variables in the first Piola–Kirchhoff stress tensor we obtain that
$\tensorsym P(\tensorsym F)=\tensorsym P(\tensorsym f^{-1})$, so that the
push-forward of the linear momentum equation is given by

$-\operatorname{\text{div}}_{\vectorsym x}\left(j\tensorsym P(\tensorsym
f^{-1})\tensorsym f^{-1}\right)=\vectorsym g_{e}\quad\text{in $\Omega$},$ (6)

for a given volumetric load $\vectorsym g_{e}$, where
$\tensorsym\sigma(\tensorsym f)\coloneqq j\tensorsym P(\tensorsym
f^{-1})\tensorsym f^{-1}$ is known as the Cauchy stress tensor. Boundary
conditions are given as in the Lagrangian case, with the difference that
Neumann and Robin boundary conditions are written using spatial quantities as
$\tensorsym\sigma\vectorsym n=\tensorsym t$.

#### Mass conservation.

We recast mass conservation in terms of its natural Eulerian variable
$\phi=\Phi/J=j\Phi$, and denote the Eulerian fluid velocity $\vectorsym
u=J\tensorsym F^{-T}\vectorsym U$ to write mass conservation as

$\frac{d\phi}{dt}+\operatorname{\text{div}}_{\vectorsym x}\vectorsym
u=\frac{1}{\rho_{f}}\theta\quad\text{in $\Omega$},$ (7)

for a given source/sink term $\theta$. See [15] for a derivation of this
equation. As before, the fluid velocity is computed from Darcy’s law

$\vectorsym u\coloneqq-\tensorsym k\operatorname{\nabla}_{\vectorsym x}p.$

Boundary conditions are given as in the Lagrangian case. The choice of the
initial conditions will be discussed in Section 2.4 as it is one of the main
contributions of this work.

#### Solid incompressibility.

Pushing forward (3) to $\Omega$, we obtain the following Eulerian form of
solid incompressibility:

$1-\phi=j\left(1-\phi_{0}\right)\qquad\text{in $\Omega$}.$ (8)

This is incorporated in our formulation by modifying the Helmholtz potential
as

$\tilde{\Psi}=\Psi+\lambda\left(1-\phi-j\left(1-\phi_{0}\right)\right),$

where $\lambda$ is the Lagrange multiplier associated to the constraint. With
it, the Cauchy stress tensor becomes

$\widetilde{\tensorsym\sigma}=\tensorsym\sigma+\lambda\tensorsym I,$

and the pressure remains the same. Note that the stress generated by the
Lagrange multiplier is significantly simpler in the Eulerian formulation.

###  2.3 Constitutive modeling

Even though our formulation is valid independently of the mechanical material
model used, we will consider a material related to cardiac modeling, since
cardiac applications are our main interest, and are also numerically
challenging. In virtue of this, the Helmholtz potential $\Psi$, as in [11],
will be separated additively into a mechanical contribution $\Psi_{M}$ and
porous contributions $\Psi_{P}$ as

$\Psi(\tensorsym F,\Phi)=\Psi_{M}(\tensorsym F)+\Psi_{P}(\Phi),$

where the term $\Psi_{P}$ acts as a barrier function such that $\Phi>0$.
Similarly to [14], the inequality $\phi<1$ (or equivalently $\Phi<J$) is
naturally satisfied using solid incompressibility. Indeed, using equation (8)

$\phi=1-j\phi_{s,0}\qquad\quad\text{(equivalently $\Phi=J-\phi_{s,0}$)},$

and assuming an arbitrarily small, but positive, constant $\epsilon$ such that
$\phi=1-\epsilon$, we obtain

$\epsilon=j\phi_{s,0}=J^{-1}\phi_{s,0}.$

Thus, $\epsilon\to 0$ if $j\to 0$, i.e. when $J\to\infty$, which is
nonphysical and already characterized by the potential $\Psi_{M}$. In this
work, we consider an Usyk energy of ventricular muscle [30] for $\Psi_{M}$,
given by

$\Psi_{M}(\tensorsym F)=C\left(\exp(Q(\overline{\tensorsym
E}))-1\right)+\frac{B}{2}(J-1)\log J,$

where

$\displaystyle Q(\overline{\tensorsym
E})=b_{\mathrm{ff}}\overline{E}_{\mathrm{ff}}^{2}+b_{\mathrm{ss}}\overline{E}_{\mathrm{ss}}^{2}+b_{\mathrm{nn}}\overline{E}_{\mathrm{nn}}^{2}+b_{\mathrm{fs}}\left(\overline{E}_{\mathrm{fs}}^{2}+\overline{E}_{\mathrm{sf}}^{2}\right)+b_{\mathrm{fn}}\left(\overline{E}_{\mathrm{fn}}^{2}+\overline{E}_{\mathrm{nf}}^{2}\right)+b_{\mathrm{sn}}\left(\overline{E}_{\mathrm{sn}}^{2}+\overline{E}_{\mathrm{ns}}^{2}\right),$
$\displaystyle\overline{E}_{\text{ab}}=\overline{\tensorsym
E}\vectorsym{a}\cdot\vectorsym{b},\qquad
a,b\in\\{f,s,n\\},\qquad\overline{E}=\frac{1}{2}\left(\overline{\tensorsym
F}^{T}\overline{\tensorsym F}-\tensorsym I\right),\qquad\overline{\tensorsym
F}=J^{1/3}\tensorsym F,$

where $(\vectorsym f,\vectorsym s,\vectorsym n)$ is a locally orthonormal
basis that represents the orientation of the cardiac fibers,
$C=880\,\texttt{Pa}$ and $B=5\times 10^{4}\,\texttt{Pa}$. For the porous
energy $\Psi_{P}$ we consider the arterial Bruinsma et al. energy [5],
upscaled in [13] to be given by

$\Psi_{P}(\Phi)=\frac{q_{1}}{q_{3}}\exp\left(q_{3}\Phi\right)+q_{2}\Phi\left(\log\left(q_{3}\Phi\right)-1\right)$

where the values of the parameters $q_{i}$ for arteries, capillaries, and
veins have been computed in [5]. Still, as observed in [11], the energy
$\Psi_{P}$ must be normalized such that it is minimized at the reference
values of $\Phi$, i.e. the pressure at rest yields the reference pressures.
Therefore, we make the dependence of the energy explicit on the reference
porosities as

$\widetilde{\Psi}_{P}(\Phi,\Phi_{0})\coloneqq\Psi_{P}(\Phi)-\frac{\partial\Psi_{P}}{\partial\Phi}(\tensorsym
I,\Phi_{0})\Phi.$ (9)

By considering also the solid incompressibility constraint, we obtain the
following expressions for the first Piola–Kirchhoff stress tensor and the
pressure:

$\displaystyle\tensorsym P(\tensorsym F)$
$\displaystyle=\frac{\partial\Psi_{M}}{\partial\tensorsym F}(\tensorsym
F)+\lambda J\tensorsym F^{-T},$ $\displaystyle p(\Phi,\Phi_{0})$
$\displaystyle=\frac{\partial\Psi_{P}}{\partial\Phi}(\Phi)-\frac{\partial\Psi_{P}}{\partial\Phi}(\Phi_{0})-\lambda.$

We make the following very important remarks.

######  Remark 1.

The first Piola–Kirchhoff tensor and the pressure seem decoupled from one
another. Indeed, it is the Lagrange multiplier that yields the interaction
between fluid and solid phases, and it can be interpreted as the pressure
existing at the interface between both phases.

######  Remark 2.

The pressure contribution from the reference porosity $\Phi_{0}$ has a flipped
sign with respect to the pressure obtained from the energy. This flipped sign,
together with the negative sign obtained from the time-reversal strategy
proposed, are what yield that the nonlinear porous media equation behaves as a
heat equation instead of a reversed heat equation.

######  Remark 3.

The pressure is considered such that it is gives a reference porosity of 0 at
$\Phi_{0}$ for simplicity. It could of course be considered that the reference
porosity is non-zero, i.e.

$p(\Phi,\Phi_{0})=\frac{\partial\Psi_{P}}{\partial\Phi}(\Phi)-\frac{\partial\Psi_{P}}{\partial\Phi}(\Phi_{0})+p_{\texttt{ref}}-\lambda,$

for some reference pressure $p_{\texttt{ref}}$.

###  2.4 Formulation of the reference configuration model

In this section we extend the elastic reference configuration model to
poroelasticity. The extension is based on two main points: (i) the notion of
mechanical equilibrium can be understood as a steady state solution in the
context of porous media, and (ii) the resulting steady state equation can
present extreme nonlinearities that are very challenging to solve
computationally. This motivates considering another related time dependent
problem instead of the steady state problem directly. We now provide the
details of each point separately.

#### Equilibrium.

The reference configuration problem departs from the assumption that we know
not only the forces acting of a system, but also its _resulting_ solution. The
time dependent nature of system (7) makes it unclear at exactly which time we
observe the system, so we assume that the solution given is _stationary_ ,
i.e. it is the solution $(\hat{\vectorsym d},\phi_{0})$ of (6) and the
following stationary mass conservation problem:

$-\operatorname{\text{div}}_{\vectorsym x}\tensorsym
k\operatorname{\nabla}_{\vectorsym x}p(\tensorsym
f^{-1},j\phi,\phi_{0})=\frac{1}{\rho_{f}}\theta\qquad\text{on}\quad\Omega,$
(10)

as seen from (7). We highlight that in this model the given porosity is
$\phi$, and the unknown is the reference one $\phi_{0}$. The way in which
$\theta$ is considered here is problem dependent. For example, if it varies in
time in the forward model, then its initial value could be considered. If it
is instead an autonomous function of the variables (as we do in the numerical
tests section), i.e. $\theta=\theta(\tensorsym f^{-1},j\phi,\phi_{0})$, then
it can be left as it is.

#### Time dependence.

As seen from Section 2.3, the pressure function can be severely nonlinear, and
indeed all of our preliminary attempts at solving problem (10) resulted in
diverging nonlinear iterations. For this reason, and to allow for more
flexible continuation methods, we provide a way to reconsider time dependence
in (10) in a way that yields the reference configuration solution. Equation
(7) can be rewritten as

$\frac{\partial(\phi-\phi_{0})}{\partial
t}-\operatorname{\text{div}}_{\vectorsym x}\tensorsym
k\operatorname{\nabla}_{\vectorsym x}p=\frac{1}{\rho_{f}}\theta,$

as $\phi_{0}$ is constant in time and the quantity $\phi-\phi_{0}$ can be
regarded as the variation of porosity. Inverting the variables, we consider a
given (time constant) porosity $\phi=\overline{\phi}$ and write the reference
mass conservation problem as

$\displaystyle-\frac{\partial\phi_{0}}{\partial
t}-\operatorname{\text{div}}_{\vectorsym x}\tensorsym
k\operatorname{\nabla}_{\vectorsym x}p(\tensorsym
f^{-1},j\overline{\phi},\phi_{0})$ $\displaystyle=\frac{1}{\rho_{f}}\theta$
$\displaystyle\text{on}\quad\Omega,$ (11) $\displaystyle\phi_{0}(0)$
$\displaystyle=\overline{\phi}$ $\displaystyle\text{on}\quad\Omega.$

Naturally, a steady-state solution of problem (11) provides a solution to
(10). If we think of (11) as a heat equation for a moment, the reference mass
conservation equation behaves as a reversed heat equation, which is severely
ill-posed. This is not the case here because the normalization factor
$\frac{\partial\Psi_{P}}{\partial\Phi}(\Phi_{0})=\frac{\partial\Psi_{P}}{\partial\Phi}(\phi_{0})$
shown in Section 2.3 has a negative sign.

#### Complete reference configuration model.

Putting everything together, we obtain the weak formulation of (6), (11) and
(8) assuming homogeneous Dirichlet or Neumann boundary conditions as: Find an
inverse displacement $\hat{\vectorsym d}$ in $V^{d}$, a spatial reference
porosity $\phi_{0}$ in $V^{\phi}$, and a Lagrange multiplier $\lambda$ in
$V^{\lambda}$ such that

$\displaystyle\int_{\Omega}\left(j\tensorsym P(\tensorsym f^{-1})\tensorsym
f^{-T}\right):\operatorname{\nabla}_{\vectorsym x}\vectorsym
d^{*}\,dx+\int_{\Omega}\lambda\operatorname{\text{div}}_{\vectorsym
x}\vectorsym d^{*}\,dx$ $\displaystyle=\int_{\Omega}\vectorsym
g_{e}\cdot\vectorsym d^{*}\,dx$ $\displaystyle\forall\vectorsym d^{*}\in
V^{d},$ (12)
$\displaystyle-\int_{\Omega}\frac{d\phi_{0}}{dt}\phi^{*}\,dx+\int_{\Omega}\tensorsym
k\operatorname{\nabla}_{\vectorsym x}\widetilde{p}(\tensorsym
f^{-1},j\overline{\phi},\phi_{0}):\operatorname{\nabla}_{\vectorsym
x}\phi^{*}\,dx$ $\displaystyle=\int_{\Omega}\rho_{f}^{-1}\theta\phi^{*}\,dx$
$\displaystyle\forall\phi^{*}\in V^{\phi},$
$\displaystyle\int_{\Omega}j\left(1-\phi_{0}\right)\lambda^{*}\,dx$
$\displaystyle=\int_{\Omega}\left(1-\overline{\phi}\right)\lambda^{*}\,dx$
$\displaystyle\forall\lambda^{*}\in V^{\lambda},$

where $\tensorsym P=\frac{\partial\Psi}{\partial\tensorsym F}$,
$p=\frac{\partial\Psi}{\partial\Phi}$, and $\phi_{0}(0)=\overline{\phi}$.

##  3 Numerical solution strategy

In this section we provide numerical details for solving the reference
configuration problem. All computations were performed using the Firedrake
library [29]. We briefly provide the discretization choices used to solve both
(5) and (12):

#### Space discretization.

We use the inf-sup stable finite element space
$\mathbb{P}_{2}\times\mathbb{P}_{1}\times\mathbb{P}_{1}$ [7], where higher
order elements are used for the displacement. Firedrake ships Gauss-Jacobi
quadrature rules by default, which we have set to a degree of 6 to handle the
nonlinearities adequately. Anything lower than that exhibits numerical
artifacts in our experience.

#### Time discretization.

Given the quasi-static nature of most equations, we only need to provide a
time-discretization for the porous media equations. These equations are
parabolic (in virtue of the convex potential) [32], so we use a simple
backward Euler scheme [20], and thus treat all variables implicitly.

#### Iterative solvers.

The nonlinear system arising from the implicit discretization is solved using
Newton’s method [34], where the linear system is solved at each time step
using the direct solver available in the MUMPS library [1], which includes an
MPI-parallel implementation. Better iterative solvers for nonlinear
poroelasticity based on Krylov Subspace methods can be found in [6], but their
extension to the fully nonlinear case we considered has not yet been
developed.

#### Continuation.

The model considered is highly nonlinear, and thus a naive use of Newton’s
method yields diverging iterations. To avoid this, we consider a continuation
strategy for the source term, such that for a given source term $\theta^{*}$
(possibly time dependent) we consider the linear ramp given by
$\theta(t)=\min\\{t/t_{\text{ramp}},1\\}\theta^{*}(t)$, so that the source
term is fully active at $t=t_{\text{ramp}}$.

#### Convergence criterion.

To compute the reference solutions as stationary states of equation (12), we
solve the equations until the normalized residual of the steady state mass
conservation equation (10) is below a given threshold. We thus set

$R_{i}(t)\coloneqq\int_{\Omega}\tensorsym k\operatorname{\nabla}_{\vectorsym
x}p\cdot\operatorname{\nabla}_{\vectorsym
x}\phi_{i}^{*}-\rho_{f}^{-1}\theta(t)\phi_{i}^{*}\,dx,$

where $\phi_{i}^{*}$ stands for the $i$-th test function in the discrete
space, and use as a convergence criterion the following:

$\left|\vectorsym R(t)\right|_{\ell^{2}}\leq\texttt{tol}R_{0},$

where $[\vectorsym R]_{i}=R_{i}$,
$|(\cdot)|_{\ell^{2}}\coloneqq\sqrt{\sum_{i}(\cdot)^{2}_{i}}$, and
$R_{0}\coloneqq|R(t_{\text{ramp}})|_{\ell^{2}}$ is a normalization factor. We
assume a tolerance of $\texttt{tol}=10^{-6}$ was used unless otherwise stated.

###  3.1 Accelerated stationary state computation

In this section we show that the computation of a stationary state can be
accelerated through a fixed-point acceleration technique known as Anderson
acceleration [3]. The method changes a fixed point iteration
$x^{k}=g(x^{k-1})$ with a memory dependent algorithm given by

$x^{k}=\sum_{i\in\\{1,...,m\\}}\alpha_{i}g(x^{k-i}),$

where the $\alpha_{i}$ weights are optimal for a given norm, and $m$ is known
as the depth parameter. More details on how to compute the weights, as well as
the theoretical foundations of this technique, can be found in [35]. We
highlight that this method is equivalent to the GMRES method in the linear
case, and to the multisecant method in the nonlinear case.

To apply Anderson acceleration to the problem of finding a stationary state,
we formulate a fixed-point operator in two steps. First, we denote with
$\mathbb{T}:V^{i}\to V^{d}\times V^{\phi}\times V^{\lambda}$ the operator
given by

$(\hat{\vectorsym d}^{k},\phi^{k},\lambda^{k})=\mathbb{T}(\phi^{k-1}),$

where $(\cdot)^{k}$ is the discrete solution at instant $t^{k}=k\Delta t$.
Then we denote with $\Pi_{\phi}$ the projection

$\phi=\Pi_{\phi}(\hat{\vectorsym d},\phi,\lambda).$

With these operators, we note that the problem of finding a stationary state
to equation (12) can be recast as finding a fixed-point of the operator
$\mathbb{S}:V^{\phi}\to V^{\phi}$ given by

$\mathbb{S}=\Pi_{\phi}\circ\mathbb{T}.$

This operator is amenable to acceleration, and we compare the number of
iterations it takes to compute the stationary state with and without Anderson
acceleration in Section 5. Naturally, all computations performed here can be
trivially extended to the forward problem.

##  4 Consistent mixed models

In this section we propose two mixed models associated with nonlinear
poroelasticity that can circumvent the primal inconsistency problem, detailed
in Section 2.1. The first one is the commonly used mixed formulation given by
considering the fluid velocity as a variable $\vectorsym
u\coloneqq-\tensorsym\kappa\operatorname{\nabla}p$, which is useful also for
extended models in which the Brinkman effect is included in the fluid inertia.
The second one is a simpler formulation in which the pressure becomes a
variable $\mu\coloneqq p(\tensorsym F,\Phi)$. Given that the main difficulty
is the treatment of the mass conservation equation, throughout this section we
focus on the following equation:

$\frac{d\phi}{dt}-\operatorname{\text{div}}\left(\tensorsym
k\operatorname{\nabla}p(\tensorsym F,\phi)\right)=\theta\qquad\text{in
$\Omega$},$ (13)

for given functions $\tensorsym F$ and $\theta$. We highlight that this
problem is independent of the previous models, so we consider $\Omega$ to be a
general domain in this section, and we avoid unnecessary notation regarding
the differentiation variables. The boundary conditions for each formulation
will be different, thus, to highlight how to convert from one formulation to
the other one, we consider mixed boundary conditions on Dirichlet and Neumann
boundaries, $\Gamma_{D}$ and $\Gamma_{N}$ respectively, with

$\phi=\phi_{D}\quad\text{on $\Gamma_{D}$},\quad\text{and}\quad\tensorsym
k\operatorname{\nabla}p\cdot\vectorsym n=\bar{p}\quad\text{on $\Gamma_{N}$},$

where it holds that
$\overline{\partial\Omega}=\overline{\Gamma_{D}}\cup\overline{\Gamma_{N}}$.

###  4.1 Mixed formulation in pressure

Starting from equation (13), we define the variable $\mu\coloneqq p(\tensorsym
F,\phi)$ so that we can rewrite the problem as the following system of
equations:

$\displaystyle\frac{d\phi}{dt}-\operatorname{\text{div}}\tensorsym
k\operatorname{\nabla}\mu$ $\displaystyle=\theta$ $\displaystyle\text{in
$\Omega$},$ (14) $\displaystyle\mu-p(\tensorsym F,\phi)$ $\displaystyle=0$
$\displaystyle\text{in $\Omega$}.$

From the integration by parts of the second order operator against a test
function $\mu^{*}$:

$-\int_{\Omega}\operatorname{\text{div}}\left(\tensorsym
k\operatorname{\nabla}\mu\right)\mu^{*}\,dx=\int_{\Omega}\tensorsym
k\operatorname{\nabla}\mu\cdot\operatorname{\nabla}\mu^{*}\,dx-\int_{\partial\Omega}\left(\tensorsym
k\operatorname{\nabla}\mu\cdot\vectorsym n\right)\mu^{*}\,ds$

we can see that the Dirichlet boundary conditions for this model are given by

$\mu=\mu_{D}\quad\text{on $\Gamma_{D}$},$

which in virtue of the definition of $\mu$, correspond to imposing pressure
values on the boundary instead of porosity values as in the primal formulation
(13). The Neumann boundary condition is instead given by

$\tensorsym k\operatorname{\nabla}\mu\cdot\vectorsym n=\bar{p}\quad\text{on
$\Gamma_{N}$},$

which in virtue of the definition of $\mu$, correspond to the same Neumann
boundary conditions used in (13), thus justifying the use of the same function
$\bar{p}$. We highlight that these boundary conditions are better from a
physical point of view, as typically in real applications the boundary values
correspond to pressure instead of porosity.

The weak formulation of (14) is given by finding $\phi$ in $V^{\phi}$, and
$\mu$ in $V^{\mu}$ such that

$\displaystyle\int_{\Omega}\frac{d\phi}{dt}\mu^{*}\,dx+\int_{\Omega}\tensorsym
k\operatorname{\nabla}\mu\cdot\operatorname{\nabla}\mu^{*}\,dx$
$\displaystyle=\int_{\Omega}\theta\mu^{*}\,dx$ $\displaystyle\forall\mu^{*}\in
V^{\mu},$ (15)
$\displaystyle\int_{\Omega}\mu\phi^{*}\,dx-\int_{\Omega}p(\tensorsym
F,\phi)\phi^{*}\,dx$ $\displaystyle=0$ $\displaystyle\forall\phi^{*}\in
V^{\phi}.$

This formulation does not present a gradient on the nonlinear function $p$,
thus is avoids the primal inconsistency problem.

######  Remark 4.

The choice of switching the test functions is purely theoretical. To see this,
consider the simplest case where $p(\tensorsym F,\phi)=\phi$ and an implicit
time discretization with time-step $\Delta t=1$. If we use the test functions
$(\mu^{*},\phi^{*})=(\mu,-\phi)$ in (15) and add both equations, we obtain

$\int_{\Omega}|\operatorname{\nabla}\mu|^{2}\,dx+\int_{\Omega}\phi^{2}\,dx,$

which hints on better theoretical properties for this formulation, as well as
not requiring an inf-sup condition among the spaces $V^{\phi}$ and $V^{\mu}$.
This in particular allows for a robust lowest order discretization.

###  4.2 Mixed formulation in fluid velocity

Proceeding as before, we define the function $\vectorsym u\coloneqq-\tensorsym
k\operatorname{\nabla}p$, which converts (13) into the following system of
equations:

$\displaystyle\frac{d\phi}{dt}+\operatorname{\text{div}}\vectorsym u$
$\displaystyle=\theta$ $\displaystyle\text{in $\Omega$},$ (16)
$\displaystyle\tensorsym k^{-1}\vectorsym u+\operatorname{\nabla}p$
$\displaystyle=0$ $\displaystyle\text{in $\Omega$}.$

Again, we deduce the boundary conditions from the integration by parts of the
$\operatorname{\nabla}p$ term against a test function $\vectorsym u^{*}$:

$\int_{\Omega}\vectorsym
u^{*}\cdot\operatorname{\nabla}p\,dx=-\int_{\Omega}p\operatorname{\text{div}}\vectorsym
u^{*}\,dx+\int_{\partial\Omega}p\vectorsym u^{*}\cdot\vectorsym n\,dx.$

In this case, the Dirichlet boundary conditions correspond to

$\vectorsym u\cdot\vectorsym n=-\bar{p}\quad\text{on $\Gamma_{D}$},$

and the Neumann ones correspond to

$p=p_{N}\quad\text{ on $\Gamma_{N}$}.$

Interestingly, the boundary conditions of this formulation correspond to the
ones from (14), but with exchanged roles (i.e. Dirichlet conditions become
Neumann conditions and vice-versa). This phenomenon is typical in mixed
methods [16], but the nonlinearity of $p$ breaks that relationship with
respect to the original primal formulation (13).

Its weak formulation is given by finding $\vectorsym u$ in $V^{u}$ and $\phi$
in $V^{\phi}$ such that

$\displaystyle\int_{\Omega}\frac{d\phi}{dt}\phi^{*}\,dx+\int_{\Omega}\phi^{*}\operatorname{\text{div}}\vectorsym
u\,dx$ $\displaystyle=\int_{\Omega}\theta\phi^{*}\,dx$
$\displaystyle\forall\phi^{*}\in V^{\phi},$ (17)
$\displaystyle\int_{\Omega}\tensorsym k^{-1}\vectorsym u\cdot\vectorsym
u^{*}\,dx-\int_{\Omega}p(\tensorsym F,\phi)\operatorname{\text{div}}\vectorsym
u^{*}\,dx$ $\displaystyle=0$ $\displaystyle\forall\vectorsym u^{*}\in V^{u}.$

Similarly to (15), this formulation does not present a gradient on the
nonlinear function $p$, and thus it also avoids the primal inconsistency
error.

######  Remark 5.

The mixed $\vectorsym u$ formulation presents two main advantages: (i) it
yields the physically significant variable $\vectorsym u$ with optimal
convergence rates, and (ii) it allows for generalized Darcy models such as the
ones considering a Brinkman term, or also fluid inertia. Still, this comes
with the disadvantage of being a more difficult problem. Indeed, considering
the test functions $(\vectorsym u^{*},\phi^{*})=(\vectorsym u,\phi)$ as
before, we obtain

$\int_{\Omega}|\vectorsym u|^{2}\,dx+\int_{\Omega}|\phi|^{2}\,dx,$

but the presence of the divergence operator shows we require $\vectorsym u$ to
be in $H(\operatorname{\text{div}};\Omega)$. This can be overcome using the
theory of saddle point problems, which also shows that inf-sup stable finite
element spaces need to be chosen for the $V^{u}\times V^{\phi}$ pair.

######  Remark 6.

The mixed models follow the same logic used to define the fixed-point operator
in Section 3.1, meaning that they are also amenable to acceleration.

We show the main features of all three models in Table 1. We emphasize the
difference in terms of the variables and the boundary conditions, where again
we highlight the following: (i) The mixed models present more physically
adequate boundary conditions, as the quantity that is typically measured is
the pressure, not the porosity, (ii) the role of the boundary conditions
between the mixed models is exchanged, and (iii) the mixed models do not
require second order derivatives of the displacement.

Model | Variables | Dirichlet BC | Neumann BC | Consistent
---|---|---|---|---
Primal | $\phi$ | $\phi=\bar{\phi}$ | $\operatorname{\nabla}p\cdot\vectorsym n=\bar{\vectorsym u}$ | ✗
Mixed-$p$ | $\phi,\mu$ | $\mu=\bar{p}$ | $\operatorname{\nabla}p\cdot\vectorsym n=\bar{\vectorsym u}$ | ✓
Mixed-$u$ | $\phi,\vectorsym u$ | $\vectorsym u\cdot\vectorsym n=\bar{\vectorsym u}$ | $p=\bar{p}$ | ✓
Table 1: Comparison of variables and boundary conditions required for the
primal, pressure mixed and velocity mixed models.

##  5 Numerical tests

In this section we provide numerical tests that support our theory. The models
have been implemented in Firedrake [29] following the numerical solution
strategy detailed in Section 3, and the visualizations are performed with
Paraview [4]. The tests considered, and their scope, are as follows:

Test 1:

We compute the reference configuration of a 2D geometry with a physically
motivated nonlinear source term. The scope of this test is to show that the
loaded reference configuration yields the given current configuration, and
also recovers the given porosity.

Test 2:

We repeat Test 1 for a 3D geometry.

Test 3:

We verify the impact of using Anderson acceleration to compute the desired
stationary state in both 2D and 3D tests, so we perform acceleration for
varying levels of memory and assess its performance by measuring the number of
iterations required for convergence.

Test 4:

We verify that the mixed models are also able to yield a satisfactory solution
to the reference configuration problem by comparing the solution obtained by
all 3 models (primal, pressure mixed, and velocity mixed).

Test 5:

We compute the reference configuration of a realistic geometry using
physically accurate parameters.

In both square and slab tests, we use the canonical basis as fiber directions.
Instead, for the left ventricle simulation we use physically accurate fiber
directions, computed numerically as in [8]. More details on this are provided
in the corresponding test.

###  5.1 Test 1: Validation in 2D

In this test, we validate the performance of the reference configuration
formulation. The way to do this is by loading the problem with a pressure
driven source term that models the inflow of fluid according to a pressure
difference:

$\theta=-\beta(p(\tensorsym F,\Phi)-p_{a}),$

where $p_{a}=10^{4}\,\texttt{Pa}$ and $\beta=10^{-4}\,\frac{1}{\texttt{s
Pa}}$. For the porous media potential we used the arterial pressure model
described in Section 2.3 where $q_{1}=1.333\,\texttt{Pa}$,
$q_{2}=550\,\texttt{Pa}$, and $q_{3}=10$, with a given porosity of
$\overline{\phi}=0.1$. Finally, we considered an isotropic permeability given
by the scalar $k=2\cdot\,10^{-7}\texttt{m}^{2}\texttt{(s\,Pa)}^{-1}$. We
considered a square geometry $\Omega$ of side length $1\,\texttt{cm}$, a time
step of $\Delta t=0.01\,s$, and a ramp time given by $t_{\text{ramp}}=0.1\,s$.
As boundary conditions, we used Dirichlet conditions that allowed for sliding,
i.e. $\vectorsym d_{x}=0$ on $\\{x=0\\}$ and $\vectorsym d_{y}=0$ on
$\\{y=0\\}$, and homogeneous Neumann conditions elsewhere.

In Figure 1, we show the deformed solution in two snapshots: In the left
figure, we show the computed reference configuration with respect to the
initial geometry. In the right one, we show the solution obtained by solving
the forward problem on the computed reference configuration, where the
accuracy of this method can be appreciated. Additionally, we computed the
average porosity $\frac{1}{|\Omega|}\int_{\Omega}\phi\,dx$ and displayed its
evolution on the current configuration in Figure 2, where we plotted the
porosity for both reference configuration and forward problems. Here, it is
possible to observe that the reference configuration deflates and loses a
significant amount of fluid, given by the 84% of the initial value as it
arrives at an average of 0.014. The recovered porosity presents a small error
of roughly a 1.1% with respect to the correct value, which we attribute to
approximation error.

(a)

(b)

Figure 1: (a) the computed reference configuration with respect to the initial
geometry (black contour). (b) the forward problem solution as computed
starting from the reference configuration computed in (a).
$0$$100$$200$$300$$400$$500$$600$$700$$2\cdot 10^{-2}$$4\cdot 10^{-2}$$6\cdot
10^{-2}$$8\cdot 10^{-2}$$0.1$IterationsRef. conf.ForwardGiven
($\overline{\phi}$) Figure 2: Evolution of the average spatial porosity
$\frac{1}{|\Omega|}\int_{\Omega}\phi\,dx$ along the iterations of both the
reference configuration and forward problems. The given value
$\overline{\phi}$ is shown for comparison.

###  5.2 Test 2: Validation in 3D

Similar to test 1, we compute the reference configuration for a 3D slab
geometry given by
$\Omega=(0,5\,\texttt{cm})\times(0,1\,\texttt{cm})\times(0,1\,\texttt{cm})$.
The only difference with respect to the 2D test is the relative tolerance
$\texttt{tol}=10^{-5}$. In Figure 3, we show the deformed solution in two
snapshots as before: reference configuration and forward problems in the left
and right figures respectively. As in the previous case, the forward solution
presents an excellent match with respect to the initial geometry. The
reference porosities are displayed in Figure 4, where the results obtained are
analogous to the ones from Test 1.

(a)

(b)

Figure 3: (a) the computed reference configuration with respect to the initial
geometry (black contour). (b) the forward problem solution as computed
starting from the reference configuration in (a).
$-10$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$$110$$120$$130$$2\cdot
10^{-2}$$4\cdot 10^{-2}$$6\cdot 10^{-2}$$8\cdot 10^{-2}$$0.1$IterationsRef.
conf.ForwardGiven ($\overline{\phi}$) (a) Phase 1 Figure 4: Evolution of the
average spatial porosity $\frac{1}{|\Omega|}\int_{\Omega}\phi\,dx$ along the
iterations of both the reference configuration and forward problems in both
phases for the brick 3D geometry. The given value $\overline{\phi}$ is shown
for comparison.

###  5.3 Test 3: Acceleration test

In this section, we test the impact of using Anderson acceleration to reduce
the iterations required for convergence in both Tests 1 and 2. The number of
iterations required for both the reference configuration and forward problems
is shown in Table 2, with respect to different levels of acceleration. For the
sake of comparison, we used a tolerance of $\texttt{tol}=10^{-5}$ for both 2D
and 3D tests.

Acceleration is not used immediately, but only when the ramping has finished
in order to accelerate with respect to the desired dynamics, which happens
after 10 timesteps. As seen from Table 2, we make the following observations:
(i) both formulations require roughly 100 iterations for convergence without
acceleration, (ii) acceleration works even using only 1 previous vector, and
(iii) all accelerated formulations take roughly less than 80% the number of
iterations it takes to solve the problem, where using only one vector yields
the lowest iteration count. We conclude that acceleration is an effective
strategy to reduce the overall computational cost of computing the reference
configuration as a stationary state of the nonlinear porous media equations.

| AA(0) | AA(1) | AA(2) | AA(5)
---|---|---|---|---
Ref. conf. 2D | 98 | 15 | 17 | 20
Forward 2D | 122 | 14 | 15 | 18
Ref. conf. 3D | 98 | 15 | 17 | 20
Forward 3D | 122 | 14 | 15 | 18
Table 2: Test 3, number of time steps required to converge to a stationary
state using Anderson Acceleration for both forward and reference configuration
models.

###  5.4 Test 4: Validation of mixed models

This test computes the reference configuration of the 2D square geometry as in
Test 1, using the mixed models presented in Section 4. The evolution of the
average porosity in all three models is shown in Figure 5, where all models
yield the same solution. Indeed, the average reference porosity obtained by
the primal, mixed in pressure, and mixed in velocity models is given by
$0.140767$, $0.140767$, and $0.140597$ respectively. We chose to prioritize
stability in order to have an accurate comparison of the methods, so we used
second order elements for the fluid velocity in the mixed velocity model. In
addition, both primal and mixed pressure models took the same amount of
iterations to converge (680), while the mixed velocity model took 842
iterations to converge.

$-50$$0$$50$$100$$150$$200$$250$$300$$350$$400$$450$$500$$550$$600$$650$$700$$750$$2\cdot
10^{-2}$$4\cdot 10^{-2}$$6\cdot 10^{-2}$$8\cdot 10^{-2}$$0.1$IterationsRef.
conf.ForwardGiven ($\overline{\phi}$) (a) Primal

$-50$$0$$50$$100$$150$$200$$250$$300$$350$$400$$450$$500$$550$$600$$650$$700$$750$$2\cdot
10^{-2}$$4\cdot 10^{-2}$$6\cdot 10^{-2}$$8\cdot 10^{-2}$$0.1$IterationsRef.
conf.ForwardGiven ($\overline{\phi}$) (b) Mixed $p$

$0$$100$$200$$300$$400$$500$$600$$700$$800$$900$$2\cdot 10^{-2}$$4\cdot
10^{-2}$$6\cdot 10^{-2}$$8\cdot 10^{-2}$$0.1$IterationsRef. conf.ForwardGiven
($\overline{\phi}$) (c) Mixed $u$

Figure 5: Evolution of the average spatial porosity
$\frac{1}{|\Omega|}\int_{\Omega}\phi\,dx$ along the iterations of both the
reference configuration and forward problems. We show the evolution obtained
using the (a) primal formulation, (b) the mixed pressure formulation and (c)
the mixed velocity formulation.

###  5.5 Test 5: Realistic perfusion model

In this test, we compute the reference configuration associated to a perfusion
model of an idealized left ventricular geometry. For this, the modifications
to the model considered are: (i) boundary conditions for the mechanics are
divided into an endocardial pressure of $1.5\,\texttt{kPa}$ (physiological
end-diastolic pressure) and friction due to the pericardium on the epicardium
as in [26], (ii) the source term $\theta$ considers both an arterial source
and a venous sink of $5\,\texttt{kPa}$ and $1\,\texttt{kPa}$ respectively,
(iii) the passive material is the anisotropic one from [31], with the fiber,
sheet and normal directions computed using an accurate $H^{1}$ conforming
approximation proposed in [8], and (iv) we replace the solid incompressibility
with a quasi-incompressible law of the same magnitude of the geometry quasi-
incompressibility as in [11]. The geometry, the computed reference
configuration and the forward solution are displayed in Figure 6, together
with their corresponding reference porosities in the caption. The solution
displays a very good match with respect to the given geometry, and the average
stationary porosity yields a small error of a 2% (0.098 instead of 0.1) due to
numerical accuracy. Additionally, the average stationary porosity obtained
with the primal model is 0.094, showing that not using a consistent mixed
formulation can yield in this simple test an additional 4% error. We expect
this number to increase in more complex scenarios.

We note that this test was significantly more difficult than the others, as it
would frequently result in diverging nonlinear iterations. To circumvent this
difficulty, instead of increasing our ramp slowly and then computing a
stationary state, we computed the stationary state for each ramp level before
increasing it. A natural improvement for this was Anderson acceleration, so in
Figure 7 we show the number of fixed point iterations accumulated at each ramp
level for 1, 2, and 5 previous vectors, and in Table 3 we show the total
number of iterations (rightmost value at Figure 7). We note that iterations
are reduced by up to roughly a 63% and an 80% of the ones obtained without
acceleration ($AA(0)$) for the forward and reference configuration problems
respectively, and interestingly, the optimal number of previous vectors
depends on the formulation. Still, considering also the previous tests, we
believe that a depth of 2 is a good general choice.

(a) $\overline{\phi}=0.1$

(b) $\phi_{h,0}=0.075$

(c) $\phi_{h}=0.098$

Figure 6: The average Eulerian steady steate porosities in the caption for the
(a) given geometry, (b) computed reference configuration with the given
geometry in a traslucent blue, and (c) the deformed geometry from the given
loads with the given geometry in a traslucent blue. The ’h’ subscript shows
that the function is a FEM approximation.

$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$0$$50$$100$$150$$200$Time$AA(0)$$AA(1)$$AA(2)$$AA(5)$
(a) Ref. conf.

$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$20$$40$Time$AA(0)$$AA(1)$$AA(2)$$AA(5)$
(b) Forward

Figure 7: Number of accumulated time iterations incurred to solve the (a) reference configuration and (b) forward problems at each ramp step. In (b), the $AA(2)$ and $AA(5)$ lines overlap. | $AA(0)$ | $AA(1)$ | $AA(2)$ | $AA(5)$
---|---|---|---|---
Ref. conf. | 187 | 39 | 37 | 56
Forward | 52 | 22 | 19 | 19
Table 3: Total number of fixed point iterations incurred when solving both
reference configuration and forward models in the left ventricle test.

##  6 Discussion

This work proposes a novel model for computing the reference configuration of
a fully nonlinear poroelastic material, where the main ingredient is the
identification of the given porous solution as the steady state configuration
of a porous media problem with unknown initial conditions. The proposed model
is consistent with fundamental thermodynamics principles, which renders it
applicable for all applications where nonlinear poroelasticity is relevant. We
have tested the model using parameters and laws coming from cardiac
applications, which involve the combination of logarithmic and exponential
nonlinearities that render the model computationally challenging.

One interesting discovery in this work is the _primal inconsistency_ of the
primal formulation, which can be avoided using mixed formulations. We proposed
two approaches, one that yields a natural energy norm and another one that
yields a saddle point problem. The saddle point formulation is much more
challenging, but it can be useful when considering more general fluid models.
Our main interest is to consider Darcy-Brinkman models for cardiac ablation,
which has been shown to be relevant in applications [33]. However solved, the
computation of a stationary state is an expensive problem, which results in
simulations with hundreds of iterations. To circumvent this difficulty, we
have shown that Anderson acceleration yields an extremely successful strategy
to reduce iteration counts, with reductions of up to an 80%.

Our future work mainly focuses on two topics: the formulation of robust
reference configuration models for multi-phase porous media, which we observed
to be much more difficult to solve, and the use of the proposed models to
accurately model the interplay between blood flow and deformation during
cardiac ablation as in [27], by extending the large deformations model
proposed in [22].

##  7 Acknowledgments

NB has been funded by the ANID Postdoctoral 3230326, by CMM BASAL FB2100005
and partially supported by Johannes Kepler University. AP and LGG were
partially supported by the State of Upper Austria.

## References

  * ADLK [00] PR Amestoy, IS Duff, J-Y L’Excellent, and J Koster.  Mumps: a general purpose distributed memory sparse solver.  In International Workshop on Applied Parallel Computing, pages 121–130. Springer, 2000.
  * AKL [10] D Algranati, GS Kassab, and Y Lanir.  Mechanisms of myocardium-coronary vessel interaction.  American Journal of Physiology-Heart and Circulatory Physiology, 298(3):H861–H873, 2010.
  * And [65] DG Anderson.  Iterative procedures for nonlinear integral equations.  Journal of the ACM (JACM), 12(4):547–560, 1965.
  * Aya [15] U Ayachit.  The paraview guide: a parallel visualization application.  Kitware, Inc., 2015.
  * BADS [88] P Bruinsma, T Arts, J Dankelman, and JAE Spaan.  Model of the coronary circulation based on pressure dependence of coronary resistance and compliance.  Basic research in cardiology, 83:510–524, 1988.
  * BDOVMRB [22] NA Barnafi, LM De Oliveira Vilaca, MC Milinkovitch, and R Ruiz-Baier.  Coupling chemotaxis and growth poromechanics for the modelling of feather primordia patterning.  Mathematics, 10(21):4096, 2022.
  * BGVL+ [22] NA Barnafi, B Gómez-Vargas, WJ Lourenço, RF Reis, BM Rocha, M Lobosco, R Ruiz-Baier, and RW dos Santos.  Finite element methods for large-strain poroelasticity/chemotaxis models simulating the formation of myocardial oedema.  Journal of Scientific Computing, 92(3):1–40, 2022.
  * BO [23] NA Barnafi and A Osses.  Modeling of cardiac fibers as oriented liquid crystals.  arXiv preprint arXiv:2305.01199, 2023.
  * BRR [23] NA Barnafi, F Regazzoni, and D Riccobelli.  Reconstructing relaxed configurations in elastic bodies: Mathematical formulation and numerical methods for cardiac modeling.  arXiv preprint arXiv:2312.11477, 2023.
  * BV [87] J Bear and A Verruijt.  Theory and applications of transport in porous media.  Springer, 1987.
  * BWDGD+ [22] NA Barnafi Wittwer, S Di Gregorio, L Dede’, P Zunino, C Vergara, and AM Quarteroni.  A multiscale poromechanics model integrating myocardial perfusion and the epicardial coronary vessels.  SIAM Journal on Applied Mathematics, 82(4):1167–1193, 2022.
  * CGSMVC [10] D Chapelle, J-F Gerbeau, J Sainte-Marie, and IE Vignon-Clementel.  A poroelastic model valid in large strains with applications to perfusion in cardiac modeling.  Computational Mechanics, 46:91–101, 2010.
  * CLM+ [12] AN Cookson, J Lee, C Michler, R Chabiniok, E Hyde, DA Nordsletten, M Sinclair, M Siebes, and NP Smith.  A novel porous mechanical framework for modelling the interaction between coronary perfusion and myocardial mechanics.  Journal of biomechanics, 45(5):850–855, 2012.
  * CM [14] D Chapelle and P Moireau.  General coupling of porous flows and hyperelastic formulations - from thermodynamics principles to energy balance and compatible time schemes.  European Journal of Mechanics, B/Fluids, 46:82–96, 2014.
  * Cou [04] O Coussy.  Poromechanics.  John Wiley & Sons, 2004.
  * Gat [14] GN Gatica.  A simple introduction to the mixed finite element method.  Theory and Applications. Springer Briefs in Mathematics. Springer, London, 2014.
  * GM [96] S Govindjee and PA Mihalic.  Computational methods for inverse finite elastostatics.  Computer Methods in Applied Mechanics and Engineering, 136(1-2):47–57, 1996.
  * GREW [09] MW Gee, CH Reeps, HH Eckstein, and WA Wall.  Prestressing in finite deformation abdominal aortic aneurysm simulation.  Journal of biomechanics, 42(11):1732–1739, 2009.
  * HAVCR [92] JM Huyghe, T Arts, DH Van Campen, and RS Reneman.  Porous medium finite element model of the beating left ventricle.  American Journal of Physiology-Heart and Circulatory Physiology, 262(4):H1256–H1267, 1992.
  * HNW [93] E Hairer, SP Nørsett, and G Wanner.  Solving ordinary differential equations i. nonstiff problems, 1993.
  * LA [16] A LoMauro and A Aliverti.  Physiology masterclass: Extremes of age: newborn and infancy.  Breathe, 12(1):65–68, 2016.
  * MGGG [22] L Molinari, L Gerardo-Giorda, and A Gizzi.  A transversely isotropic thermo-hyperelastic constitutive model of myocardial tissue with a three-state cell death dynamics for cardiac radiofrequency ablation.  Journal of the Mechanics and Physics of Solids, 161:104810, 2022.
  * PDBSR [18] M Peirlinck, M De Beule, P Segers, and N Rebelo.  A modular inverse elastostatics approach to resolve the pressure-induced stress state for in vivo imaging based cardiovascular modeling.  Journal of the mechanical behavior of biomedical materials, 85:124–133, 2018.
  * PGA+ [23] Argyrios Petras, Matthias AF Gsell, Christoph M Augustin, Jairo Rodriguez-Padilla, Alexander Jung, Marina Strocchi, Frits W Prinzen, Steven A Niederer, Gernot Plank, and Edward J Vigmond.  Mechanoelectric effects in healthy cardiac function and under left bundle branch block pathology.  Computers in Biology and Medicine, 156:106696, 2023.
  * PGC [22] C Patte, M Genet, and D Chapelle.  A quasi-static poromechanical model of the lungs.  Biomechanics and Modeling in Mechanobiology, 21(2):527–551, 2022.
  * PHW+ [19] MR Pfaller, JM Hörmann, M Weigl, A Nagler, R Chabiniok, C Bertoglio, and WA Wall.  The importance of the pericardium for cardiac biomechanics: from physiology to computational modeling.  Biomechanics and modeling in mechanobiology, 18:503–529, 2019.
  * PLG+ [19] A Petras, M Leoni, JM Guerra, J Jansson, and L Gerardo-Giorda.  A computational model of open-irrigated radiofrequency catheter ablation accounting for mechanical properties of the cardiac tissue.  International Journal for Numerical Methods in Biomedical Engineering, 35(11):e3232, 2019.
  * RGH [17] MK Rausch, M Genet, and JD Humphrey.  An augmented iterative method for identifying a stress-free reference configuration in image-based biomechanical modeling.  Journal of biomechanics, 58:227–231, 2017.
  * RHM+ [16] F Rathgeber, DA Ham, L Mitchell, M Lange, F Luporini, ATT McRae, G-T Bercea, GR Markall, and PHJ Kelly.  Firedrake: automating the finite element method by composing abstractions.  ACM Transactions on Mathematical Software (TOMS), 43(3):1–27, 2016.
  * ULM [02] TP Usyk, IJ LeGrice, and AD McCulloch.  Computational model of three-dimensional cardiac electromechanics.  Computing and Visualization in Science, 4(4):249–257, 2002.
  * UMM [00] TP Usyk, R Mazhari, and AD McCulloch.  Effect of laminar orthotropic myofiber architecture on regional stress and strain in the canine left ventricle.  Journal of Elasticity, 61:143–164, 2000.
  * Váz [07] JL Vázquez.  The porous medium equation: mathematical theory.  Oxford University Press, 2007.
  * WDR [23] P Wongchadakul, AK Datta, and P Rattanadecho.  Tissue poromechanical deformation effects on steam pop likelihood in 3-d radiofrequency cardiac ablation.  Journal of Biological Engineering, 17(1):52, 2023.
  * WN [99] S Wright and J Nocedal.  Numerical optimization.  Springer Science, 35(67-68):7, 1999.
  * WN [11] HF Walker and P Ni.  Anderson acceleration for fixed-point iterations.  SIAM Journal on Numerical Analysis, 49(4):1715–1735, 2011.