remarkRemark hypothesisHypothesis claimClaim Multiple 3D solutions in topology
optimizationI. P. A. Papadopoulos and P. E. Farrell

# Preconditioners for computing multiple solutions in three-dimensional fluid
topology optimization ††thanks: Submitted DATE. This manuscript was prepared
whilst the first author was at the University of Oxford and was supported by
the EPSRC Centre for Doctoral Training in Partial Differential Equations:
Analysis and Applications [grant number EP/L015811/1] and The MathWorks, Inc.
Revisions were supported by the EPSRC grant “Spectral element methods for
fractional differential equations, with applications in applied analysis and
medical imaging” [grant number EP/T022132/1] and the Leverhulme Trust Research
Project “Constructive approximation theory on and inside algebraic curves and
surfaces” [grant number RPG-2019-144]. The second author is supported by the
Engineering and Physical Sciences Research Council [grant numbers EP/R029423/1
and EP/W026163/1].

Ioannis P. A. Papadopoulos Department of Mathematics, Imperial College London,
London, UK
<EMAIL_ADDRESS>Patrick E. Farrell Mathematical
Institute, University of Oxford, Oxford, UK<EMAIL_ADDRESS>
###### Abstract

Topology optimization problems generally support multiple local minima, and
real-world applications are typically three-dimensional. In previous work [I.
P. A. Papadopoulos, P. E. Farrell, and T. M. Surowiec, _Computing multiple
solutions of topology optimization problems_ , SIAM Journal on Scientific
Computing, (2021)], the authors developed the deflated barrier method, an
algorithm that can systematically compute multiple solutions of topology
optimization problems. In this work we develop preconditioners for the linear
systems arising in the application of this method to Stokes flow, making it
practical for use in three dimensions. In particular, we develop a nested
block preconditioning approach which reduces the linear systems to solving two
symmetric positive-definite matrices and an augmented momentum block. An
augmented Lagrangian term is used to control the innermost Schur complement
and we apply a geometric multigrid method with a kernel-capturing relaxation
method for the augmented momentum block. We present multiple solutions in
three-dimensional examples computed using the proposed iterative solver.

###### keywords:

topology optimization, multiple solutions, deflated barrier method,
preconditioning, multigrid

65F08, 65K10, 65N55, 35Q35, 49M15, 90C26

##  1 Introduction

Topology optimization has become a prominent tool in engineering design. The
objective is to find the optimal distribution of a continuum within a given
domain that minimizes a problem-specific cost functional, without prior
knowledge of the shape or topology of the solution [11].

Topology optimization problems are generally nonconvex and can support
multiple local minima. A major challenge in topology optimization is to
identify multiple local minimizers, so that the best one (in performance,
manufacturability, or aesthetics) may be chosen. Continuation in model
parameters is often used to promote convergence to “better” local minima.
However, Stolpe and Svanberg [65] show that such methods can be ineffective
even on the simplest of topology optimization models.

In recent work, Papadopoulos et al. [51] developed an algorithm, called the
deflated barrier method, that can systematically discover _multiple_
stationary points of topology optimization problems formulated using a density
approach. The deflated barrier method combines the techniques of barrier
methods [33, 34, 36, 60, 61, 67, 70], primal-dual active set solvers [12], and
deflation [26, 27]. In one example [51, Fig. 1], the authors found 42
stationary points of a fluid topology optimization problem constrained by the
Navier–Stokes equations in a rectangular domain with five small holes removed.
However, an application of the deflated barrier method requires the solution
of many linear systems similar to those solved in an all-at-once (simultaneous
analysis and design, SAND) approach. In [51], a sparse LU factorization was
used to solve these linear systems, which is not generally affordable for
three-dimensional problems.

In this work, we develop preconditioners for the linear systems arising in the
deflated barrier method when applied to the model proposed by Borrvall and
Petersson [15]. The goal of this model is to minimize the power dissipation of
a fluid in Stokes flow, subject to a volume constraint restricting the
proportion of the domain that the fluid can occupy. The preconditioners we
develop make it feasible to identify multiple solutions to the
Borrvall–Petersson problem in three dimensions with the deflated barrier
method. A number of strategies have been proposed for the solution of three-
dimensional topology optimization of fluid flow in other contexts [2, 3, 4,
10, 20, 21, 25, 42, 53, 58, 74, 76].

We choose a discontinuous Galerkin $H(\mathrm{div})$-conforming finite element
discretization for the velocity and pressure pair [16, 17] due to its
pressure-robustness [43] and its crisp characterization of the kernel of the
divergence operator. We will show that block preconditioning can reduce the
linear systems arising in the deflated barrier method to ones that resemble
the systems arising in the discretization of the Stokes–Brinkman equations
[25]. Then, we apply modern block preconditioning, pioneered by Wathen and
coworkers [22, 64, 69] to solve the inner linear systems. The inner linear
systems are themselves solved with block preconditioning, using an augmented
Lagrangian term to control the innermost Schur complement [35]. Finally, we
develop a geometric multigrid method for the augmented momentum block with a
vertex-star patch relaxation that captures the kernel of the augmented
Lagrangian term [29, 31, 40, 62]. The multigrid scheme also requires a
characterization of the active set on all levels of the mesh hierarchy, which
we discuss.

We introduce the Borrvall–Petersson model in Section 2 and the deflated
barrier method in Section 3. In Section 4 we detail the discontinuous Galerkin
finite element discretization and derive the linear systems that arise in the
subproblems of the deflated barrier method. In Section 5 we develop the
preconditioning strategies for these linear systems. In Section 6 we
investigate a number of two- and three-dimensional examples. All the examples
support multiple solutions which are discovered by the deflated barrier
method, and solved with our proposed iterative methods.

##  2 Topology optimization of Stokes flow

Given a volume constraint on a fluid in a fixed Lipschitz domain
$\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, the Borrvall–Petersson model
attempts to minimize the energy lost by the flow due to viscous dissipation,
whilst maximizing the flow velocities at the applied body force. More
precisely, the objective is to find $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ that
minimizes

(BP) $\displaystyle
J(\boldsymbol{u},\rho)\coloneqq\frac{1}{2}\int_{\Omega}\left(\alpha(\rho)|\boldsymbol{u}|^{2}+\nu|\nabla\boldsymbol{u}|^{2}-2\boldsymbol{f}\cdot\boldsymbol{u}\right)\text{d}x,$

where $\boldsymbol{u}:\Omega\to\mathbb{R}^{d}$ denotes the velocity of the
fluid, $\rho:\Omega\to\mathbb{R}$ is the material distribution of the fluid
and

$\displaystyle H^{1}_{\boldsymbol{g}}(\Omega)^{d}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in
H^{1}(\Omega)^{d}:\boldsymbol{v}|_{\partial\Omega}=\boldsymbol{g}\\},$
$\displaystyle H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in
H^{1}_{\boldsymbol{g}}(\Omega)^{d}:\mathrm{div}(\boldsymbol{v})=0\;\text{a.e.\
in}\;\Omega\\},$ $\displaystyle C_{\gamma}$
$\displaystyle\coloneqq\left\\{\eta\in L^{\infty}(\Omega):0\leq\eta\leq
1\;\text{a.e.},\;\;\int_{\Omega}\eta\;\text{d}x\leq\gamma|\Omega|\right\\}.$

Here, $H^{1}(\Omega)$ denotes the Sobolev space $W^{1,2}(\Omega)$, equipped
with the inner product

(1) $\displaystyle(u,v)_{H^{1}(\Omega)}=\int_{\Omega}uv+\nabla u\cdot\nabla
v\,\mathrm{d}x,$

which induces the norm
$\|v\|_{H^{1}(\Omega)}\coloneqq(v,v)_{H^{1}(\Omega)}^{1/2}$.
$L^{\infty}(\Omega)$ denotes the vector space of essentially bounded
measurable functions equipped with the essential supremum norm, i.e.
$\|\eta\|_{L^{\infty}(\Omega)}\coloneqq\inf\\{c\geq 0:|\eta(x)|\leq
c\;\text{for a.e.}\;x\in\Omega\\}$. Furthermore, $\boldsymbol{f}\in
L^{2}(\Omega)^{d}$ is a body force, $\nu>0$ is the (constant) viscosity, and
$\gamma\in(0,1)$ is the volume fraction. The restriction,
$|_{\partial\Omega}$, is to be understood in the boundary trace sense [24, Ch.
5.5]. Moreover, the boundary data $\boldsymbol{g}\in
H^{1/2}(\partial\Omega)^{d}$ and $\boldsymbol{g}=\boldsymbol{0}$ on
$\Gamma\subset\partial\Omega$, with $\mathcal{H}^{d-1}(\Gamma)>0$, i.e.
$\Gamma$ has nonzero Hausdorff measure on the boundary. Borrvall and Petersson
introduced the inverse permeability term, $\alpha$, which models the influence
of the material distribution on the flow. For values of $\rho$ close to one,
$\alpha(\rho)$ is small, permitting fluid flow; for values of $\rho$ close to
zero, $\alpha(\rho)$ is very large, restricting fluid flow. The function
$\alpha$ satisfies the following properties:

  1. (A1)

$\alpha:[0,1]\to[\underline{\alpha},\overline{\alpha}]$ with
$0\leq\underline{\alpha}<\overline{\alpha}<\infty$;

  2. (A2)

$\alpha$ is monotonically decreasing;

  3. (A3)

$\alpha(0)=\overline{\alpha}$ and $\alpha(1)=\underline{\alpha}$;

  4. (A4)

$\alpha$ is twice continuously differentiable,

  5. (A5)

$\alpha$ is strongly convex,

generating an operator also denoted $\alpha:C_{\gamma}\to
L^{\infty}(\Omega;[\underline{\alpha},\overline{\alpha}])$. Typically, in the
literature, $\alpha$ takes the form [15, 25]

(2)
$\displaystyle\alpha(\rho)=\bar{\alpha}\left(1-\frac{\rho(q+1)}{\rho+q}\right),$

where $q>0$ is a penalty parameter, so that
$\lim_{q\to\infty}\alpha(\rho)=\bar{\alpha}(1-\rho)$. Borrvall and Petersson
[15, Sec. 3.2] remark that as $q\to\infty$ the material distribution tends to
a 0-1 solution.

The following existence theorem is due to Borrvall and Petersson [15, Th.
3.1].

######  Theorem 2.1.

Suppose that $\Omega\subset\mathbb{R}^{d}$ is a Lipschitz domain, with $d=2$
or $d=3$ and $\alpha$ satisfies properties (A1)–(A5). Then, there exists a
pair $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)\times C_{\gamma}$ that minimizes
$J$ (as defined in Eq. BP).

##  3 The deflated barrier method

In this section, we summarize how the deflated barrier method finds multiple
solutions of fluid topology optimization problems.

###  3.1 Forming the barrier functional

The first step in the algorithm is to form a barrier-like objective from the
original objective functional. This is both to aid convergence and to provide
an opportunity for deflation to discover multiple local minima (described in
Section 3.2). Consider the minimization problem: find the minimizers,
$(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$, of

(BF) $\displaystyle
J_{\mu}^{\epsilon_{\text{log}}}(\boldsymbol{u},\rho)\coloneqq
J(\boldsymbol{u},\rho)-\mu\int_{\Omega}\log(\rho+\epsilon_{\text{log}})+\log(1+\epsilon_{\text{log}}-\rho)\mathrm{d}x,$

where $\mu\geq 0$ is the barrier parameter and $0<\epsilon_{\text{log}}\ll 1$.
We note that the box constraints imposed by the barrier-like terms are never
active as $0\leq\rho\leq 1$ a.e. in $\Omega$. A similar result to one shown by
Evgrafov [25, Sec. 4] guarantees the existence of a minimizer to Eq. BF [50,
Prop. 4.2]. As with the original optimization problem, Eq. BP, the minimizer
is not necessarily unique.

The deflated barrier method targets the first-order optimality conditions of
this problem, which we derive next. Define the following forms to be used in
stating the first-order optimality conditions:

(3) $\displaystyle c(\rho,\eta;\boldsymbol{u},\lambda)$
$\displaystyle\coloneqq\int_{\Omega}\left[\frac{1}{2}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}-\frac{\mu}{\rho+\epsilon_{\text{log}}}+\frac{\mu}{1+\epsilon_{\text{log}}-\rho}+\lambda\right]\eta\;\mathrm{d}x,$
(4) $\displaystyle a(\boldsymbol{u},\boldsymbol{v};\rho)$
$\displaystyle\coloneqq\int_{\Omega}\alpha(\rho)\boldsymbol{u}\cdot\boldsymbol{v}+\nu\nabla\boldsymbol{u}:\nabla\boldsymbol{v}\;\mathrm{d}x,\quad
l(\boldsymbol{v})\coloneqq\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}\;\mathrm{d}x,$
(5) $\displaystyle b(\boldsymbol{u},p)$
$\displaystyle\coloneqq-\int_{\Omega}p\nabla\cdot\boldsymbol{u}\;\mathrm{d}x,\quad
d(\lambda,\zeta;\rho)\coloneqq-\zeta\int_{\Omega}(\gamma-\rho)\mathrm{d}x.$

The proof of the following proposition, concerning the first-order optimality
conditions of Eq. BF, follows from the result in [50, Prop. 2.5].

######  Proposition 3.1.

Fix $\mu\geq 0$, $0<\epsilon_{\mathrm{log}}\ll 1$ and suppose that
$(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ is a strict
local or global minimizer of Eq. BF. Then, there exist unique Lagrange
multipliers $p\in L^{2}_{0}(\Omega)$ and $\lambda\in\mathbb{R}$, such that,
for all $(\eta,\boldsymbol{v},q,\zeta)\in C_{[0,1]}\times
H^{1}_{0}(\Omega)^{d}\times L^{2}_{0}(\Omega)\times\mathbb{R}$, the following
necessary first-order optimality conditions are satisfied:

(6) $\displaystyle c(\rho,\eta-\rho;\boldsymbol{u},\lambda)$
$\displaystyle\geq 0,$ (7) $\displaystyle
a(\boldsymbol{u},\boldsymbol{v};\rho)+b(\boldsymbol{v},p)$
$\displaystyle=l(\boldsymbol{v}),$ (8) $\displaystyle b(\boldsymbol{u},q)$
$\displaystyle=0,$ (9) $\displaystyle d(\lambda,\zeta;\rho)$
$\displaystyle=0.$

Here $C_{[0,1]}\coloneqq\cup_{\gamma\in[0,1]}C_{\gamma}$, i.e. we relax the
volume constraint on the variation in the material distribution.

Under suitable conditions, Eq. 6–Eq. 9 can be reformulated as a
complementarity problem. After a suitable discretization, throughout the
deflated barrier method, the nonlinear complementarity problems are solved
using the Benson–Munson (BM) algorithm [12], a Newton-like algorithm that can
enforce the true box constraints on $\rho$. A thorough description of this
strategy is given in [51, Sec. 3.1] and Section 4.2. We emphasize that the
barrier-like terms are _not_ used to enforce box constraints, which are
handled by BM, and the enlarged box constraints
$[-\epsilon_{\text{log}},1+\epsilon_{\text{log}}]$ are never active for any
$\mu\geq 0$. The barrier-like terms are included solely to aid nonlinear
convergence via continuation, and to find multiple solutions via deflation.

###  3.2 Deflation

Deflation is a mechanism to systemically discover multiple solutions of a
nonlinear system with a Newton-like algorithm [26, 27]. Let $Z$ and $Y$ be
Banach spaces. Consider the nonlinear system $F(z)=0,$ $F:Z\to Y$ that has the
solutions $z_{1},\dots,z_{l}$. Given a suitable initial guess, a Newton-like
algorithm might converge to the solution $z_{1}$. Now, deflation modifies the
nonlinear system in order to remove the known solution whilst retaining
undiscovered solutions. This is done via a deflation operator
$\mathcal{M}(z;z_{1}):Y\to Y$ that ensures that $\mathcal{M}(z;z_{1})F(z)$ has
a root if the original problem $F(z)=0$ has an unknown solution, and ensures
that a Newton-like method applied to the newly deflated system does not
converge to $z=z_{1}$. In this work, the following shifted deflation operator
[26] is used

(10)
$\displaystyle\mathcal{M}(\boldsymbol{z};\boldsymbol{z}_{1})=\left(\frac{1}{\|\rho-\rho_{1}\|^{2}_{L^{2}(\Omega)}}+1\right)\mathcal{I},$

where $\boldsymbol{z}=(\rho,\boldsymbol{u},p,\lambda)$ and $\mathcal{I}:Y\to
Y$ is the identity operator. Deflation can be very efficiently implemented;
the BM update of the deflated system can be expressed as a (nonlinear) scaling
of the BM update of the _undeflated_ system evaluated at the same iterate [26,
Sec. 3]. This implementation detail is crucial for this work as it implies
that preconditioning strategies for the undeflated BM systems can be
immediately applied to the deflated systems. A discussion of how the nonlinear
scaling is cheaply computed is given in [51, Sec. 3.2].

###  3.3 Prediction, continuation, and deflation

The deflated barrier method [51, Alg. 3.1] is an iterative procedure
consisting of three steps. The algorithm is initialized with barrier parameter
$\mu=\mu_{0}$ and the idea is to follow branches (a generalization of the
central path) of solutions as $\mu\to 0$. At iteration $k$ of the deflated
barrier method, the barrier functional subproblem at $\mu=\mu_{k}$ is
constructed. In the _prediction step_ , a cheap initial guess is computed for
the barrier functional subproblem from the solution on the same branch at
$\mu=\mu_{k-1}$. In the _continuation step_ , the BM solver is used to find a
solution to the first-order optimality conditions Eq. 6–Eq. 9 at
$\mu=\mu_{k}$, using the initial guess computed in the prediction step, whilst
deflating away all other known solutions at $\mu=\mu_{k}$. In the _deflation
step_ , we search for new branches of solutions at $\mu=\mu_{k}$, using
solutions discovered at $\mu=\mu_{k-1}$ as initial guesses. The deflation step
is terminated if a pre-determined number of branches have already been found
or the BM solver fails to converge (by reaching a specified number of
iterations without converging). The algorithm terminates when the first-order
optimality conditions have been solved (for multiple solutions) at a barrier
parameter value of zero.

##  4 Discretization and linearization

###  4.1 Finite element discretization

The linear systems that arise in the deflated barrier method will be tied to
our choice of finite element discretization. In this work, we choose a
discontinuous Galerkin (DG) $k$th-order
$\mathrm{BDM}_{k}\times\mathrm{DG}_{k-1}$ Brezzi–Douglas–Marini finite element
discretization for the velocity and the pressure [16, 17]. Let
$\mathcal{T}_{h}$ denote a shape-regular triangulation of the domain $\Omega$.
We define the spaces $\boldsymbol{X}_{\mathrm{BDM}_{k}}$ and
$\boldsymbol{X}_{\mathrm{BDM}_{k}^{\boldsymbol{g}}}$ as:

(11) $\displaystyle\boldsymbol{X}_{\mathrm{BDM}_{k}}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in
L^{2}(\Omega)^{d}:\,\nabla\cdot\boldsymbol{v}\in
L^{2}(\Omega),\boldsymbol{v}|_{T}\in\mathcal{P}_{k}(T)^{d}\,\forall\,T\in\mathcal{T}_{h}\\},$
(12) $\displaystyle\boldsymbol{X}_{\mathrm{BDM}_{k}^{\boldsymbol{g}}}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in\boldsymbol{X}_{\mathrm{BDM}_{k}}:(\boldsymbol{v}-\boldsymbol{g})\cdot\boldsymbol{n}=0\;\text{on}\;\partial\Omega\\},$

where $\mathcal{P}_{k}$ denotes the set of polynomials of order $k$. The
degrees of freedom for two dimensions are given in [17, Sec. 2] and for three
dimensions in [16, Sec. 2]. Similarly, $X_{\mathrm{DG}_{k-1}}$ denotes the set
of discontinuous piecewise $(k-1)$-th order polynomials:

(13) $X_{\mathrm{DG}_{k-1}}\coloneqq\\{v\in
L^{2}(\Omega):v|_{T}\in\mathcal{P}_{k-1}(T)\,\forall\,T\in\mathcal{T}_{h}\\}.$

Some finite element methods for fluid flow, such as the Taylor–Hood finite
element pair, do not satisfy the incompressibility constraint
$\nabla\cdot\boldsymbol{u}=0$ pointwise. Failure to satisfy the
incompressibility constraint pointwise has been observed to support
instabilities that result in nonphysical solutions [43, 46]. In the BDM finite
element pair, $\nabla\cdot\boldsymbol{X}_{\mathrm{BDM}_{k}}\subset
X_{\mathrm{DG}_{k-1}}$. Therefore, for any solution $\boldsymbol{u}_{h}$
satisfying the incompressibility constraint Eq. 8, we have that
$\|\nabla\cdot\boldsymbol{u}_{h}\|_{L^{2}(\Omega)}=0$. Hence, the solution is
pointwise divergence-free.

The finite element space for the material distribution is denoted by
$C_{[0,1],h}$ and is defined by

(14) $\displaystyle C_{[0,1],h}\coloneqq\\{\eta_{h}\in
X_{\mathrm{DG}_{0}}:0\leq\eta_{h}\leq 1\\},$

where $X_{\mathrm{DG}_{0}}$ is the set of piecewise constant finite element
functions. As the test functions of $\mathrm{DG}_{0}$ take the value of one in
their respective element and are zero elsewhere, in practice we discretize the
material distribution with $\mathrm{DG}_{0}$ functions and allow the
Benson–Munson strategy described below to handle the box constraints.

Since the discontinuous velocity space is not $H^{1}$-conforming, we use an
interior penalty to penalize jumps across edges and faces. We now describe the
DG discretization of Eq. 6–Eq. 9 as found in [37, Sec. 7.1] and [49, Sec. 3].
The interior penalty only arises in the discretization of Eq. 7.

Let the set $\mathcal{F}_{h}$ denote the set of all facets of the
triangulation $\mathcal{T}_{h}$ and $h_{F}$ represent the diameter of each
facet $F\in\mathcal{F}_{h}$. We split the set of facets into the union
$\mathcal{F}_{h}=\mathcal{F}^{i}_{h}\cup\mathcal{F}^{\partial}_{h}$ where
$\mathcal{F}^{i}_{h}$ is the subset of interior facets and
$\mathcal{F}^{\partial}_{h}$ collects all Dirichlet boundary facets
$F\subset\partial\Omega$. For every facet $F\in\mathcal{F}_{h}$, we assign a
unit norm vector $\boldsymbol{n}_{F}$, where if
$F\in\mathcal{F}^{\partial}_{h}$, then $\boldsymbol{n}_{F}$ is the outer unit
normal vector $\boldsymbol{n}$. If $F\in\mathcal{F}^{i}_{h}$, then
$F=\overline{\partial K^{+}}\cap\overline{\partial K^{-}}$ for two elements
$K^{-},K^{+}\in\mathcal{T}_{h}$ and $\boldsymbol{n}_{F}$ points in an
arbitrary but fixed direction. Let
$\boldsymbol{\phi}\in(X_{\mathrm{DG}_{k}})^{d}$ and
$\boldsymbol{\mathsf{\Phi}}\in(X_{\mathrm{DG}_{k}})^{d\times d}$ be any
piecewise vector- or matrix-valued function, with traces from within the
interior of $K^{\pm}$ denoted by $\boldsymbol{\phi}^{\pm}$ and
$\boldsymbol{\mathsf{\Phi}}^{\pm}$, respectively. We define the jump
$[\\![\cdot]\\!]_{F}$ and the average $\\{\\!\\!\\{\cdot\\}\\!\\!\\}_{F}$
operators across interior facets $F\in\mathcal{F}^{i}_{h}$ by

(15)
$\displaystyle[\\![\boldsymbol{\phi}]\\!]_{F}=\boldsymbol{\phi}^{+}\otimes\boldsymbol{n}_{F}^{+}+\boldsymbol{\phi}^{-}\otimes\boldsymbol{n}_{F}^{-}\quad\text{and}\quad\\{\\!\\!\\{\boldsymbol{\mathsf{\Phi}}\\}\\!\\!\\}_{F}=\frac{1}{2}\left(\boldsymbol{\mathsf{\Phi}}^{+}+\boldsymbol{\mathsf{\Phi}}^{-}\right).$

If $F\in\mathcal{F}^{\partial}_{h}$, we set
$[\\![\boldsymbol{\phi}]\\!]_{F}=\boldsymbol{\phi}\otimes\boldsymbol{n}_{F}$
and
$\\{\\!\\!\\{\boldsymbol{\mathsf{\Phi}}\\}\\!\\!\\}_{F}=\boldsymbol{\mathsf{\Phi}}$.
We note that, for any $F\in\mathcal{F}^{\partial}_{h}$,
$\int_{F}|[\\![\boldsymbol{\phi}]\\!]_{F}|^{2}\,\mathrm{d}s=\int_{F}|\boldsymbol{\phi}|^{2}\,\mathrm{d}s$.
Finally, for a sufficiently large penalization parameter $\sigma>0$, we define
the broken form $a_{h}(\boldsymbol{u},\boldsymbol{v};\rho)$ by

(16)
$\displaystyle\begin{split}a_{h}(\boldsymbol{u},\boldsymbol{v};\rho)&\coloneqq\sum_{K\in\mathcal{T}_{h}}\int_{K}\alpha(\rho)\boldsymbol{u}\cdot\boldsymbol{v}+\nu\nabla\boldsymbol{u}:\nabla\boldsymbol{v}\;\mathrm{d}x+\nu\sum_{F\in\mathcal{F}_{h}}\sigma
h_{F}^{-1}\int_{F}[\\![\boldsymbol{u}]\\!]_{F}:[\\![\boldsymbol{v}]\\!]_{F}\mathrm{d}s\\\
&-\nu\sum_{F\in\mathcal{F}_{h}}\int_{F}\\{\\!\\!\\{\nabla\boldsymbol{u}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{v}]\\!]_{F}\mathrm{d}s-\nu\sum_{F\in\mathcal{F}_{h}}\int_{F}[\\![\boldsymbol{u}]\\!]_{F}:\\{\\!\\!\\{\nabla\boldsymbol{v}\\}\\!\\!\\}_{F}\mathrm{d}s,\end{split}$

and the linear functional $l_{h}$ as

(17) $\displaystyle l_{h}(\boldsymbol{v})$
$\displaystyle\coloneqq\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}\;\mathrm{d}x+\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\sigma
h_{F}^{-1}\int_{F}[\\![\boldsymbol{g}]\\!]_{F}:[\\![\boldsymbol{v}]\\!]_{F}\;\mathrm{d}s-\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}[\\![\boldsymbol{g}]\\!]_{F}:\\{\\!\\!\\{\nabla\boldsymbol{v}\\}\\!\\!\\}_{F}\mathrm{d}s.$

Then, given a barrier parameter $\mu$ and a penalization parameter $\sigma>0$,
the discretized deflated barrier subproblem is to find
$(\rho_{h},\boldsymbol{u}_{h},p_{h},\lambda_{h})\in
C_{\gamma,h}\times\boldsymbol{X}_{\mathrm{BDM}^{\boldsymbol{g}}_{k}}\times
X_{\mathrm{DG}_{k-1}}\times\mathbb{R}$ such that, for all
$(\eta_{h},\boldsymbol{v}_{h},q_{h},\zeta_{h})\in
C_{[0,1],h}\times\boldsymbol{X}_{\mathrm{BDM}^{\boldsymbol{0}}_{k}}\times(X_{\mathrm{DG}_{k-1}}\backslash\mathbb{R})\times\mathbb{R}$,
we have:

(18) $\displaystyle
c(\rho_{h},\eta_{h}-\rho_{h};\boldsymbol{u}_{h},\lambda_{h})$
$\displaystyle\geq 0,$ (19) $\displaystyle
a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h};\rho_{h})+b(\boldsymbol{v}_{h},p_{h})$
$\displaystyle=l_{h}(\boldsymbol{v}_{h}),$ (20) $\displaystyle
b(\boldsymbol{u}_{h},q_{h})$ $\displaystyle=0,$ (21) $\displaystyle
d(\lambda_{h},\zeta_{h};\rho_{h})$ $\displaystyle=0.$

Let the broken $H^{1}$-norm
$\|\cdot\|^{2}_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}$ be defined as

(22)
$\displaystyle\begin{split}\|\boldsymbol{v}\|^{2}_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}&\coloneqq\|\boldsymbol{v}\|^{2}_{L^{2}(\Omega)}+\sum_{K\in\mathcal{T}_{h}}\|\nabla\boldsymbol{v}\|^{2}_{L^{2}(K)}\\\
&\indent+\sum_{F\in\mathcal{F}^{i}_{h}}\int_{F}h_{F}^{-1}|[\\![\boldsymbol{v}]\\!]_{F}|^{2}\mathrm{d}s+\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}h_{F}^{-1}|[\\![\boldsymbol{v}-\boldsymbol{g}]\\!]_{F}|^{2}\mathrm{d}s.\end{split}$

Building on previous work [52], it was shown by Papadopoulos [49] that, for
every isolated minimizer $(\boldsymbol{u},\rho,p)$ of Eq. BP, there exists a
sequence of discretized solutions
$(\boldsymbol{u}_{h},\rho_{h},p_{h},\lambda_{h})$ to Eq. 18–Eq. 21, such that,
as $h\to 0$,
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$, $\|\rho-\rho_{h}\|_{L^{s}(\Omega)}\to 0$, $s\in[1,\infty)$, and
$\|p-p_{h}\|_{L^{2}(\Omega)}\to 0$. By extrapolating known results of a BDM
discretization for the Stokes and Stokes–Brinkman equations [16, 17, 44] we
expect a first-order $\mathrm{BDM}_{1}$ discretization to converge at a rate
of $\mathcal{O}(h)$ in the broken $H^{1}$-norm for the velocity, the
$L^{2}$-norm for the pressure, and the $L^{2}$-norm for the material
distribution. Numerical evidence for these rates of convergence is given in
[49, Fig. 2].

######  Remark 4.1.

The value of $\alpha(\rho)$ in the momentum equation Eq. 7 can range between 0
and $\bar{\alpha}$ where in practice $\bar{\alpha}\sim\mathcal{O}(10^{4})$.
Hence, an element that is robust to the transition between Stokes and Darcy
flow may have better accuracy, e.g. the Mardal–Tai–Winther (MTW) finite
element [47]. The MTW finite element also produces pointwise divergence-free
solutions and the characterization of the kernel of the divergence operator is
known [47, Sec. 4.2]. In particular, it has been shown that vertex-star patch
relaxation is also effective for multigrid cycles involving MTW finite
elements [8, Sec. 7.2]. Hence, we believe the preconditioner described below
would still be effective. However, unlike BDM finite elements with an interior
penalty, there are currently no convergence results for an MTW finite element
discretization to minimizers of the Borrvall–Petersson problem.

###  4.2 The Benson–Munson linear system

The BM solver [12] attempts to find a solution of a complementarity problem
via linearizations of the residual constrained to the inactive set. First, the
discrete Newton system is formed and the active set is defined. The active set
contains the degrees of freedom that satisfy a strict complementarity
condition in the primal and residual vectors. Next, the rows and columns of
the Jacobian in the Newton system associated with the active set degrees of
freedom are set to those of the identity. Finally, the rows on the right-hand
side vector associated with the active set degrees of freedom are fixed to
zero. Once the update, $\delta\boldsymbol{\mathbf{z}}$, of this modified
system is computed, the new iterate
$\boldsymbol{\mathbf{z}}^{k+1}=\boldsymbol{\mathbf{z}}^{k}+\delta\boldsymbol{\mathbf{z}}$
is component-wise projected onto the box constraints.

Let $\boldsymbol{\mathbf{f}}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and consider the
box constraints
$\boldsymbol{\mathbf{a}},\boldsymbol{\mathbf{b}}\in\mathbb{R}^{n}$ where
$\boldsymbol{\mathbf{a}}_{i}<\boldsymbol{\mathbf{b}}_{i}$ for all
$i=1,\dots,n$. Consider the mixed complementarity problem given by

(23) either
$\displaystyle\;\;\boldsymbol{\mathbf{a}}_{i}<\boldsymbol{\mathbf{z}}_{i}<\boldsymbol{\mathbf{b}}_{i}\;\;$
$\displaystyle\text{and}\;\;\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}})_{i}=0,$
(24) or
$\displaystyle\;\;\boldsymbol{\mathbf{a}}_{i}=\boldsymbol{\mathbf{z}}_{i}\;\;$
$\displaystyle\text{and}\;\;\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}})_{i}\geq
0,$ (25) or
$\displaystyle\;\;\boldsymbol{\mathbf{z}}_{i}=\boldsymbol{\mathbf{b}}_{i}\;\;$
$\displaystyle\text{and}\;\;\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}})_{i}\leq
0.$

We define
$\hat{\boldsymbol{\mathbf{f}}}:\\{\boldsymbol{\mathbf{z}}\in\mathbb{R}^{n}:\boldsymbol{\mathbf{a}}_{i}\leq\boldsymbol{\mathbf{z}}_{i}\leq\boldsymbol{\mathbf{b}}_{i},\;i=1,\dots,n\\}\to\mathbb{R}^{n}$
as follows, for $i=1,\dots,n$,

(26)
$\displaystyle[\hat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{z}})]_{i}\coloneqq\begin{cases}\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}})_{i}&\text{if}\;\;\boldsymbol{\mathbf{a}}_{i}<\boldsymbol{\mathbf{z}}_{i}<\boldsymbol{\mathbf{b}}_{i},\\\
\min\\{\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}})_{i},0\\}&\text{if}\;\;\boldsymbol{\mathbf{z}}_{i}=\boldsymbol{\mathbf{a}}_{i},\\\
\max\\{\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}})_{i},0\\}&\text{if}\;\;\boldsymbol{\mathbf{z}}_{i}=\boldsymbol{\mathbf{b}}_{i}.\end{cases}$

We denote the component-wise projection onto the true box constraints by
$\pi$, i.e.

(27)
$\displaystyle\pi(\boldsymbol{\mathbf{z}})_{i}=\begin{cases}\boldsymbol{\mathbf{a}}_{i}&\text{if}\;\boldsymbol{\mathbf{z}}_{i}<\boldsymbol{\mathbf{a}}_{i},\\\
\boldsymbol{\mathbf{z}}_{i}&\text{if}\;\boldsymbol{\mathbf{a}}_{i}\leq\boldsymbol{\mathbf{z}}_{i}\leq\boldsymbol{\mathbf{b}}_{i},\\\
\boldsymbol{\mathbf{b}}_{i}&\text{if}\;\boldsymbol{\mathbf{z}}_{i}>\boldsymbol{\mathbf{b}}_{i}.\end{cases}$

Consider a matrix $\boldsymbol{\mathsf{A}}\in\mathbb{R}^{n\times n}$ and the
subset of indices $S\subset\\{1,\dots,n\\}$. Then, the matrix
$\boldsymbol{\mathsf{A}}_{S,S}\in\mathbb{R}^{n\times n}$ is defined by zeroing
the rows and columns of $\boldsymbol{\mathsf{A}}$ in
$\\{1,\dots,n\\}\backslash S$ and placing a 1 on the diagonal of the zeroed
rows and columns. Similarly for any column vector
$\boldsymbol{\mathbf{x}}\in\mathbb{R}^{n}$, the vector
$\boldsymbol{\mathbf{x}}_{S}\in\mathbb{R}^{n}$ is constructed by zeroing the
rows in the set $\\{1,\dots,n\\}\backslash S$ of $\boldsymbol{\mathbf{x}}$.

######  Remark 4.2.

In the original BM solver found in [12, Sec. 3.2], the rows and columns are
eliminated rather than zeroed, resulting in a reduced matrix. However, the BM
solver as described in Algorithm 1 and [12, Sec. 3.2] are equivalent.

Algorithm 1 Benson–Munson solver [12, Sec. 3.2]

1:Initialize:

2: $k\leftarrow 0$ $\triangleright$ Initial iteration number
$\boldsymbol{\mathbf{a}},\boldsymbol{\mathbf{b}}$ $\triangleright$ Box
constraints
$\boldsymbol{\mathbf{z}}^{0}$ $\triangleright$ Feasible initial guess
tol $\triangleright$ Approximate solve tolerance

3:while
$\|\hat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{z}}^{k})\|_{2}>\mathrm{tol}$
do

4: Define the active set $\mathcal{A}^{k}$ and the inactive set
$\mathcal{I}^{k}$ as:  $\displaystyle\mathcal{A}^{k}$
$\displaystyle\coloneqq\\{i:\boldsymbol{\mathbf{z}}^{k}_{i}=\boldsymbol{\mathbf{a}}_{i}\;\text{and}\;\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}}^{k})_{i}>0\\}\cup\\{i:\boldsymbol{\mathbf{z}}^{k}_{i}=\boldsymbol{\mathbf{b}}_{i}\;\text{and}\;\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}}^{k})_{i}<0\\},$
$\displaystyle\mathcal{I}^{k}$
$\displaystyle\coloneqq\\{i\\}_{i=1}^{n}\backslash\mathcal{A}^{k}.$

5: Solve the BM linear system:
$\boldsymbol{\mathbf{f}}^{\prime}(\boldsymbol{\mathbf{z}}^{k})_{\mathcal{I}^{k},\mathcal{I}^{k}}\delta\boldsymbol{\mathbf{z}}=-\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{z}}^{k})_{\mathcal{I}^{k}}$.

6:
$\boldsymbol{\mathbf{z}}^{(k+1)}\leftarrow\pi(\boldsymbol{\mathbf{z}}^{k}+\beta^{k}\delta\boldsymbol{\mathbf{z}})$,
where $\beta^{k}$ is a (possibly adaptive) linesearch.

7: $k\leftarrow k+1$

8:end while

When applied to a linear elliptic problem, and if the active and inactive sets
of the two algorithms coincide, the updates computed by BM and Hintermüller et
al.’s primal-dual active set strategy [39] are identical [51, Th. A.1]. Under
suitable conditions, the primal-dual active set strategy can be shown to be a
semismooth Newton method [55, 56, 66]. Although the problems we consider here
do not satisfy the assumptions made in Theorem A.1 of [51], we do observe
local superlinear convergence in practice.

We now derive the linear systems that arise when using the BM active set
strategy to solve the nonlinear system Eq. 18–Eq. 21. Denote the basis
functions of the finite element spaces of the material distribution, the
velocity, the pressure, and $\mathbb{R}$ by
$\eta_{i},\boldsymbol{\phi}_{i},\psi_{i}$, and $r$, respectively and the
number of degrees of freedom by $n_{\rho},n_{\boldsymbol{u}}$, $n_{p}$, and 1,
respectively, so that the total number of degrees of freedom of the system is
given by $n=n_{\rho}+n_{\boldsymbol{u}}+n_{p}+1$. Consider the finite element
BM iterate
$\boldsymbol{z}^{k}_{h}=(\rho^{k}_{h},\boldsymbol{u}^{k}_{h},p^{k}_{h},\lambda^{k}_{h})$.
Let $\boldsymbol{\mathbf{f}}^{k}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ denote the
nonlinear residual induced by the complementarity reformulation of Eq. 18–Eq.
21. In the following, we drop the superscript iteration number k for clarity.
Let $\boldsymbol{\mathbf{z}}$ denote the discrete coefficient vector of
$\boldsymbol{z}_{h}$. In the context of Algorithm 1, the box constraints take
the values $\boldsymbol{\mathbf{a}}_{i}=0$ and $\boldsymbol{\mathbf{b}}_{i}=1$
for all degrees of freedom associated with $\rho_{h}$ and
$\boldsymbol{\mathbf{a}}_{i}=-\infty$, $\boldsymbol{\mathbf{b}}_{i}=+\infty$,
otherwise. The BM linear system, as solved on line 4 of Algorithm 1, is the
following:

(28)
$\displaystyle\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}\delta\boldsymbol{\mathbf{z}}=\begin{pmatrix}\boldsymbol{\mathsf{C}}_{\mu}&\boldsymbol{\mathsf{D}}^{\top}&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{E}}^{\top}\\\
\boldsymbol{\mathsf{D}}&\boldsymbol{\mathsf{A}}&\boldsymbol{\mathsf{B}}^{\top}&\boldsymbol{\mathsf{0}}\\\
\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{B}}&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{0}}\\\
\boldsymbol{\mathsf{E}}&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{0}}\end{pmatrix}\begin{pmatrix}\delta\boldsymbol{\mathbf{\rho}}\\\
\delta\boldsymbol{\mathbf{u}}\\\ \delta\boldsymbol{\mathbf{p}}\\\
\delta\boldsymbol{\mathbf{\lambda}}\end{pmatrix}=-\begin{pmatrix}\boldsymbol{\mathbf{f}}_{\rho}\\\
\boldsymbol{\mathbf{f}}_{\boldsymbol{u}}\\\ \boldsymbol{\mathbf{f}}_{p}\\\
\boldsymbol{\mathbf{f}}_{\lambda}\end{pmatrix}=-\boldsymbol{\mathbf{f}},$

where
$\delta\boldsymbol{\mathbf{\rho}},\delta\boldsymbol{\mathbf{u}},\delta\boldsymbol{\mathbf{p}}$
and $\delta\boldsymbol{\mathbf{\lambda}}$ denote the discrete coefficient
vector BM updates for $\rho,\boldsymbol{u},p$ and $\lambda$, and
$\boldsymbol{\mathbf{f}}_{\rho}$, $\boldsymbol{\mathbf{f}}_{\boldsymbol{u}}$,
$\boldsymbol{\mathbf{f}}_{p}$, and $\boldsymbol{\mathbf{f}}_{\lambda}$ are the
corresponding blocks of the nonlinear residual with the active set rows,
$i\in\mathcal{A}$, in $\boldsymbol{\mathbf{f}}_{\rho}$ zeroed. The entries of
$\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ are given by

(29)
$\displaystyle[\boldsymbol{\mathsf{A}}]_{ij}=a_{h}(\boldsymbol{\phi}_{j},\boldsymbol{\phi}_{i};\rho_{h})\text{
and }[\boldsymbol{\mathsf{B}}]_{ij}=b(\boldsymbol{\phi}_{j},\psi_{i}).$

Furthermore, if $j\in\mathcal{I}$, then

(30)
$\displaystyle[\boldsymbol{\mathsf{D}}]_{ij}=\int_{\Omega}(\alpha^{\prime}(\rho_{h})\boldsymbol{u}_{h}\cdot\boldsymbol{\phi}_{i})\eta_{j}\;\mathrm{d}x,\;\;\;[\boldsymbol{\mathsf{E}}]_{ij}=-r\int_{\Omega}\eta_{j}\;\mathrm{d}x.$

Otherwise if $j\in\mathcal{A}$, then $[\boldsymbol{\mathsf{D}}]_{ij}=0$ and
$[\boldsymbol{\mathsf{E}}]_{ij}=0$ for all $i$. Finally if $i,j\in\mathcal{I}$
then

(31)
$\displaystyle[\boldsymbol{\mathsf{C}}_{\mu}]_{ij}=\int_{\Omega}\left[\frac{1}{2}\alpha^{\prime\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}+\frac{\mu}{(\rho_{h}+\epsilon_{\text{log}})^{2}}+\frac{\mu}{(1+\epsilon_{\text{log}}-\rho_{h})^{2}}\right]\eta_{i}\eta_{j}\mathrm{d}x.$

Otherwise, if $i\in\mathcal{A}$ or $j\in\mathcal{A}$, then
$[\boldsymbol{\mathsf{C}}_{\mu}]_{ij}=\delta_{ij}$, where $\delta_{ij}$ is the
Kronecker delta.

######  Remark 4.3.

$\boldsymbol{\mathsf{E}}$ is a row vector of size $1\times n_{\rho}$.

In the remainder of this subsection we discuss the invertibility of the matrix
in Eq. 28 and its subblocks.

######  Proposition 4.4.

The matrix $\boldsymbol{\mathsf{A}}\in\mathbb{R}^{n_{\boldsymbol{u}}\times
n_{\boldsymbol{u}}}$ is symmetric, and provided the penalization parameter
$\sigma>0$ is sufficiently large, then it also positive-definite.

######  Proof 4.5.

Symmetry is realized by swapping the indices $i$ and $j$ of the basis
functions in their respective definitions and noting that the resulting
integrals are equal. Positive-definiteness of $\boldsymbol{\mathsf{A}}$, for
sufficiently large $\sigma>0$, follows from $\alpha(\rho)\geq 0$ and [40, Sec.
3.3].

######  Proposition 4.6.

The matrix $\boldsymbol{\mathsf{C}}_{\mu}$ is symmetric positive semi-
definite. Moreover, if either $\mu>0$ or $|\boldsymbol{u}_{h}|>0$ a.e., then
$\boldsymbol{\mathsf{C}}_{\mu}$ is symmetric positive-definite.

######  Proof 4.7.

The symmetry of $\boldsymbol{\mathsf{C}}_{\mu}$ is realized by swapping the
indices $i$ and $j$ of the basis functions in its definition and noting that
the resulting integrals are equal. Consider the unmodified matrix
$\hat{\boldsymbol{\mathsf{C}}}_{\mu}$, defined by Eq. 31 for all $i$ and $j$.
Pick an arbitrary function $\eta_{h}\in X_{\mathrm{DG}_{0}}$ with discrete
coefficient vector $\boldsymbol{\mathbf{\eta}}\in\mathbb{R}^{n_{\rho}}$. We
note that

(32)
$\boldsymbol{\mathbf{\eta}}^{\top}\hat{\boldsymbol{\mathsf{C}}}_{\mu}\boldsymbol{\mathbf{\eta}}=\int_{\Omega}\left[\frac{1}{2}\alpha^{\prime\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}+\frac{\mu}{(\rho_{h}+\epsilon_{\text{log}})^{2}}+\frac{\mu}{(1+\epsilon_{\text{log}}-\rho_{h})^{2}}\right]|\eta_{h}|^{2}\mathrm{d}x.$

All the terms in the integral Eq. 32 are non-negative. Hence,
$\hat{\boldsymbol{\mathsf{C}}}_{\mu}$ is positive semi-definite.

Assumption (A5) implies that $\alpha^{\prime\prime}(\rho_{h})>0$. If $\mu>0$,
then the rational expressions are strictly greater than zero as
$0\leq\rho_{h}\leq 1$. Otherwise, if $|\boldsymbol{u}_{h}|>0$ a.e., then the
term $\alpha^{\prime\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}>0$ a.e. Hence,
if either $\mu>0$ or $|\boldsymbol{u}_{h}|>0$ a.e. the right-hand side of Eq.
32 is equal to zero if and only if $\eta_{h}=0$, which is true if and only if
$\boldsymbol{\mathbf{\eta}}=\boldsymbol{\mathbf{0}}$. Therefore, if either
$\mu>0$ or $|\boldsymbol{u}_{h}|>0$ a.e.,
$\boldsymbol{\mathbf{\eta}}^{\top}\hat{\boldsymbol{\mathsf{C}}}_{\mu}\boldsymbol{\mathbf{\eta}}\geq
0$ with equality if and only if
$\boldsymbol{\mathbf{\eta}}=\boldsymbol{\mathbf{0}}$. Hence,
$\hat{\boldsymbol{\mathsf{C}}}_{\mu}$ is symmetric positive-definite.

Since the discretization for $\rho$ is piecewise constant,
$\hat{\boldsymbol{\mathsf{C}}}_{\mu}$ is a diagonal matrix and, therefore, all
diagonal entries must be positive. The procedure of zeroing rows and columns
associated with the BM active set and replacing the diagonal entry with a one
will result in a diagonal matrix with non-negative (positive if either $\mu>0$
or $|\boldsymbol{u}_{h}|>0$ a.e.) diagonal entries. We conclude that
$\boldsymbol{\mathsf{C}}_{\mu}$ must be symmetric positive semi-definite and
if either $\mu>0$ or $\boldsymbol{u}_{h}>0$ a.e., then it is symmetric
positive-definite.

######  Remark 4.8.

For problems with a pure Dirichlet boundary condition,
$\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}$ will have a
nullspace of at least one dimension associated with the fact that the pressure
is only unique up to a constant [22, Ch. 4]. Since we will use an FGMRES
Krylov method on the outermost level of the linear solve, this particular
nullspace can be handled by FGMRES without the method breaking down [22, Ch.
9.3.5]. To show that the nullspace is indeed the space of constants with
respect to the pressure, we assume that $\boldsymbol{\mathsf{B}}$ has been
modified so that it is of full row rank, e.g. by picking a pressure degree of
freedom and fixing it to zero.

######  Proposition 4.9.

Suppose that either $\mu>0$ or $|\boldsymbol{u}_{h}|>0$ a.e. Moreover, assume
that
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is symmetric positive-definite and $\boldsymbol{\mathsf{B}}$ has been modified
so that it has full row rank. Then, the matrix
$\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}$, as defined in Eq.
28, is invertible.

######  Proof 4.10.

If either $\mu>0$ or $|\boldsymbol{u}_{h}|>0$ a.e. then
$\boldsymbol{\mathsf{C}}_{\mu}$ is symmetric positive-definite by Proposition
4.6 and, therefore, invertible. Thus the matrix
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is well-defined. By assumption
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is symmetric positive-definite. Therefore, since $\boldsymbol{\mathsf{A}}$ is
symmetric positive-definite by Proposition 4.4, we have that the following
matrix is also symmetric positive-definite [75, Th. 1.12]:

(33)
$\displaystyle\boldsymbol{\mathsf{G}}\coloneqq\begin{pmatrix}\boldsymbol{\mathsf{C}}_{\mu}&\boldsymbol{\mathsf{D}}^{\top}\\\
\boldsymbol{\mathsf{D}}&\boldsymbol{\mathsf{A}}\end{pmatrix}.$

We define
$\boldsymbol{\mathsf{B}}_{0}\coloneqq\begin{pmatrix}\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{B}}\end{pmatrix}\in\mathbb{R}^{n_{p}\times(n_{\rho}+n_{\boldsymbol{u}})}$
and
$\boldsymbol{\mathsf{E}}_{0}\coloneqq\begin{pmatrix}\boldsymbol{\mathsf{E}}&\boldsymbol{\mathsf{0}}\end{pmatrix}\in\mathbb{R}^{1\times(n_{\rho}+n_{\boldsymbol{u}})}$
and re-block $\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}$ as
follows:

(34)
$\displaystyle\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}=\begin{pmatrix}\boldsymbol{\mathsf{G}}&\boldsymbol{\mathsf{B}}_{0}^{\top}&\boldsymbol{\mathsf{E}}^{\top}_{0}\\\
\boldsymbol{\mathsf{B}}_{0}&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{0}}\\\
\boldsymbol{\mathsf{E}}_{0}&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{0}}\end{pmatrix}.$

The matrix on the right-hand side of Eq. 34 is of double-saddle point type.
$\boldsymbol{\mathsf{E}}_{0}$ is a nonzero row vector and thus must have full
row rank. Moreover, by assumption $\boldsymbol{\mathsf{B}}$ has been modified
so that it has full row rank. Hence, $\boldsymbol{\mathsf{B}}_{0}$ also has
full row rank. Consider any vector
$\boldsymbol{\mathbf{x}}=(\boldsymbol{\mathbf{x}}_{1}^{\top}\;\;\;\boldsymbol{\mathbf{x}}_{2}^{\top})^{\top}\in\mathbb{R}^{n_{\rho}+n_{\boldsymbol{u}}}$.
We note that
$\boldsymbol{\mathsf{E}}_{0}^{\top}\boldsymbol{\mathbf{x}}^{\top}=\boldsymbol{\mathsf{E}}^{\top}\boldsymbol{\mathbf{x}}_{1}^{\top}$
and
$\boldsymbol{\mathsf{B}}_{0}^{\top}\boldsymbol{\mathbf{x}}^{\top}=\boldsymbol{\mathsf{B}}^{\top}\boldsymbol{\mathbf{x}}_{2}^{\top}$.
Hence,
$\mathrm{range}(\boldsymbol{\mathsf{B}}_{0}^{\top})\cap\mathrm{range}(\boldsymbol{\mathsf{E}}_{0}^{\top})=\\{\boldsymbol{0}\\}$.
Thus $\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}$ is invertible
[5, Prop. 2.3].

######  Remark 4.11.

The proof of the invertibility of
$\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}$ (provided the
pressure null space has been eliminated) relied heavily on the assumption that
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is symmetric positive-definite. Although the symmetry is guaranteed, positive-
definiteness is not. Numerically, we found that
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is often positive-definite for large ranges of $\mu$, choices of $\alpha$, and
BM iterates $\rho_{h,k}$ and $\boldsymbol{u}_{h,k}$. Despite this, there are
choices where
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
has negative eigenvalues. In practice we have not encountered a case where
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is singular.

######  Remark 4.12.

A full proof of the invertibility of
$\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}$ in the case where
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is indefinite or singular is beyond the scope of this work. However, we note
the following. The first-order optimality conditions Eq. 6–Eq. 9 can be
rewritten as a semismooth system of equations. Suppose, after discretization,
we denote this system
$F:\mathbb{R}^{n_{\boldsymbol{u}}+3n_{\rho}+n_{p}+1}\to\mathbb{R}^{n_{\boldsymbol{u}}+3n_{\rho}+n_{p}+1}$
where the extra $2n_{\rho}$ degrees of freedom are associated with the
Lagrange multipliers that enforce the box constraints on $\rho$. Then,
provided one can show that $F$ is locally Lipschitz continuous, $F$ is
semismooth and $F^{\prime}(\boldsymbol{\mathbf{z}}_{*})$ is invertible at
$\boldsymbol{\mathbf{z}}_{*}$ where $\boldsymbol{\mathbf{z}}_{*}$ satisfies
$F(\boldsymbol{\mathbf{z}}_{*})=0$, then one may be able to invoke a
semismooth analogue of the Rall–Rheinboldt theory [27]. Hence, provided the
iterate $\boldsymbol{\mathbf{z}}_{k}$ is sufficiently close to
$\boldsymbol{\mathbf{z}}_{*}$, then $F^{\prime}(\boldsymbol{\mathbf{z}}_{k})$
is invertible. However, we note that although the BM solver has been shown to
be almost equivalent to a semismooth Newton method for linear elliptic
problems in [51, App. A], that theory does not extend to the systems we are
solving here. Hence the invertibility of
$F^{\prime}(\boldsymbol{\mathbf{z}}_{k})$ would not necessarily imply the
invertibility of $\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p,\lambda}$.

##  5 Preconditioning

In this section, we develop a preconditioner for solving Eq. 28. As discussed
in Section 3.2, preconditioning strategies that are robust for the undeflated
system can also be used to compute solutions of the deflated systems. On the
outermost level of the deflated barrier method, we perform continuation in the
barrier parameter $\mu$. Next, at a given $\mu$, we use the BM solver to find
a solution of Eq. 18–Eq. 21. A direct sparse LU factorization of the matrix in
Eq. 28 is infeasible on fine meshes of three-dimensional problems. Thus we
turn to preconditioning techniques to reduce the cost of each inner linear
solve. The preconditioning is made difficult by the saddle point nature of the
matrix in Eq. 28 and the barrier-like terms in
$\boldsymbol{\mathsf{C}}_{\mu}$. In the following subsections we introduce a
nested block preconditioning method for solving Eq. 28, where the Schur
complements are controlled with an augmented Lagrangian term. As outermost
Krylov solver, we use a preconditioned FGMRES method [59].

###  5.1 Block preconditioning

Consider the well-posed linear system

(35) $\displaystyle\begin{pmatrix}\mathbb{A}&\mathbb{B}\\\
\mathbb{C}&\mathbb{D}\end{pmatrix}\begin{pmatrix}\boldsymbol{\mathbf{x}}\\\
\boldsymbol{\mathbf{y}}\end{pmatrix}=\begin{pmatrix}\boldsymbol{\mathbf{c}}\\\
\boldsymbol{\mathbf{d}}\end{pmatrix},$

where $\mathbb{A}\in\mathbb{R}^{n_{1}\times n_{1}}$ is invertible,
$\mathbb{B}\in\mathbb{R}^{n_{1}\times n_{2}}$,
$\mathbb{C}\in\mathbb{R}^{n_{2}\times n_{1}}$ and
$\mathbb{D}\in\mathbb{R}^{n_{2}\times n_{2}}$. Then, under suitable conditions
[13, §3.2] the inverse of the matrix in Eq. 35 admits a full block
factorization of the form

(36) $\displaystyle\begin{pmatrix}\mathbb{A}&\mathbb{B}\\\
\mathbb{C}&\mathbb{D}\end{pmatrix}^{-1}=\begin{pmatrix}I&-\mathbb{A}^{-1}\mathbb{B}\\\
0&I\end{pmatrix}\begin{pmatrix}\mathbb{A}^{-1}&0\\\
0&\mathbb{S}^{-1}\end{pmatrix}\begin{pmatrix}I&0\\\
-\mathbb{C}\mathbb{A}^{-1}&I\end{pmatrix},$

where $\mathbb{S}=\mathbb{D}-\mathbb{C}\mathbb{A}^{-1}\mathbb{B}$.
Preconditioners for Eq. 35 can be found by developing cheap approximations to
$\mathbb{A}^{-1}$ and $\mathbb{S}^{-1}$ and substituting them into Eq. 36 [48,
69].

The subspace spanned by the volume constraint Lagrange multiplier $\lambda$ is
one-dimensional and can be handled by at most one iteration of a Krylov
subspace solver or via block preconditioning. Experimentally, we found that a
full block preconditioner of the real block performed best. Writing the
density-momentum-pressure block,

(37)
$\displaystyle\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p}\coloneqq\begin{pmatrix}\boldsymbol{\mathsf{C}}_{\mu}&\boldsymbol{\mathsf{D}}^{\top}&\boldsymbol{\mathsf{0}}\\\
\boldsymbol{\mathsf{D}}&\boldsymbol{\mathsf{A}}&\boldsymbol{\mathsf{B}}^{\top}\\\
\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{B}}&\boldsymbol{\mathsf{0}}\end{pmatrix},$

we choose $\mathbb{A}=\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p}$,
$\mathbb{B}=\boldsymbol{\mathsf{E}}^{\top}$,
$\mathbb{C}=\boldsymbol{\mathsf{E}}$, and
$\mathbb{D}=\boldsymbol{\mathsf{0}}$. The Schur complement,

(38)
$\displaystyle\mathbb{S}=\boldsymbol{\mathsf{S}}_{0}\coloneqq-\boldsymbol{\mathsf{E}}\boldsymbol{\mathsf{H}}_{\rho,\boldsymbol{u},p}^{-1}\boldsymbol{\mathsf{E}}^{\top},$

is a $1\times 1$ matrix and can be inverted by taking its reciprocal. Hence,
the difficulty now lies in solving linear systems involving Eq. 37. Since we
have only decreased the size of the linear system by one dimension, an LU
factorization is still infeasible and we consider block preconditioners for
Eq. 37. We summarize the initial components of the solver in Fig. 1.

[width=]figures/preconditioning/flowchart-beginning

Figure 1: Main components of the deflated barrier method solver. The remainder
of this section focuses on developing preconditioners for the item in blue.

We precondition (37) by taking the Schur complement with respect to the
momentum-pressure block. This approach was utilized by Evgrafov for
preconditioning the linear systems arising in a similar solver [25, Sec. 5].
In the notation of Eq. 35, $\mathbb{A}=\boldsymbol{\mathsf{C}}_{\mu}$,
$\mathbb{B}=(\boldsymbol{\mathsf{D}}^{\top}\;\;\boldsymbol{\mathsf{0}})$,
$\mathbb{C}=(\boldsymbol{\mathsf{D}}\;\;\boldsymbol{\mathsf{0}})^{\top}$ and,

(39)
$\displaystyle\mathbb{D}=\boldsymbol{\mathsf{Q}}\coloneqq\begin{pmatrix}\boldsymbol{\mathsf{A}}&\boldsymbol{\mathsf{B}}^{\top}\\\
\boldsymbol{\mathsf{B}}&\boldsymbol{\mathsf{0}}\end{pmatrix}.$

We know that $\boldsymbol{\mathsf{C}}_{\mu}$ is invertible by Proposition 4.6.
Hence, we write $\mathbb{A}^{-1}=\boldsymbol{\mathsf{C}}_{\mu}^{-1}$. The
matrix $\boldsymbol{\mathsf{Q}}$ resembles the linear system that arises in
the discretization of the Stokes–Brinkman equations. The Schur complement is
given by

(40)
$\displaystyle\mathbb{S}=\boldsymbol{\mathsf{S}}_{1}\coloneqq\boldsymbol{\mathsf{Q}}-\begin{pmatrix}\boldsymbol{\mathsf{D}}\\\
\boldsymbol{\mathsf{0}}\end{pmatrix}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\begin{pmatrix}\boldsymbol{\mathsf{D}}^{\top}&\boldsymbol{\mathsf{0}}\end{pmatrix}=\begin{pmatrix}\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}&\boldsymbol{\mathsf{B}}^{\top}\\\
\boldsymbol{\mathsf{B}}&\boldsymbol{\mathsf{0}}\end{pmatrix}.$

The reason we use a $\mathrm{DG}_{0}$ piecewise constant discretization for
the material distribution is to ensure that $\boldsymbol{\mathsf{S}}_{1}$ is
sparse. Since the material distribution is discretized with $\mathrm{DG}_{0}$
finite elements, $\boldsymbol{\mathsf{C}}_{\mu}$ is a diagonal matrix,
$\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is still sparse and hence $\boldsymbol{\mathsf{S}}_{1}$ is also sparse.

Thus far, the application of block preconditioning has reduced the solution of
the full matrix Eq. 28 to the following:

  1. 1.

Apply an outer FGMRES solver;

  2. 2.

Apply the reciprocal of $\boldsymbol{\mathsf{S}}_{0}\in\mathbb{R}$;

  3. 3.

Invert the diagonal matrix $\boldsymbol{\mathsf{C}}_{\mu}$;

  4. 4.

Apply the inverse of the $2\times 2$ block matrix
$\boldsymbol{\mathsf{S}}_{1}$, which is the same size as the matrix that
arises in a discretized pure Stokes problem.

For now we assume that $\boldsymbol{\mathsf{S}}_{1}$ is invertible and defer
discussions of invertibility to later.

We must now develop solvers for $\boldsymbol{\mathsf{S}}_{1}$ as given in Eq.
40. One option is to use a direct solver. However, we can further reduce the
computational work with another application of block preconditioning. Consider
taking the inner Schur complement in $\boldsymbol{\mathsf{S}}_{1}$ with
respect to the pressure block. Using the notation of Eq. 35,
$\mathbb{B}=\boldsymbol{\mathsf{B}}^{\top}$,
$\mathbb{C}=\boldsymbol{\mathsf{B}}$, and
$\mathbb{A}=\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$.
The innermost Schur complement takes the form

(41)
$\displaystyle\mathbb{S}=\boldsymbol{\mathsf{S}}_{2}\coloneqq-\boldsymbol{\mathsf{B}}(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top})^{-1}\boldsymbol{\mathsf{B}}^{\top}.$

This time, $\boldsymbol{\mathsf{S}}_{2}$ is dense, and we employ an augmented
Lagrangian approach.

We first recall some known results for the Stokes and Stokes–Brinkman
equations. Let $-\Delta_{h}$ denote the negative discretized Laplacian matrix.
In the context of the incompressible Stokes equations,
$\boldsymbol{\mathsf{B}}(-\Delta_{h})^{-1}\boldsymbol{\mathsf{B}}^{\top}$ is
spectrally equivalent to the viscosity-scaled pressure mass matrix
$\nu^{-1}\boldsymbol{\mathsf{M}}_{p}$ [63, 68], see also [22, Th. 5.22].
$\boldsymbol{\mathsf{M}}_{p}$ is a sparse mass matrix and can be cheaply
factorized or solved with a multigrid method. Therefore, a good approximation
to the Schur complement of the pure Newtonian Stokes problem is given by
$\nu^{-1}\boldsymbol{\mathsf{M}}_{p}$. The idea is that the momentum block can
then be solved with a direct solver, multigrid methods, or other alternative
solvers. Unfortunately, in the context of the Stokes–Brinkman equations, Popov
[54] noted that the presence of the Brinkman term
$\alpha(\rho_{h})\boldsymbol{u}_{h}$ in the momentum block
$\boldsymbol{\mathsf{A}}$ renders the approximation given by
$\nu^{-1}\boldsymbol{\mathsf{M}}_{p}$ ineffective. Popov proposed a Schur
complement preconditioning technique based on incomplete LU factorization, but
such factorizations do not generally yield mesh-independent preconditioners.
An alternative is a preconditioning scheme utilized by Borrvall and Petersson
in their original paper [15, Sec. 2.6] based on the work of Cahouet and
Chabard [19]. However, during numerical experiments, we found that an
augmented Lagrangian approach performed better for the problems in this work,
see [50, Ch. 5.3.1].

We now propose an augmented Lagrangian strategy to control the second Schur
complement $\boldsymbol{\mathsf{S}}_{2}$. The essential idea is to add a term
to the momentum equation that does not change the exact solution, but does
change the Schur complement $\boldsymbol{\mathsf{S}}_{2}$; in particular, by
scaling the augmentation appropriately, $\boldsymbol{\mathsf{S}}_{2}$ can be
almost perfectly approximated with a scaled pressure mass matrix. The cost of
this approach is that it makes the augmented momentum block more difficult to
solve. The augmented Lagrangian approach has been shown to be robust for a
variety of difficult saddle-point systems such as the stationary Navier–Stokes
equations at high Reynolds number [31], implicitly-constituted anisothermal
non-Newtonian flow [32], and magnetohydrodynamics [45]. Hence, this approach
has potential for extension to different fluid topology optimization problems.

There are two possible augmented Lagrangian approaches: continuous and
discrete. These approaches are mathematically equivalent for exactly
divergence-free elements such as the one employed here. We choose to introduce
the method in the discrete setting. Post-discretization, the matrix
$\boldsymbol{\mathsf{A}}$ in Eq. 28 is modified by adding an augmented
Lagrangian term

(42)
$\displaystyle\boldsymbol{\mathsf{A}}_{\gamma_{d}}\coloneqq\boldsymbol{\mathsf{A}}+\gamma_{d}\boldsymbol{\mathsf{B}}^{\top}\boldsymbol{\mathsf{M}}_{p}^{-1}\boldsymbol{\mathsf{B}},$

where $\gamma_{d}\gg 0$, and the right-hand side of Eq. 28 is modified so that
the solution of linear system remains unchanged (since
$\boldsymbol{\mathsf{B}}\boldsymbol{\mathbf{u}}$ is known). In particular, if
the current velocity iterate is divergence-free, then
$\gamma_{d}\boldsymbol{\mathsf{B}}^{\top}\boldsymbol{\mathsf{M}}_{p}^{-1}\boldsymbol{\mathsf{B}}\delta\boldsymbol{\mathbf{u}}=\boldsymbol{\mathbf{0}}$
and no modification to the right-hand side is required. While it does not
change the solution, the addition of the augmented Lagrangian term influences
the nature of the inner Schur complement. In particular,
$\boldsymbol{\mathsf{S}}_{1}$ becomes

(43)
$\displaystyle\boldsymbol{\mathsf{S}}_{1,\gamma_{d}}=\begin{pmatrix}\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}&\boldsymbol{\mathsf{B}}^{\top}\\\
\boldsymbol{\mathsf{B}}&\boldsymbol{\mathsf{0}}\end{pmatrix},$

and $\boldsymbol{\mathsf{S}}_{2}$ becomes

(44)
$\displaystyle\begin{split}\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}=-\boldsymbol{\mathsf{B}}(\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top})^{-1}\boldsymbol{\mathsf{B}}^{\top}.\end{split}$

######  Remark 5.1.

As noted in Remark 4.12, it is possible that
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is indefinite and, therefore, it is unclear if
$\boldsymbol{\mathsf{S}}_{1,\gamma_{d}}$ is invertible. If we were to
encounter such a case, then the BM linear system solve would fail. Here, the
deflated barrier method would repeatedly halve the step size in the barrier
parameter $\mu$ (keeping the initial guess the same). We did not observe this
behaviour in our numerical experiments. Henceforth, we assume that
$\boldsymbol{\mathsf{S}}_{1,\gamma_{d}}$ is invertible. We note that if
$\boldsymbol{\mathsf{S}}_{1,\gamma_{d}}$ is invertible, then
$\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}$ is invertible. This is seen by
applying Sylvester’s law of inertia to the Schur complement decomposition Eq.
36 [22, Ch. 4].

The action of $\boldsymbol{\mathsf{S}}^{-1}_{2,\gamma_{d}}$ is required during
the solve. As already mentioned, $\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}$ is
dense and, hence, it is expensive to assemble and apply its inverse. In the
next proposition we show that the action of
$\boldsymbol{\mathsf{S}}^{-1}_{2,\gamma_{d}}$ can be approximated with
$-\gamma_{d}\boldsymbol{\mathsf{M}}^{-1}_{p}$ with increasing accuracy as
$\gamma_{d}\to\infty$. In practice, values of
$\gamma_{d}\sim\mathcal{O}(10^{5})$ were sufficient to approximate the action
of $\boldsymbol{\mathsf{S}}^{-1}_{2,\gamma_{d}}$ with
$-\gamma_{d}\boldsymbol{\mathsf{M}}^{-1}_{p}$ to tolerances of
$\mathcal{O}(10^{-10})$.

######  Proposition 5.2.

Suppose that
$\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$,
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$,
$\boldsymbol{\mathsf{S}}_{2}$ and $\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}$ are
invertible and $\boldsymbol{\mathsf{B}}$ has been modified to ensure it has
full rank (see Remark 4.8). Consider the eigenvalue problem

(45)
$\displaystyle\gamma_{d}\boldsymbol{\mathsf{M}}_{p}^{-1}\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}\boldsymbol{\mathbf{q}}_{i}=\lambda_{i}\boldsymbol{\mathbf{q}}_{i},\;\;i=1,\dots,n_{p}.$

Then, $\lambda_{i}\to-1$, $i=1,\dots,n_{p}$ as $\gamma_{d}\to\infty$.

######  Proof 5.3.

We note that $\boldsymbol{\mathsf{M}}_{p}$ is symmetric positive-definite and,
therefore, invertible. Our assumptions satisfy the requirements of [14, Lem.
4.1] and thus

(46)
$\displaystyle-\boldsymbol{\mathsf{S}}^{-1}_{2,\gamma_{d}}=-\boldsymbol{\mathsf{S}}^{-1}_{2}+\gamma_{d}\boldsymbol{\mathsf{M}}_{p}^{-1}.$

By left multiplying Eq. 45 by $\gamma_{d}^{-1}\boldsymbol{\mathsf{M}}_{p}$
then $\boldsymbol{\mathsf{S}}^{-1}_{2,\gamma_{d}}$ and applying Eq. 46 we see
that

(47)
$\displaystyle\lambda_{i}\left[\boldsymbol{\mathsf{S}}_{2}^{-1}\gamma_{d}^{-1}\boldsymbol{\mathsf{M}}_{p}-\boldsymbol{\mathsf{I}}\right]\boldsymbol{\mathbf{q}}_{i}=\boldsymbol{\mathbf{q}}_{i}.$

As $\boldsymbol{\mathsf{M}}_{p}^{-1}\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}$ is
invertible, then $\lambda_{i}\neq 0$ for all $i=1,\dots,n_{p}$. Hence, by
rearranging Eq. 47 we observe that

(48)
$\displaystyle\boldsymbol{\mathsf{S}}_{2}^{-1}\boldsymbol{\mathsf{M}}_{p}\boldsymbol{\mathbf{q}}_{i}=\gamma_{d}(\lambda_{i}^{-1}+1)\boldsymbol{\mathbf{q}}_{i}.$

Since $\boldsymbol{\mathsf{S}}_{2}^{-1}\boldsymbol{\mathsf{M}}_{p}$ is
invertible, its eigenvalues (denoted $\delta_{i}$) are nonzero, and satisfy

(49)
$\displaystyle\delta_{i}=\gamma_{d}(\lambda_{i}^{-1}+1),\;\;i=1,\dots,n_{p}.$

Therefore, by taking the limit $\gamma_{d}\to\infty$ we conclude that
$\lambda_{i}\to-1$ for all $i=1,\dots,n_{p}$.

If assembled naïvely, the triple matrix product
$\boldsymbol{\mathsf{B}}^{\top}\boldsymbol{\mathsf{M}}_{p}^{-1}\boldsymbol{\mathsf{B}}$,
as it occurs in the augmented Lagrangian term, is expensive to compute.
However, it can be verified that the augmented Lagrangian term
$\gamma_{d}\boldsymbol{\mathsf{B}}^{\top}\boldsymbol{\mathsf{M}}_{p}^{-1}\boldsymbol{\mathsf{B}}$
corresponds to augmenting the weak form
$a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h};\rho_{h})$ in Eq. 19 by

(50)
$\displaystyle\gamma_{d}\int_{\Omega}\Pi(\nabla\cdot\boldsymbol{u}_{h})\Pi(\nabla\cdot\boldsymbol{v}_{h})\mathrm{d}x,$

where $\Pi$ is the projection onto the discretized pressure space. The
projection is the identity for the BDM-DG pair. Therefore, assembling
$\boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{M}}_{p}^{-1}\boldsymbol{\mathsf{B}}^{\top}$
is equivalent to assembling the matrix associated with the bilinear form
$\int_{\Omega}\nabla\cdot\boldsymbol{\phi}_{i}\nabla\cdot\boldsymbol{\phi}_{j}\mathrm{d}x$,
where $\boldsymbol{\phi}_{i}$, $i=1,\dots,n_{\boldsymbol{u}}$, are the basis
functions of the velocity finite element space.

With the proposed nested block preconditioning, we have reduced solving linear
systems involving Eq. 28 into the following steps:

  1. 1.

Apply an outer FGMRES solver;

  2. 2.

Apply the reciprocal of $\boldsymbol{\mathsf{S}}_{0}\in\mathbb{R}$;

  3. 3.

Invert the diagonal matrix $\boldsymbol{\mathsf{C}}_{\mu}$;

  4. 4.

Factorize and solve the block-diagonal pressure mass matrix
$\boldsymbol{\mathsf{M}}_{p}$;

  5. 5.

Apply the action of the inverse of the augmented momentum block
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$.

Factorizing $\boldsymbol{\mathsf{M}}_{p}$, and
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
with a direct solver such as MUMPS [6] is faster than factorizing the full
matrix in Eq. 28. We note that $\boldsymbol{\mathsf{M}}_{p}$ only needs to be
factorized once at the start of the algorithm. With ideal inner solvers, most
of the computational time during the run of the deflated barrier method is
spent on factorizing
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
at each BM iteration.

In Fig. 2, we summarize the block preconditioning strategy for solving linear
systems involving the density-momentum-pressure block Eq. 37. In the next
section we develop a specialized geometric multigrid scheme to efficiently
solve linear systems involving
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
(highlighted in pink in Fig. 2) in order to reduce the computational time
further when the problem is discretized on a fine mesh.

[width=]figures/preconditioning/flowchart-augmented-lagrangian

Figure 2: The preconditioning strategy to solve the density-momentum-pressure
block Eq. 37. We develop a geometric multigrid scheme for the item in pink in
Section 5.2.

###  5.2 A specialized multigrid scheme for
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$

As already mentioned, the tradeoff for using an augmented Lagrangian term to
control the Schur complement $\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}$ is that
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
becomes difficult to solve, due to the semi-definite term with a large
coefficient $\gamma_{d}\gg 0$. In the past two decades, there has been
progress on specialized multigrid schemes, based on the pioneering work of
Schöberl [62], to handle the effects of the augmented Lagrangian term in
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}$. Strategies based on Schöberl’s work
have been shown to be extremely effective in parameter-robust preconditioning
of the three-dimensional incompressible Navier–Stokes equations [30, 31],
Oseen–Frank models of cholesteric liquid crystals [71], implicitly-constituted
non-Newtonian incompressible flow [28, 32], and magnetohydrodynamics [45].
Related multigrid schemes have also been analyzed in the context of the
$\boldsymbol{H}(\mathrm{div};\Omega)$ and
$\boldsymbol{H}(\mathrm{curl};\Omega)$ Riesz maps [7] and, more relevant to
our problem, an $\boldsymbol{H}(\mathrm{div};\Omega)$-conforming
discretization of the Stokes equations [40].

In this work, we construct a mesh hierarchy and construct our multilevel
hierarchy geometrically. The coarse-level operators are induced from Eq.
18–Eq. 21 via rediscretization. Moreover, in all our examples, we choose a
direct solver for the coarse-level solver. Schöberl’s analysis gives
sufficient conditions on the relaxation method and transfer operators to
achieve robustness in the context of multigrid cycles applied to symmetric
positive-definite problems augmented with a parameter-dependent positive semi-
definite term. The first is that the relaxation method must stably capture the
kernel of the semi-definite term. The second requirement is that the
prolongation operator must have a continuity constant that is independent of
$\gamma_{d}$. As noted by Hong et al. [40, Sec. 1], in a nested mesh
hierarchy, an exactly divergence-free function on the coarse-grid will be
divergence-free on the fine-grid. Therefore, in our context, the natural
prolongation operator suffices thanks to our choice of discretization, and we
only discuss the relaxation method in this work. The kernel of the semi-
definite term involving $\gamma_{d}$ is

(51)
$\displaystyle\mathcal{N}_{h}=\\{\boldsymbol{w}_{h}\in\boldsymbol{X}_{\mathrm{BDM}_{k}}:(\nabla\cdot\boldsymbol{w}_{h},\nabla\cdot\boldsymbol{v}_{h})_{L^{2}(\Omega)}=0\;\text{for
all}\;\boldsymbol{v}_{h}\in\boldsymbol{X}_{\mathrm{BDM}_{k}}\\},$

i.e. all functions with divergence zero. For large $\gamma_{d}$,
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}$ becomes increasingly singular. Common
relaxation methods like Jacobi and Gauss-Seidel do not offer
$\gamma_{d}$-robust smoothing and yield ineffective multigrid cycles. To
understand the degradation of Jacobi and Gauss–Seidel as
$\gamma_{d}\to\infty$, it is fruitful to view the relaxation method as a
subspace correction method [72, 73]. Consider the space decomposition

(52)
$\displaystyle\boldsymbol{X}_{\mathrm{BDM}_{k}}=\sum_{i}\boldsymbol{X}_{i},$

where the sum is not necessarily direct. A subspace correction method solves
for an approximation of the error in each subspace, and combines them
(additively or multiplicatively) to form an updated guess for the solution. In
the classical Jacobi and Gauss–Seidel iterations, the decomposition, Eq. 52,
is given by $\\{\boldsymbol{X}_{i}\\}=\\{\boldsymbol{\phi}_{i}\\}$ where
$\boldsymbol{\phi}_{i}$, $i=1,\dots,n_{\boldsymbol{u}}$, are the velocity
basis functions. The difference between Jacobi and Gauss-Seidel is whether the
updates are applied additively (Jacobi) or multiplicatively (Gauss–Seidel).

A sufficient condition for the subspace correction method induced by the
decomposition Eq. 52 to be robust in $\gamma_{d}$ for a symmetric positive-
definite matrix, is that the decomposition captures the kernel
$\mathcal{N}_{h}$ in the following sense [14, 62, 31, 40]:

(53)
$\displaystyle\mathcal{N}_{h}=\sum_{i}\boldsymbol{X}_{i}\cap\mathcal{N}_{h}.$

In other words, the decomposition must be sufficiently rich so that any
divergence-free velocity can be written as a combination of divergence-free
functions from the subspaces $\boldsymbol{X}_{i}$111In addition, the
decomposition must be stable, but we do not elaborate here.. Jacobi fails this
criterion, as each $\boldsymbol{\phi}_{i}$ is not divergence-free. A
decomposition satisfying Eq. 53 for the BDM-DG discretization was developed by
Hong et al. [40, Sec. 4.5], where the decomposition is the so-called _star_
patch around every vertex of the mesh. This decomposition is visualized in
Fig. 3 for a $\mathrm{BDM}_{1}$ discretization in two dimensions. The same
topological decomposition extends to higher orders and three dimensions. Since
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
is not guaranteed to be positive-definite, the theory does not guarantee
robust convergence. Nevertheless, we find that a small number of FGMRES
iterations preconditioned with the vertex-star patch iteration is very
effective as a smoother, as reported in [31] and subsequent works.

Figure 3: The patch of degrees of freedom (black dots inside the blue patch)
around a vertex (red diamond) used in the multigrid relaxation for a
$\mathrm{BDM}_{1}$ discretization in two dimensions. Each vertex-star patch
contains 12 degrees of freedom in two dimensions with this lowest-order
element pair.

#### Injecting the active set

A complication arises in the representation of
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
on the coarser levels. By first ignoring the BM active set, we note that
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
can be assembled by injecting the current finite element iterates,
$\boldsymbol{u}_{h}$ and $\rho_{h}$, to the relevant level on the mesh
hierarchy, assembling the submatrices $\boldsymbol{\mathsf{A}}$,
$\boldsymbol{\mathsf{D}}$, $\boldsymbol{\mathsf{D}}^{\top}$ and
$\boldsymbol{\mathsf{C}}_{\mu}$, applying the Dirichlet boundary conditions of
the injected velocity to the relevant rows and columns of
$\boldsymbol{\mathsf{A}}$, $\boldsymbol{\mathsf{D}}$ and
$\boldsymbol{\mathsf{D}}^{\top}$, and subtracting the triple matrix product
$\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
from $\boldsymbol{\mathsf{A}}$. The triple matrix product is sparse and cheap
to compute as $\boldsymbol{\mathsf{C}}_{\mu}$ is diagonal on all levels.
However, we found that an accurate representation of the active set on the
coarser levels is essential for the convergence of the multigrid scheme. This
is similar to other experiences reported in the literature [23, 41]. Hence,
the difficulty lies in defining the active set on the coarser levels. An
obvious choice is to use the definitions of $\mathcal{A}$ and $\mathcal{I}$ in
Algorithm 1 defined via the injected material distribution iterate. However,
in numerical experiments we found that this choice resulted in poor iteration
counts.

Consider a two-grid method with the fine-level triangulation
$\mathcal{T}_{h}$, $h=H/2$, obtained by a uniform refinement of the simplices
in coarse-level triangulation $\mathcal{T}_{H}$. As the material distribution
is discretized with $\mathrm{DG}_{0}$ elements, each degree of freedom $i$
associated with the fine-level material distribution iterate can be associated
with an element $K_{h}\in\mathcal{T}_{h}$ in the fine level and analogously
with the degrees of freedom of the coarse-level material distribution iterate
with elements in the coarse level. We say that a fine-level element
$K_{h}\in\mathcal{T}_{h}$ is in the active set $\mathcal{A}_{h}$ (written as
$K_{h}\in\mathcal{A}_{h}$) if the degree of freedom associated with $K_{h}$ is
in the active set $\mathcal{A}_{h}$. This definition naturally extends to the
coarse-level elements and active sets.

We now utilize an idea inspired by the work of Hoppe [41] and Engel and
Griebel [23] to define the coarse-level active sets. A coarse-level element,
$K_{H}\in\mathcal{T}_{H}$ containing the parent fine-level elements
$K_{h,1},\dots,K_{h,s}\in\mathcal{T}_{h}$ is defined to be in the coarse-level
active set $\mathcal{A}_{H}$ if

(54) $\displaystyle|\\{K_{h,j}\in\mathcal{A}_{h}:j=1,\dots,s\\}|\geq m,$

where $m\in[1,s]$ and $s=4$ in two dimensions and $s=8$ in three dimensions.
In other words, the coarse-level element is in the coarse-level active set if
it contains $m$ or more fine-level parent elements that are in the fine-level
active set. By starting at the finest-level active set that is defined by
Algorithm 1, we recursively define all the active sets in mesh hierarchy via
Eq. 54. Experiments revealed that a good choice for $m$ is $m=s/2$, i.e. a
coarse-level element is active if at least half of its parent fine-level
elements are active. A summary of the multigrid strategy is given in Fig. 4.

Figure 4: The multigrid solver strategy of Section 5.2 to solve
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$.

######  Remark 5.4.

The choice of Eq. 54 is more generous than utilizing the definition of the
fine-level active set directly with the injected material distribution
iterate. In particular, some coarse cells that are “borderline” between the
active set and inactive set are placed in the active set by Eq. 54 but in the
inactive set when defining the active set using the injected material
distribution. Numerically, we see that the iteration counts suffer if the
criteria for a coarse cell to be in the coarse-level active set are too
strict.

##  6 Numerical results

All examples in this chapter were implemented with the finite element software
Firedrake [57]. Block preconditioning and Krylov subspace methods were
implemented using Firedrake [57] and PETSc [9], and sparse LU factorizations
were performed with MUMPS [6]. Vertex-star patch relaxation is implemented via
the PCPATCH functionality [29] recently introduced to PETSc. The meshes were
created in Firedrake or Gmsh [38]. The uniqueness of the pressure was enforced
by orthogonalizing against the nullspace of constants in the Krylov method.
The coarsest-level correction in the multigrid scheme of Section 5.2 is
computed via an LU factorization. The BM updates are scaled with a (damped)
$l^{2}$-minimizing linesearch [18, Alg. 2] and we do not use a prediction step
in any examples. Wherever (F)GMRES is used, it is not restarted.

###  6.1 Double-pipe

The first example is the two-dimensional double-pipe problem first introduced
by Borrvall and Petersson [15, Sec. 4.5]. The double-pipe problem is posed on
a rectangular domain $\Omega=(0,3/2)\times(0,1)$ with two inlets and two
outlets fixed by the following Dirichlet boundary condition

(55)
$\displaystyle\boldsymbol{g}(x,y)=\begin{cases}\left(1-144(y-3/4)^{2},0\right)^{\top}&\text{if}\;\;2/3\leq
y\leq 5/6,x=0\;\text{or}\;3/2,\\\
\left(1-144(y-1/4)^{2},0\right)^{\top}&\text{if}\;\;1/6\leq y\leq
1/3,x=0\;\text{or}\;3/2,\\\ (0,0)^{\top}&\text{elsewhere
on}\;\partial\Omega.\end{cases}$

We choose a volume fraction of $\gamma=1/3$ and the inverse permeability
$\alpha$ is given in Eq. 2, with $\bar{\alpha}=2.5\times 10^{4}$ and $q=1/10$.
The problem supports two minima: a local minimum of two straight channels from
each inlet to its opposite outlet, and a global minimum in the shape of a
double-ended wrench. These are depicted in Fig. 5.

Figure 5: The material distribution of the straight-channel (left) and double-
ended wrench (right) solutions of the double-pipe optimization problem.

The two strategies we utilize for solving the linear systems are the
following:

  1. (aL1)

The nested block preconditioning approach of Section 5.1 on Eq. 28 with
$\gamma_{d}=10^{4}$, and an LU factorization for the augmented momentum block
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$;

  2. (aL2)

The nested block preconditioning approach of Section 5.1 on Eq. 28 with
$\gamma_{d}=10^{4}$, and the geometric multigrid method of Section 5.2 to
approximate the action of the inverse of
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$.
We fix the relaxation to 5 FGMRES iterations preconditioned with a vertex-star
patch iteration and a full multigrid cycle is used.

We opt for a first-order Brezzi–Douglas–Marini
$\mathrm{BDM}_{1}\times\mathrm{DG}_{0}$ mixed finite element discretization
for the velocity-pressure pair, with interior penalty parameter $\sigma=10$,
and a $\mathrm{DG}_{0}$ discretization for the material distribution. This
choice of discretization makes both $\boldsymbol{\mathsf{M}}_{p}$ and
$\boldsymbol{\mathsf{C}}_{\mu}$ diagonal. For all mesh sizes, we initialize
the deflated barrier method at $\mu_{0}=105$ and perform deflation immediately
to find the second branch. The first and second branches converge to the
straight channels and double-ended wrench solutions, respectively, as $\mu\to
0$. The nonlinear solves are terminated with an absolute tolerance of
$10^{-5}$. The outer FGMRES solver’s absolute and relative tolerances are both
set to $10^{-7}$. In Tables 1 and 2, we list the iteration counts for the
strategies (aL1) and (aL2) on meshes with decreasing mesh sizes. In the (aL2)
strategy, the augmented block solve is approximated to an absolute tolerance
of $10^{-8}$ or a relative tolerance of $10^{-9}$.

|  | (aL1)
---|---|---
$h$ | Dofs | BM | OK
0.0361 | 25,201 | 274 | 461 (1.68)
0.0180 | 100,401 | 626 | 1180 (1.88)
0.0090 | 400,801 | 733 | 1283 (1.75)
0.0045 | 1,601,601 | 809 | 1607 (1.99)
Table 1: The total cumulative number of iterations to compute both minimizers of the double-pipe problem over all the subproblems with the (aL1) preconditioner. BM stands for the number of Benson–Munson iterations and OK stands for the number of outer Krylov FGMRES iterations. The numbers in brackets in the OK column are the number of average Krylov iterations per BM iteration. |  |  (aL2) 2-grid
---|---|---
$h$ | Dofs | BM | OK | IK
0.0180 | 100,401 | 627 | 1270 (2.03) | 100,109 (13.14)
0.0090 | 400,801 | 735 | 1360 (1.85) | 86,745 (10.63)
Table 2: The total cumulative number of iterations to compute both minimizers
of the double-pipe problem over all the subproblems with the (aL2)
preconditioner. BM stands for the number of Benson–Munson iterations, OK
stands for the number of outer Krylov FGMRES iterations, and IK is the number
of inner Krylov FGMRES iterations preconditioned with the geometric multigrid
method of Section 5.2 to solve linear systems involving
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$.
The numbers in brackets in the OK and IK columns are the number of average
Krylov iterations per BM iteration and per
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
solve, respectively.

We see that the Krylov iterations per BM iteration are robust to the mesh size
for both preconditioning strategies. The preconditioning strategy with an LU
factorization for the augmented momentum block (aL1) averages to under 2
preconditioned FGMRES iterations per BM iteration. Similarly with (aL2), where
the augmented block solve is approximated with FGMRES preconditioned with a
2-grid multigrid cycle, the outer FGMRES iterations remain under 2
preconditioned FGMRES iterations per BM iteration on average. Moreover, the
average inner FGMRES based on the kernel-preserving multigrid scheme
iterations remain under 13 iterations per augmented momentum block solve over
all mesh sizes. In particular, the average inner FGMRES iterations decreases
on the fine mesh relative to the coarse mesh. We note that, unlike the
conforming discretizations in [51], the number of BM iterations slowly
increases with decreasing mesh size. This may be due to the discretization or
could be related to the fact that the linear systems are not being solved
exactly. We do not report the timings of the solves in this example. In two
dimensions a direct solve, using MUMPS with $\gamma_{d}=0$, of the full BM
system Eq. 28 is faster than the block preconditioning and multigrid strategy
proposed here. However, we note that in three dimensions, a direct solve often
fails due to the ill-conditioning of the system. Moreover, assembling fine-
grid matrices has large memory requirements. As the (aL2) strategy only
requires an assembly of the coarsest-grid augmented block, this allows for
finer mesh linear solves. This is further discussed in Section 6.3.

A natural question is the sensitivity of the solver with respect to the
augmented Lagrangian parameter $\gamma_{d}$. We wish to explore the following:

  1. (H1)

The effect of roundoff error as $\gamma_{d}$ increases;

  2. (H2)

The improvement of the approximation of
$-\gamma_{d}\boldsymbol{\mathsf{M}}_{p}^{-1}$ for
$\boldsymbol{\mathsf{S}}_{2,\gamma_{d}}^{-1}$;

  3. (H3)

The robustness of the vertex-star patch multigrid cycle for increasing
$\gamma_{d}$.

An unfortunate side-effect of increasing $\gamma_{d}$ is the expectation of
high roundoff error. This would manifest in the nonlinear solver stagnating at
higher residual norm values than desired. However, as $\gamma_{d}$ increases,
we expect that the Schur complement approximation by the pressure mass matrix
improves as shown in Proposition 5.2. Hence, the number of outer Krylov FGMRES
iterations per BM iteration should decrease as $\gamma_{d}\to\infty$. Finally,
the multigrid cycle should be robust to the value of $\gamma_{d}$. Hence, we
expect that the number of inner Krylov FGMRES iterations preconditioned by the
multigrid cycle to approximate the inverse of the augmented momentum block
should stay roughly constant with increasing $\gamma_{d}$. In Table 3 we test
these hypotheses.

We investigate (H1) and (H2) via the (aL1) solver (with varying $\gamma_{d}$).
We fix a mesh with mesh size $h=0.0361$ and choose the discretization used
previously. We consider the deflated barrier method subproblem at $\mu=105$
and $\mu=1$ for the first branch. The nonlinear solves are terminated when the
decrease in the residual norm stagnates. The outer FGMRES solver’s absolute
and relative tolerances are set to $10^{-15}$ and $10^{-5}$, respectively. The
middle four columns confirm our first two hypotheses. At $\gamma_{d}=10^{7}$
we achieve the optimal number of outer FGMRES iterations per BM iteration. The
approximation of the Schur complement with the pressure mass matrix in near-
perfect. However, we note that for each increase in the order of magnitude
after $\gamma_{d}=10^{2}$, the roundoff error causes the nonlinear solver to
stagnate at an order of magnitude higher. We believe a good compromising
choice is $\gamma_{d}=10^{4}$.

To investigate (H3) we utilize the (aL2) solver (with varying $\gamma_{d}$).
We fix a mesh with mesh size $h=0.0361$ for the coarse level and uniformly
refine the mesh once for the fine level. We consider the first outer Krylov
FGMRES iteration at $\mu=105$ and tabulate the average number of multigrid
preconditioned inner Krylov FGMRES iterations per augmented momentum block
solve required. The augmented momentum block solve is approximated to an
absolute tolerance of $10^{-8}$ or a relative tolerance of $10^{-9}$. We see
that with increasing $\gamma_{d}$ the iterations stay roughly constant. This
confirms the cycles are robust in $\gamma_{d}$.

| (aL1) | (aL2)
---|---|---
| $\mu=105$ | $\mu=1$ | $\mu=105$
$\gamma_{d}$ | avg. OK | res. norm | avg. OK | res. norm | avg. IK
$10^{-5}$ | 157.2 | $9.90\times 10^{-12}$ | 66.0 | $6.40\times 10^{-12}$ | 10.0
$10^{-1}$ | 133.5 | $9.93\times 10^{-12}$ | 93.9 | $6.73\times 10^{-12}$ | 11.3
$10^{0}$ | 77.8 | $1.09\times 10^{-11}$ | 40.6 | $6.92\times 10^{-12}$ | 12.8
$10^{1}$ | 30.0 | $1.15\times 10^{-11}$ | 16.3 | $8.99\times 10^{-12}$ | 10.0
$10^{2}$ | 8.23 | $3.92\times 10^{-11}$ | 4.44 | $4.73\times 10^{-11}$ | 9.50
$10^{3}$ | 3.31 | $3.75\times 10^{-10}$ | 3.12 | $4.67\times 10^{-10}$ | 9.17
$10^{4}$ | 1.85 | $3.71\times 10^{-9}$ | 2.20 | $4.57\times 10^{-9}$ | 9.50
$10^{5}$ | 1.33 | $3.70\times 10^{-8}$ | 1.67 | $4.68\times 10^{-8}$ | 9.50
$10^{6}$ | 1.25 | $3.77\times 10^{-7}$ | 1.27 | $4.57\times 10^{-7}$ | 9.50
$10^{7}$ | 1.14 | $3.73\times 10^{-6}$ | 1.20 | $4.59\times 10^{-6}$ | 9.50
$10^{8}$ | 1.18 | $3.63\times 10^{-5}$ | 1.50 | $4.92\times 10^{-5}$ | 9.50
Table 3: The sensitivity of the (aL1) and (aL2) solver strategies to the value
of $\gamma_{d}$. For the (aL1) strategy, the columns labelled avg. OK are the
average number of outer Krylov FGMRES iterations per BM iteration to solve the
deflated barrier method subproblem at $\mu=105$ and $\mu=1$ for $h=0.0361$.
The columns labelled res. norm hold the values of the smallest residual norm
before the nonlinear solver stagnates due to roundoff errors. For the (aL2)
strategy the coarse level has a mesh size of $h=0.0361$ and we perform one
uniform refinement. The column labelled avg. IK is the average number of
multigrid preconditioned inner Krylov FGMRES iterations per augmented momentum
block solve as seen in the first outer Krylov FGMRES iteration at subproblem
$\mu=105$.

###  6.2 3D cross-channel

The first three-dimensional example we consider is the cross-channel problem
as found in Sá et al. [58, Sec. 7.5]. The domain is the unit cube,
$\Omega=(0,1)^{3}$, with two circular inlets and two circular outlets that are
arranged in a cross pattern as visualized in Fig. 6. The volume fraction is
given by $\gamma=1/10$ and we use Eq. 2 as our choice of $\alpha$, with
$\bar{\alpha}=2.5\times 10^{4}$ and $q=1/10$. The Dirichlet boundary condition
on $\boldsymbol{u}$ is

(56)
$\displaystyle\boldsymbol{g}(x,y,z)=\left(1-12\pi((y-a)^{2}+(z-b)^{2}),0,0\right)^{\top},$

if $12\pi((y-a)^{2}+(z-b)^{2})\leq 1$ and $x=0$ with $a=1/2$,
$b\in\\{1/4,3/4\\}$ or $x=1$ with $a\in\\{1/4,3/4\\}$, $b=1/2$, and
$\boldsymbol{g}(x,y,z)=(0,0,0)^{\top}$ elsewhere on $\partial\Omega$.

Figure 6: Setup of the 3D cross-channel problem. This problem features a unit
cube domain with two inlets and two outlets arranged in a cross pattern.

We apply the same first-order BDM discretization, with interior penalty
penalization parameter $\sigma=10$, and run the deflated barrier method twice.
The first pass is on a $20\times 20\times 20$ mesh resulting in 391,201
degrees of freedom. The augmented Lagrangian parameter is chosen to be
$\gamma_{d}=10^{6}$. Due to the nested block preconditioning, the action of
the inverse of the augmented momentum block
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
must be applied six times per outer FGMRES iteration. On this relatively
coarse mesh, it is cheaper to factorize
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
with MUMPS at the start of each BM iteration and reuse the factorization,
rather than iteratively solve
$\boldsymbol{\mathsf{A}}_{\gamma_{d}}-\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}_{\mu}^{-1}\boldsymbol{\mathsf{D}}^{\top}$
with multigrid each time an inverse action is required, i.e. we use the
augmented Lagrangian preconditioner (aL1) for the linear systems. We
initialize the barrier parameter at $\mu_{0}=100$. The nonlinear solves are
terminated with an absolute tolerance of $10^{-5}$ and the outer FGMRES
solver’s absolute tolerance is set to $5\times 10^{-7}$ (all linear solves
converged to the absolute tolerance). A second branch of solutions is found at
$\mu=38.74$ and a third branch at $\mu=34.87$.

These coarse-mesh solutions at $\mu=0$ are interpolated onto a finer $40\times
40\times 40$ mesh resulting in 3,100,801 degrees of freedom. The deflated
barrier method is then reinitialized at $\mu_{0}=10^{-6}$; we found that the
BM solver often diverges if initialized at $\mu_{0}=0$. On this finer mesh, we
again apply the augmented Lagrangian preconditioner ($\gamma_{d}=10^{5}$). The
nonlinear solves are terminated with an absolute tolerance of $10^{-4}$ if
$\mu>0$ and $10^{-6}$ if $\mu=0$. Moreover, the outer FGMRES solver’s absolute
tolerance is set to $10^{-7}$ if $\mu>0$ and $10^{-9}$ if $\mu=0$ and its
relative tolerance is set to $10^{-7}$. Now a direct solve of the augmented
momentum block is prohibitive and we switch to the (aL2) strategy where each
approximate inverse of the augmented momentum block is solved to an absolute
or relative tolerance of $10^{-8}$ or $10^{-9}$, respectively, with the
(2-grid) multigrid scheme of Section 5.2. For the relaxation on the fine
level, we use 5 FGMRES iterations preconditioned with the vertex-star patch
relaxation.

The resulting three solutions are shown in Fig. 7. Two of these are symmetric
straight channel solutions where the inlets swap which outlet they exit from.
Their symmetry results in similar costs. A third global minimizer comes in the
form of a merged channel solution; the two channels briefly merge in the
middle of the box domain before splitting to exit via the two outlets.

The iteration counts for the initial search on a $20\times 20\times 20$ mesh
and the refinement on a $40\times 40\times 40$ mesh are given in Table 4. We
see that our preconditioner is effective in both cases. When using a direct
solve of the augmented momentum block, we average slightly more than one outer
FGMRES iterations per BM iteration. Similarly, the tolerances for the
augmented momentum block solve on the fine mesh are strict enough so that the
outer FGMRES iterations average between 2.24 and 2.55. Moreover, each
augmented momentum block solve requires an average in the range of 15.77–16.77
multigrid preconditioned FGMRES iterations to reach the prescribed tolerances.

(a) Branch 0

(b) Branch 1

(c) Branch 2

Figure 7: The material distribution of the solutions discovered by the deflated barrier method to the 3D cross-channel optimization problem discretized with 3,100,801 degrees of freedom. The power dissipation values are $J_{h}=14.51,14.62$, and $13.08$ for branches 0, 1, and 2, respectively. | Coarse mesh, $\gamma_{d}=10^{6}$ | Fine mesh, $\gamma_{d}=10^{5}$
---|---|---
Branch | BM | OK | BM | OK | IK
0 | 828 | 829 (1.00) | 49 | 110 (2.24) | 10,405 (15.77)
1 | 444 | 445 (1.00) | 52 | 121 (2.33) | 12,176 (16.77)
2 | 409 | 422 (1.03) | 53 | 135 (2.55) | 13,085 (16.15)
Table 4: Cumulative number of BM iterations, outer FGMRES iterations (OK), and
for the fine mesh, inner FGMRES iterations preconditioned with the multigrid
scheme of Section 5.2 (IK) for the 3D cross-channel problem. The bracketed
numbers in the OK and IK columns are the average number of outer FGMRES
iterations per BM iteration and average number of inner FGMRES iterations per
augmented momentum block solve, respectively. The barrier parameter is
initialized at $\mu_{0}=100$ on the coarse mesh and $\mu_{0}=10^{-6}$ on the
fine mesh.

###  6.3 3D five-holes quadruple-pipe

In [51, Sec. 4.4], it was observed that introducing holes in a rectangular
domain caused a significant increase in the number of solutions. We now extend
this idea to three dimensions and introduce the generalization of the five-
holes double-pipe problem. This problem features a box domain
$\Omega=(0,3/2)\times(0,1)\times(0,1)$ with five internal holes in the shape
of cubes, of edge length $1/10$, with centres at $(3/4,1/4,1/4)$,
$(3/4,1/4,3/4)$, $(3/4,3/4,1/4)$, $(3/4,3/4,3/4)$, and $(3/4,1/2,1/2)$. There
are four inlets and four outlets. The circular inlets of radius
$1/\sqrt{12\pi}$ are positioned on the face $x=0$ with the centres
$(y,z)=(1/4,1/4)$, $(1/4,3/4)$, $(3/4,1/4)$, and $(3/4,3/4)$. The circular
outlets of the same radius are positioned on the face $x=3/2$ with the same
centres. The domain setup is depicted in Fig. 8. We impose a parabolic
Dirichlet boundary condition on the inlets and outlets and a zero Dirichlet
boundary condition elsewhere on the boundary (including the boundary of the
five internal holes), i.e. the Dirichlet boundary condition is given by:

(57)
$\displaystyle\boldsymbol{g}(x,y,z)=\left(1-12\pi((y-a)^{2}+(z-b)^{2}),0,0\right)^{\top},$

if $12\pi((y-a)^{2}+(z-b)^{2})\leq 1$, where
$a,b\in\\{1/4,3/4\\},\;x=0\;\text{or}\;3/2$ and

(58) $\displaystyle\boldsymbol{g}(x,y,z)=(0,0,0)^{\top},$

elsewhere on $\partial\Omega$, including the boundaries of the five internal
holes. We choose a volume fraction of $\gamma=1/5$ and the inverse
permeability term $\alpha$ is given by Eq. 2, with $\bar{\alpha}=2.5\times
10^{4}$ and $q=1/10$.

Figure 8: Setup of the 3D five-holes quadruple-pipe problem. This problem
features 4 inlets and 4 outlets. The domain is a box,
$\Omega=(0,3/2)\times(0,1)\times(0,1)$, with five internal holes in the shape
of cubes with edge length $1/10$ that are centred at $(3/4,1/4,1/4)$,
$(3/4,1/4,3/4)$, $(3/4,3/4,1/4)$, $(3/4,3/4,3/4)$, and $(3/4,1/2,1/2)$.

The domain is nonconvex and has a complicated geometry due to the cuboidal
internal holes, which we do not mesh. We build the coarse level mesh with
Gmsh; the finer levels are built by uniform refinement of the coarse mesh. We
choose a first-order BDM discretization for the velocity-pressure pair, with
interior penalty penalization parameter $\sigma=10^{3}$, and run the deflated
barrier method twice. The larger choice for $\sigma$ is required to
sufficiently enforce the boundary conditions in the tangential directions. The
first run is on a (relatively) coarse mesh with 30,848 elements which results
in 256,745 degrees of freedom. The barrier parameter is initialized at
$\mu_{0}=200$ and we use the augmented Lagrangian preconditioner (aL1) for the
linear systems, with an augmented Lagrangian parameter value of
$\gamma_{d}=10^{5}$. The nonlinear solves are terminated with an absolute
tolerance of $10^{-4}$. The outer FGMRES solver’s absolute tolerance is set to
$5\times 10^{-7}$ (all linear solves converged to the absolute tolerance).

In total we find 14 coarse-grid solutions. Branches 1 and 2 are found at
$\mu=53.81$. Branches 3, 4, 5, and 6 are found at $\mu=11.39$. Branch 7 is
found at $\mu=10.25$. Branch 8 is found at $\mu=9.23$. Branch 9 is found at
$\mu=6.73$, branch 10 at $\mu=6.73$, branches 11 and 12 are found at
$\mu=6.05$, and, finally, branch 13 is found at $\mu=4.41$.

We interpolate the coarse-level solutions onto the first refinement, which
results in 2,014,113 degrees of freedom. The deflated barrier method is
initialized at $\mu_{0}=10^{-5}$, using the coarse-level solutions as initial
guesses. We damp the $l^{2}$-minimizing linesearch in the nonlinear solver by
a factor of 1/2. The nonlinear solves are terminated with an absolute
tolerance of $10^{-4}$ if $\mu>0$ and $10^{-6}$ if $\mu=0$. The outer FGMRES
solver’s absolute tolerance is set to $10^{-6}$ if $\mu>0$ and $10^{-7}$ if
$\mu=0$, and its relative tolerance is set to $10^{-5}$. We apply an augmented
Lagrangian 2-grid multigrid preconditioner (aL2) to the linear systems with 5
FGMRES iterations for the relaxation of the fine level and
$\gamma_{d}=10^{5}$. The augmented momentum block is solved to an absolute or
relative tolerance of $10^{-4}$ or $10^{-7}$, respectively. Of the original 14
solutions, the nonlinear solver successfully converges to 7 fine-grid
solutions.

The resulting iteration counts are given in Table 5, the resulting fine mesh
solutions are shown in Fig. 9, and their cross sections are shown in Fig. 10.
A further mesh refinement of branch 0 is shown in Fig. 11 as computed using
the augmented Lagrangian (aL2) with a 3-grid multigrid cycle. As expected, the
five holes obstruct the channels and prevent a large channel passing through
the centre. The best solutions found are branches 0, 1, and 11 where the
channels form one large channel and either move to the left, upwards or
downwards to avoid the middle internal hole. As in the five-holes double-pipe
example, there are remaining solutions that we have not yet computed, by
symmetry.

In Table 6 we give timings for the solvers as run on a machine with 512 GB of
RAM and 32 CPUs Intel(R) Xeon(R) CPU E5-4627 v2 @ 3.30GHz. We note that a
direct LU factorization of the full BM system Eq. 28 fails to converge due to
the ill-conditioning of the system (even with $\gamma_{d}=0$). Thus the
preconditioning strategies are necessary to find the BM updates. On the coarse
mesh, the (aL1) strategy is faster than the (aL2) strategy. Moreover, we show
that the average time taken per BM iteration drops from 277 seconds with 1 CPU
to 35 seconds when using 32 CPUs. On the fine mesh, the assembly of the full
fine-grid augmented block requires more memory than what is available on the
workstation (512 GB). Hence only the (aL2) strategy succeeded in grid-
sequencing the solutions on this machine.

| Coarse mesh | Fine mesh
---|---|---
Branch | BM | OK | BM | OK | IK
0 | 438 | 438 (1) | 41 | 109 (2.66) | 15,828 (24.2)
1 | 338 | 338 (1) | 36 | 97 (2.69) | 13,736 (23.6)
2 | 396 | 396 (1) | 49 | 130 (2.65) | 18,999 (24.4)
3 | 254 | 254 (1) | 46 | 114 (2.48) | 16,464 (24.1)
7 | 221 | 221 (1) | 53 | 132 (2.49) | 20,061 (25.3)
11 | 183 | 183 (1) | 64 | 180 (2.81) | 27,738 (25.7)
13 | 187 | 187 (1) | 60 | 169 (2.82) | 28,542 (28.1)
Table 5: Cumulative number of BM iterations, outer FGMRES iterations (OK), and for the fine mesh, inner FGMRES iterations preconditioned with the multigrid scheme of Section 5.2 (IK) for the 3D five-holes quadruple-pipe problem. The bracketed numbers in the OK and IK columns are the average number of outer FGMRES iterations per BM iteration and average number of inner FGMRES iterations per augmented momentum block solve, respectively. The barrier parameter is initialized at $\mu_{0}=200$ on the coarse mesh and $\mu_{0}=10^{-5}$ on the fine mesh. Strategy | BM | Time taken (s) | Avg. time taken (s)
---|---|---|---
Coarse mesh
(aL1) 1 CPU | 6 | 1664.02 | 277.33
(aL1) 4 CPUs | 44 | 4425.71 | 100.58
(aL1) 8 CPUs | 37 | 2310.22 | 62.44
(aL1) 16 CPUs | 42 | 1844.64 | 43.92
(aL1) 32 CPUs | 6 | 208.64 | 34.77
(aL2) 32 CPUs | 7 | 2150.24 | 307.18
Fine mesh
(aL2) 32 CPUs | 41 | 53,356.56 | 1301.38
Table 6: Time measurements of the linear solves of the BM system. On the
coarse mesh we measure the number of BM iterations and overall time taken to
find the first branch at $\mu=200$. On the fine mesh we measure the number of
BM iterations and overall time taken to grid-sequence the first solution
starting at $\mu=10^{-5}$.

(a) Branch 0

(b) Branch 1

(c) Branch 2

(d) Branch 3

(e) Branch 7

(f) Branch 11

(g) Branch 13

Figure 9: The material distribution of seven solutions to the 3D fives-holes
quadruple-pipe optimization problem discretized with 2,014,113 degrees of
freedom as discovered by the deflated barrier method and grid-sequenced once.
The resulting power dissipation values for branches 0, 1, 2, 3, 7, 11, and 13
are $J_{h}=55.02,54.73,62.27,62.22,59.56,55.28$, and $62.78$, respectively.

(a) Branch 0

(b) Branch 1

(c) Branch 2

(d) Branch 3

(e) Branch 7

(f) Branch 11

(g) Branch 13

Figure 10: The crinkled cross sections at $x=3/4$ for the discovered solutions
of the 3D five-holes quadruple-pipe. The grey regions are part of the five
cuboid holes in the box domain. The material distribution has a value of one
in the red regions and zero in the blue regions, with intermediate values for
the intermediate coloured regions. Figure 11: The material distribution of
branch 0 after grid-sequencing to a second uniform mesh refinement resulting
in 15,953,537 degrees of freedom ($J_{h}=39.11$). The augmented Lagrangian
3-grid multigrid preconditioner (aL2) was used.

Code availability: For reproducibility, the implementation of the deflated
barrier method used in this work, as well as scripts to generate the
solutions, can be found at https://github.com/ioannisPApapadopoulos/fir3dab/.
The version of the software used in this paper is archived on Zenodo [1].

##  7 Conclusions

In this work we extended the application of the deflated barrier method to
discovering multiple three-dimensional solutions of the fluid topology
optimization model of Borrvall and Petersson [15]. This was achieved by
developing preconditioners for the linear systems that arise in the deflated
barrier method. The preconditioning strategy reduces the discretized $4\times
4$-block matrix solve to an outer FGMRES method and the solve of one diagonal
matrix, one block-diagonal matrix (that can be factorized once or quickly
solved with multigrid), and an augmented momentum block. Moreover, we develop
a geometric multigrid cycle for the augmented momentum block that consists of
a specialized relaxation method and a characterization of the active set on
coarser levels. We found that the preconditioner is robust to the mesh size
and successfully computed three solutions of a 3D cross-channel problem and
seven solutions of a 3D five-holes quadruple-pipe problem.

## References

  * [1] Software used in ‘Preconditioners for computing multiple solutions in three-dimensional fluid topology optimization’, 2022, https://doi.org/10.5281/zenodo.5899634.
  * [2] N. Aage and B. S. Lazarov, Parallel framework for topology optimization using the method of moving asymptotes, Structural and Multidisciplinary Optimization, 47 (2013), pp. 493–505, https://doi.org/10.1007/s00158-012-0869-2.
  * [3] N. Aage, T. H. Poulsen, A. Gersborg-Hansen, and O. Sigmund, Topology optimization of large scale Stokes flow problems, Structural and Multidisciplinary Optimization, 35 (2007), pp. 175–180, https://doi.org/10.1007/s00158-007-0128-0.
  * [4] J. Alexandersen and C. S. Andreasen, A review of topology optimisation for fluid-based problems, Fluids, 5 (2020), p. 29, https://doi.org/10.3390/fluids5010029.
  * [5] F. P. Ali Beik and M. Benzi, Iterative methods for double saddle point systems, SIAM Journal on Matrix Analysis and Applications, 39 (2018), pp. 902–921, https://doi.org/10.1137/17M1121226.
  * [6] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications, (2001), https://doi.org/10.1137/S0895479899358194.
  * [7] D. N. Arnold, R. S. Falk, and R. Winther, Multigrid in $H(\mathrm{div})$ and $H(\mathrm{curl})$, Numerische Mathematik, 85 (2000), pp. 197–217, https://doi.org/10.1007/PL00005386.
  * [8] F. R. A. Aznaran, R. C. Kirby, and P. E. Farrell, Transformations for Piola-mapped elements, SMAI Journal of Computational Mathematics, (2022).  In press.
  * [9] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkout, W. Gropp, R. Tran Mills, T. Munson, K. Rupp, P. Sana, B. Smith, S. Zampini, H. Zhang, and H. Zhang, PETSc users manual, Tech. Report ANL-95/11 - Revision 3.11, Argonne National Laboratory, 2019, http://www.mcs.anl.gov/petsc.
  * [10] R. Behrou, R. Ranjan, and J. K. Guest, Adaptive topology optimization for incompressible laminar flow problems with mass flow constraints, Computer Methods in Applied Mechanics and Engineering, 346 (2019), pp. 612–641, https://doi.org/10.1016/j.cma.2018.11.037.
  * [11] M. P. Bendsøe and O. Sigmund, Topology Optimization, Springer Berlin Heidelberg, Berlin, Heidelberg, 2004, https://doi.org/10.1007/978-3-662-05086-6.
  * [12] S. J. Benson and T. S. Munson, Flexible complementarity solvers for large-scale applications, Optimization Methods and Software, 21 (2003), pp. 155–168, https://doi.org/10.1080/10556780500065382.
  * [13] M. Benzi, G. H. Golub, and J. Liesen, Numerical solution of saddle point problems, Acta Numerica, 14 (2005), pp. 1–137, https://doi.org/10.1017/S0962492904000212.
  * [14] M. Benzi and M. A. Olshanskii, An augmented Lagrangian-based approach to the Oseen problem, SIAM Journal on Scientific Computing, 28 (2006), pp. 2095–2113, https://doi.org/10.1137/050646421.
  * [15] T. Borrvall and J. Petersson, Topology optimization of fluids in Stokes flow, International Journal for Numerical Methods in Fluids, 41 (2003), pp. 77–107, https://doi.org/10.1002/fld.426.
  * [16] F. Brezzi, J. Douglas, R. Durán, and M. Fortin, Mixed finite elements for second order elliptic problems in three variables, Numerische Mathematik, 51 (1987), pp. 237–250, https://doi.org/10.1007/BF01396752.
  * [17] F. Brezzi, J. Douglas, and L. D. Marini, Two families of mixed finite elements for second order elliptic problems, Numerische Mathematik, 47 (1985), pp. 217–235, https://doi.org/10.1007/BF01389710.
  * [18] P. R. Brune, M. G. Knepley, B. F. Smith, and X. Tu, Composing scalable nonlinear algebraic solvers, SIAM Review, 57 (2015), pp. 535–565, https://doi.org/10.1137/130936725.
  * [19] J. Cahouet and J.-P. Chabard, Some fast 3D finite element solvers for the generalized Stokes problem, International Journal for Numerical Methods in Fluids, 8 (1988), pp. 869–895, https://doi.org/10.1002/fld.1650080802.
  * [20] V. J. Challis and J. K. Guest, Level set topology optimization of fluids in Stokes flow, International Journal for Numerical Methods in Engineering, 79 (2009), pp. 1284–1308, https://doi.org/10.1002/nme.2616.
  * [21] Y. Deng, Z. Liu, J. Wu, and Y. Wu, Topology optimization of steady Navier–Stokes flow with body force, Computer Methods in Applied Mechanics and Engineering, 255 (2013), pp. 306–321, https://doi.org/10.1016/j.cma.2012.11.015.
  * [22] H. C. Elman, D. J. Silvester, and A. J. Wathen, Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics, Oxford University Press, 2014, https://doi.org/10.1093/acprof:oso/9780199678792.001.0001.
  * [23] M. Engel and M. Griebel, A multigrid method for constrained optimal control problems, Journal of Computational and Applied Mathematics, 235 (2011), pp. 4368–4388, https://doi.org/10.1016/j.cam.2011.04.002.
  * [24] L. C. Evans, Partial Differential Equations, American Mathematical Society, 2 ed., 2010.
  * [25] A. Evgrafov, State space Newton’s method for topology optimization, Computer Methods in Applied Mechanics and Engineering, 278 (2014), pp. 272–290, https://doi.org/10.1016/j.cma.2014.06.005.
  * [26] P. E. Farrell, Á. Birkisson, and S. W. Funke, Deflation techniques for finding distinct solutions of nonlinear partial differential equations, SIAM Journal on Scientific Computing, 37 (2015), pp. A2026–A2045, https://doi.org/10.1137/140984798.
  * [27] P. E. Farrell, M. Croci, and T. M. Surowiec, Deflation for semismooth equations, Optimization Methods and Software, (2019), pp. 1–24, https://doi.org/10.1080/10556788.2019.1613655.
  * [28] P. E. Farrell and P. A. Gazca-Orozco, An augmented Lagrangian preconditioner for implicitly constituted non-Newtonian incompressible flow, SIAM Journal on Scientific Computing, 42 (2020), pp. B1329–B1349, https://doi.org/10.1137/20M1336618.
  * [29] P. E. Farrell, M. G. Knepley, L. Mitchell, and F. Wechsung, PCPATCH: Software for the topological construction of multigrid relaxation methods, ACM Transactions on Mathematical Software (TOMS), 47 (2021), pp. 1–22, https://doi.org/10.1145/3445791.
  * [30] P. E. Farrell, L. Mitchell, L. R. Scott, and F. Wechsung, A Reynolds-robust preconditioner for the Scott-Vogelius discretization of the stationary incompressible Navier-Stokes equations, The SMAI Journal of Computational Mathematics, 7 (2021), pp. 75–96, https://doi.org/10.5802/smai-jcm.72.
  * [31] P. E. Farrell, L. Mitchell, and F. Wechsung, An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier–Stokes equations at high Reynolds number, SIAM Journal on Scientific Computing, 41 (2019), pp. A3073–A3096, https://doi.org/10.1137/18M1219370.
  * [32] P. E. Farrell, P. A. G. Orozco, and E. Süli, Finite element approximation and augmented Lagrangian preconditioning for anisothermal implicitly-constituted non-Newtonian flow., Mathematics of Computation, 91 (2022), pp. 659–697, https://doi.org/10.1090/mcom/3703.
  * [33] A. V. Fiacco and G. P. McCormick, Nonlinear Programming, Society for Industrial and Applied Mathematics, 1990, https://doi.org/10.1137/1.9781611971316.
  * [34] A. Forsgren, P. E. Gill, and M. H. Wright, Interior methods for nonlinear optimization, SIAM Review, 44 (2002), pp. 525–597, https://doi.org/10.1137/S0036144502414942.
  * [35] M. Fortin and R. Glowinski, Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, vol. 15 of Studies in Mathematics and Its Applications, Elsevier, 1983.
  * [36] K. R. Frisch, The logarithmic potential method of convex programming, tech. report, University Institute of Economics, Oslo, Norway, 1955\.
  * [37] N. R. Gauger, A. Linke, and P. W. Schroeder, On high-order pressure-robust space discretisations, their advantages for incompressible high Reynolds number generalised Beltrami flows and beyond, The SMAI Journal of Computational Mathematics, 5 (2019), pp. 89–129, https://doi.org/10.5802/smai-jcm.44.
  * [38] C. Geuzaine and J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering, 79 (2009), pp. 1309–1331, https://doi.org/10.1002/nme.2579.
  * [39] M. Hintermüller, K. Ito, and K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM Journal on Optimization, 13 (2003), pp. 865–888, https://doi.org/10.1137/S1052623401383558.
  * [40] Q. Hong, J. Kraus, J. Xu, and L. Zikatanov, A robust multigrid method for discontinuous Galerkin discretizations of Stokes and linear elasticity equations, Numerische Mathematik, 132 (2016), pp. 23–49, https://doi.org/10.1007/s00211-015-0712-y.
  * [41] R. H. Hoppe, Multigrid algorithms for variational inequalities, SIAM Journal on Numerical Analysis, 24 (1987), pp. 1046–1065, https://doi.org/10.1137/0724069.
  * [42] K. E. Jensen, Topology optimization of Stokes flow on dynamic meshes using simple optimizers, Computers & Fluids, 174 (2018), pp. 66–77, https://doi.org/10.1016/j.compfluid.2018.07.011.
  * [43] V. John, A. Linke, C. Merdon, M. Neilan, and L. G. Rebholz, On the divergence constraint in mixed finite element methods for incompressible flows, SIAM Review, 59 (2017), pp. 492–544, https://doi.org/10.1137/15M1047696.
  * [44] J. Könnö and R. Stenberg, Numerical computations with $H$(div)-finite elements for the Brinkman problem, Computational Geosciences, 16 (2012), pp. 139–158, https://doi.org/10.1007/s10596-011-9259-x.
  * [45] F. Laakmann, P. E. Farrell, and L. Mitchell, An augmented lagrangian preconditioner for the magnetohydrodynamics equations at high reynolds and coupling numbers, SIAM Journal on Scientific Computing, 44 (2022), pp. B1018–B1044, https://doi.org/10.1137/21M1416539.
  * [46] A. Linke and L. G. Rebholz, Pressure-induced locking in mixed methods for time-dependent (Navier–) Stokes equations, Journal of Computational Physics, 388 (2019), pp. 350–356, https://doi.org/10.1016/j.jcp.2019.03.010.
  * [47] K. A. Mardal, X.-C. Tai, and R. Winther, A robust finite element method for Darcy–Stokes flow, SIAM Journal on Numerical Analysis, 40 (2002), pp. 1605–1631, https://doi.org/10.1137/S0036142901383910.
  * [48] M. F. Murphy, G. H. Golub, and A. J. Wathen, A note on preconditioning for indefinite linear systems, SIAM Journal on Scientific Computing, 21 (2000), pp. 1969–1972, https://doi.org/10.1137/S1064827599355153.
  * [49] I. P. Papadopoulos, Numerical analysis of a discontinuous Galerkin method for the Borrvall–Petersson topology optimization problem, SIAM Journal on Numerical Analysis, 60 (2022), pp. 2538–2564, https://doi.org/10.1137/21M1438943.
  * [50] I. P. A. Papadopoulos, Computing multiple solutions of topology optimization problems, PhD thesis, University of Oxford, 2021, https://ora.ox.ac.uk/objects/uuid:455f5465-af68-49c7-8e11-96909cdb82a1.
  * [51] I. P. A. Papadopoulos, P. E. Farrell, and T. M. Surowiec, Computing multiple solutions of topology optimization problems, SIAM Journal on Scientific Computing, 43 (2021), pp. A1555–A1582, https://doi.org/10.1137/20M1326209.
  * [52] I. P. A. Papadopoulos and E. Süli, Numerical analysis of a topology optimization problem for Stokes flow, Journal of Computational and Applied Mathematics, 412 (2022), p. 114295, https://doi.org/10.1016/j.cam.2022.114295.
  * [53] G. Pingen, A. Evgrafov, and K. Maute, Topology optimization of flow domains using the lattice Boltzmann method, Structural and Multidisciplinary Optimization, 34 (2007), pp. 507–524, https://doi.org/10.1007/s00158-007-0105-7.
  * [54] P. Popov, Preconditioning of linear systems arising in finite element discretizations of the Brinkman equation, in Large-Scale Scientific Computing. LSSC 2011. Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 2012, pp. 381–389, https://doi.org/10.1007/978-3-642-29843-1_43.
  * [55] L. Qi, Convergence analysis of some algorithms for solving nonsmooth equations, Mathematics of Operations Research, 18 (1993), pp. 227–244, https://doi.org/10.1287/moor.18.1.227.
  * [56] L. Qi and J. Sun, A nonsmooth version of Newton’s method, Mathematical Programming, 58 (1993), pp. 353–367, https://doi.org/10.1007/BF01581275.
  * [57] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. McRae, G.-T. Bercea, G. R. Markall, and P. H. Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Transactions on Mathematical Software, 43 (2016), pp. 1–27, https://doi.org/10.1145/2998441.
  * [58] L. F. N. Sá, R. C. R. Amigo, A. A. Novotny, and E. C. N. Silva, Topological derivatives applied to fluid flow channel design optimization problems, Structural and Multidisciplinary Optimization, 54 (2016), pp. 249–264, https://doi.org/10.1007/s00158-016-1399-0.
  * [59] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM Journal on Scientific Computing, 14 (1993), pp. 461–469, https://doi.org/10.1137/0914028.
  * [60] A. Schiela, Barrier methods for optimal control problems with state constraints, SIAM Journal on Optimization, 20 (2007), pp. 1002–1031, https://doi.org/10.1137/070692789.
  * [61] A. Schiela and M. Weiser, Superlinear convergence of the control reduced interior point method for PDE constrained optimization, Computational Optimization and Applications, 39 (2008), pp. 369–393, https://doi.org/10.1007/s10589-007-9057-5.
  * [62] J. Schöberl, Multigrid methods for a parameter dependent problem in primal variables, Numerische Mathematik, 84 (1999), pp. 97–119, https://doi.org/10.1007/s002110050465.
  * [63] D. Silvester and A. Wathen, Fast iterative solution of stabilised Stokes systems Part II: Using general block preconditioners, SIAM Journal on Numerical Analysis, 31 (1994), pp. 1352–1367, https://doi.org/10.1137/0731070.
  * [64] M. Stoll and A. J. Wathen, Preconditioning for active set and projected gradient methods as semi-smooth Newton methods for PDE-constrained optimization with control constraints, (2009).
  * [65] M. Stolpe and K. Svanberg, On the trajectories of penalization methods for topology optimization, Structural and Multidisciplinary Optimization, 21 (2001), pp. 128–139, https://doi.org/10.1007/s001580050177.
  * [66] M. Ulbrich, Semismooth Newton methods for operator equations in function spaces, SIAM Journal on Optimization, 13 (2003), pp. 805–841, https://doi.org/10.1137/s1052623400371569.
  * [67] M. Ulbrich and S. Ulbrich, Primal-dual interior-point methods for PDE-constrained optimization, Mathematical Programming, 117 (2009), pp. 435–485, https://doi.org/10.1007/s10107-007-0168-7.
  * [68] A. Wathen and D. Silvester, Fast iterative solution of stabilised Stokes systems. Part I: Using simple diagonal preconditioners, SIAM Journal on Numerical Analysis, 30 (1993), pp. 630–649, https://doi.org/10.1137/0730031.
  * [69] A. J. Wathen, Preconditioning, Acta Numerica, 24 (2015), https://doi.org/10.1017/S0962492915000021.
  * [70] M. Weiser, T. Gänzler, and A. Schiela, A control reduced primal interior point method for a class of control constrained optimal control problems, Computational Optimization and Applications, 41 (2008), pp. 127–145, https://doi.org/10.1007/s10589-007-9088-y.
  * [71] J. Xia, P. E. Farrell, and F. Wechsung, Augmented Lagrangian preconditioners for the Oseen–Frank model of nematic and cholesteric liquid crystals, BIT Numerical Mathematics, 61 (2021), pp. 607–644, https://doi.org/10.1007/s10543-020-00838-9.
  * [72] J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Review, 34 (1992), pp. 581–613, https://doi.org/10.1137/1034116.
  * [73] J. Xu, The method of subspace corrections, Journal of Computational and Applied Mathematics, 128 (2001), pp. 335–362, https://doi.org/10.1016/S0377-0427(00)00518-5.
  * [74] K. Yaji, T. Yamada, M. Yoshino, T. Matsumoto, K. Izui, and S. Nishiwaki, Topology optimization using the lattice Boltzmann method incorporating level set boundary expressions, Journal of Computational Physics, 274 (2014), pp. 158–181, https://doi.org/10.1016/j.jcp.2014.06.004.
  * [75] F. Zhang, The Schur complement and its applications, vol. 4, Springer Science & Business Media, 2006, https://doi.org/10.1007/b105056.
  * [76] S. Zhou and Q. Li, A variational level set method for the topology optimization of steady-state Navier–Stokes flow, Journal of Computational Physics, 227 (2008), pp. 10178–10195, https://doi.org/10.1016/j.jcp.2008.08.022.