# Skorokhod energy of planar domains

Mrabet Becher Institut Préparatoire aux Etudes des Ingénieurs de Monastir,
Tunisia. Maher Boudabra King Fahd University of Petroleum and Minerals, KSA.
Fethi Haggui Institut Préparatoire aux Etudes des Ingénieurs de Monastir,
Tunisia.

###### Abstract

In this work, we introduce the Skorokhod energy of a simply connected domain.
We show that among all domains solving the planar Skorokhod embedding problem,
Gross solution generates the domain with the minimal Skorokhod energy.

##  1 Introduction and main results

In a pretty recent paper appeared in $2019$ [8], the author R. Gross
considered an interesting planar version of the Skorokhod problem, which was
originally formulated in dimension one in $1961$(we refer the reader to [14]
for a concise survey of the linear version). The statement of the planar
Skorokhod embedding problem is as follows : Given a distribution $\mu$ with
zero mean and finite second moment, is there a simply connected domain $U$
(containing the origin) such that if $Z_{t}=X_{t}+Y_{t}i$ is a standard planar
Brownian motion, then $X_{\tau}=\Re(Z_{\tau})$ has the distribution $\mu$,
where $\tau$ is the exit time from $U$. Gross provided an affirmative answer
with a smart and explicit construction of the domain. In addition, the
generated domain has an exit time of finite average. Gross has also observed
that the solution domain is not unique without further conditions. The core
tool of Gross idea is the following theorem.

######  Theorem 1 (Lévy’s theorem).

Let $f:U\rightarrow\mathbb{C}$ be a non constant analytic function and $Z_{t}$
be a planar Brownian motion running inside $U$. Then there is a planar
Brownian motion $W_{t}$ such that $f(Z_{t})=W_{\sigma(t)}$ with

$\sigma(t)=\int_{0}^{t}|f^{\prime}(Z_{s})|^{2}ds$

and $t\in[0,\tau_{U}]$.

Theorem 1 is commonly known as the conformal invariance principle. That is,
the image of a planar Brownian motion under the action of non constant
analytic function is yet another planar Brownian motion, except it runs in a
different speed. With the same notation, it is equivalent to say that
$f(Z_{\sigma^{-1}(t)})$ is planar Brownian motion. The random time
$\sigma(\tau_{U}^{-})$ is referred to as the projection of the exit time
$\tau_{U}$. It is not always the case that $\sigma(\tau_{U}^{-})$ is the exit
time of $f(U)$ unless $f$ is proper which is the case when $f$ is univalent.
In particular, the law of $f(Z_{\tau_{U}})$ is obtainable from that of
$Z_{\tau_{U}}$. The idea of Gross is based on this fact. In order to solve the
problem, it is enough to construct a univalent map acting on the unit disc
such that $\Re(f(e^{ti}))$ samples as $\mu$. As $Z_{\tau_{\mathbb{D}}}$ is
uniformly distributed on the unit circle, then the distribution of
$\Re(f(e^{ti}))$ is the same as the distribution of the real part of a stopped
planar Brownian motion on the boundary of $f(\mathbb{D})$.

In [4], the two authors, Boudabra and Markowsky gave geometric conditions
which, when imposed, the domain becomes unique. Moreover, the frame of treated
distributions has been drastically enlarged to cover all of finite
$p^{th}$moments where $p$ belongs to $(1,+\infty)$. Later on, the two authors
published a second paper which generate a new category of domains that solve
the Skorokhod embedding problem. As in their first paper, the new solution
applies to any distribution of finite $p^{th}$moment ($p>1$). Furthermore, a
uniqueness criterion is given. We shall use their terminology “$\mu$-domain”
to tag any simply connected domain that solves the problem. The differences
between the two solutions are summarized as follows :

  * •

Gross construction :

    * –

$U$ is symmetric over the real line.

    * –

$U$ is $\Delta$-convex, i.e the segment joining any point of $U$ and its
symmetric one over the real axis remains inside $U$.

    * –

If $\mu(\\{x\\})>0$ then $\partial U$ contains a vertical line segment
(possibly infinite).

  * •

Boudabra-Markowsky construction :

    * –

$U$ is $\Delta^{\infty}$-convex, i.e the upward half ray at any point of $U$
remains inside $U$.

    * –

If $\mu(\\{x\\})>0$ then $\partial U$ contains a downward half ray at $x$.

If $(a,b)$ is not supported by $\mu$, i.e $(a,b)$ is null set for $\mu$ 111 It
means that the c.d.f of $\mu$ is constant on $(a,b)$., then the two solutions
would contains the vertical strip $(a,b)\times(-\infty,+\infty)$. As per the
uniqueness, the authors essentially showed that $\Delta^{\infty}$-convexity or
symmetry and $\Delta$-convexity guarantee uniqueness. We provide the following
two examples illustrating the differences of the two methods.

Figure 1.1: For the uniform distribution on $(-1,1)$, the left domain is
Boudabra-Markowsky solution while Gross solution is on the right.

Figure 1.2: The domains generated from the shifted-scaled arcsine distribution
$d\mu(x)=\frac{1}{\pi\sqrt{1-x^{2}}}dx$. Gross method generates the unit disc
while Boudabra-Markowsky method generates the right domain.

Surprisingly, the second solution turns out to solve completely a separate
optimization problem posed and partially answered in [12]. The question
therein was to seek the $\mu$-domain having the lowest or the largest
principal Dirichlet eigenvalue. The authors in [12] gave the $\mu$-domain of
minimal principal Dirichlet eigenvalue for $\mu$ being the uniform
distribution on $(-1,1)$. Their argument leans on the geometry of the domain
they have constructed. in [3], Boudabra and Markowsky noticed that this
argument applies perfectly to their solution. In particular, they recovered
the $\mu$-domain provided in [12]. The question about the largest principal
Dirichlet eigenvalue is still open. Throughout this work, the following
notation will be adopted :

  * •

$\|f\|_{p}=\left(\frac{1}{2\pi}\int_{0}^{2\pi}|f(t)|^{p}dt\right)^{\frac{1}{p}}$
where $p>0$.

  * •

$L_{2\pi}^{p}$ : the space of $2\pi$-periodic functions $f$ such that
$\|f\|_{p}$ is finite.

  * •

$W_{2\pi}^{1,p}$ the Sobolev space for the $2\pi$-periodic functions.

  * •

The word “domain” means a simply connected domain.

  * •

$(Z_{t})_{t}$ refers to a standard planar Brownian motion.

  * •

$\tau_{U}$ is the exit time of $(Z_{t})_{t}$ from a domain $U$.

In the literature and in the context of physics, the kinetic energy of a
function $f$, refers to the energy associated with the motion of particles
described by that function. It is a concept commonly encountered in classical
mechanics, where the kinetic energy of an object is determined by its mass and
velocity. In mathematical terms, the kinetic energy of a (real valued)
$2\pi$-periodic function $f$, which we denote by $\mathfrak{E}(f)$, is
proportional to the average of the squared norm of its velocity. That is

$\mathfrak{E}(f):={\textstyle\frac{1}{2}}{\textstyle\frac{1}{2\pi}}\int_{0}^{2\pi}|f^{\prime}(t)|^{2}dt={\textstyle\frac{1}{2}}\|f^{\prime}\|_{2}^{2}$

where $f^{\prime}$ is the distributional derivative of $f$. In particular, the
energy $\mathfrak{E}(f)$ is infinite for any function $f$ with jumps. To see
this, remark that if $f$ has a discontinuity at $t_{0}$ then
$f^{\prime}(t_{0})=\delta_{t_{0}}\notin L_{2\pi}^{2}$. From Fourier series
theory, $f$ has a Fourier expansion

$f(\theta)=\sum_{n=1}^{+\infty}a_{n}\cos(n\theta)+b_{n}\sin(n\theta)$

and hence

$f^{\prime}(\theta)=\sum_{n=1}^{+\infty}-na_{n}\sin(n\theta)+nb_{n}\cos(n\theta)$

in the distribution sense. Parseval’s identity applied to $f^{\prime}$ yields

$\mathfrak{E}(f)={\textstyle\frac{1}{4}}\sum_{n=1}^{+\infty}n^{2}(a_{n}^{2}+b_{n}^{2}).$

When the kinetic energy of $f$ is finite then the Fourier coefficients decay
faster than $\frac{1}{n}$ at infinity. In this context, we introduce the
Skorokhod energy of a domain $U$ containing the origin. Before that, we state
the following lemma.

######  Lemma 2.

Let $f$ be a univalent function mapping $\mathbb{D}$ onto a domain $U$ and
fixing the origin. Then the quantity $\mathfrak{E}(\Re(f(e^{ti})))$ is
independent of $f$.

###### Proof.

The existence of a map $f$ is guaranteed by the famous Riemann mapping
theorem. Moreover, any such two functions $f$ and $g$ are related by some
rotation, i.e $f(z)=g(e^{\theta i}z)$. The result follows by a standard change
of variable. ∎

Note that if $f(z)=\sum_{n=1}^{+\infty}c_{n}z^{n}$ then

$\mathfrak{E}(\Re(f(e^{ti})))=\frac{1}{4}\sum_{n=1}^{+\infty}n^{2}|c_{n}|^{2}.$

Note that $f(e^{ti})$ is interpreted as the radial limit function of $f$. In
the second section, we will see that such a limit does always exist for
univalent functions. Now we give the definition of the Skorokhod energy of a
domain $U$.

######  Definition 3.

The Skorokhod energy of $U$, denoted by $\Lambda(U)$, is the quantity
$\mathfrak{E}(\Re(f(e^{ti})))$ where $f$ is any univalent function mapping
$\mathbb{D}$ onto $U$ and fixing the origin.

By virtue of lemma 2, the Skorokhod energy of $U$ is well defined since it
does not depend on the choice of the mapping function. For example, the
Skorokhod energy of a centered disc with radius $r$ is $\frac{r^{2}}{4}$. The
next result provides a criterion to test the finiteness of the Skorokhod
energy.

######  Proposition 4.

A finite energy implies finite perimeter and area.

A consequence of proposition 4 is that unbounded domains or bounded ones with
irregular boundaries are of infinite Skorokhod energy. In particular, the
$\mu$-domains generated by the technique in [3] have infinite Skorokhod energy
as they are unbounded. Besides, all $\mu$-domains generated from unbounded
distributions are of infinite Skorokhod energy. Note the property of having a
finite Skorokhod energy is not a size-hereditary property, i.e if $\Lambda(U)$
is finite and $V\subset U$ then $\Lambda(V)$ is not necessarily finite.

Now we are almost ahead of our main result. Let $U$ be a domain and run
$(Z_{t})_{t}$ be a standard planar Brownian motion running inside $U$ and let
$\mu_{U}$ be the distribution of $\Re(Z_{\tau_{U}})$, the real part of $Z_{t}$
when it hits the boundary of $U$. Let $G$ be the Gross domain, the
$\mu_{U}$-domain obtained by applying Gross method to $\mu_{U}$. In
particular, both domains $U$ and $G$ are $\mu_{U}$-domains.

######  Theorem 5.

Let $U$ be a domain and let $G$ be the Gross domain generated by $\mu_{U}$.
Then $G$ has less Skorokhod energy than $U$, i.e

$\Lambda(G)\leq\Lambda(U).$

Theorem 5 enables us to define Skorokhod energy of a distribution $\mu$ to be
$\Lambda(G)$ where $G$ is the Gross domain generated by $\mu$. Now, we provide
some explained examples:

  1. 1.

By the uniqueness theorem in [4], the unit disc is the Gross domain of the
scaled-shifted arc-sine law given by
$d\mu(x)=\frac{dx}{\pi\sqrt{1-x^{2}}}\mathbf{1}_{\\{x\in(-1,1)\\}}$. That is,
the Skorokhod energy of $\mu$ is $\Lambda(\mathbb{D})=\frac{1}{4}$.

  2. 2.

Let $\mu$ be the uniform distribution on $(-1,1)$. From [8], the corresponding
Gross $\mu$-domain is the image of the unit disc by the map

$\alpha(z)=-\frac{8}{\pi^{2}}\sum_{n=1}^{+\infty}\frac{z^{2n-1}}{(2n-1)^{2}}.$

Therefore the Skorokhod energy of $\mu$ is

$\frac{8}{\pi^{2}}\sum_{n=1}^{+\infty}\frac{1}{(2n-1)^{2}}=\frac{8}{\pi^{2}}\frac{\pi^{2}}{8}=1.$

  3. 3.

Let $\mathscr{P}_{m}$ be the regular polygon of vertices $e^{\nicefrac{{2\pi
ki}}{{m}}}$, $k=0,...,m-1$ where $m$ is an even integer greater than $2$. A
conformal map from the unit disc onto $\mathscr{P}_{m}$ and fixing the
$e^{\nicefrac{{2\pi ki}}{{m}}}$’s, is given by

$f_{m}(z)=\frac{m}{\beta(\nicefrac{{1}}{{m}},\nicefrac{{(m-2)}}{{m}})}\sum_{n=0}^{\infty}\frac{(\nicefrac{{2}}{{m}})_{n}z^{mn+1}}{n!(mn+1)}$

where $(x)_{n}:=x(x+1)\cdots(x+n-1)$ and $\beta$ refers to the beta function
(See [13]). As $\mathscr{P}_{m}$ fulfills the uniqueness criterion then it is
the Gross domain generated from the the distribution of
$\Re(Z_{\tau_{\mathscr{P}_{m}}})$. The energy of such distribution is infinite
as

$\sum_{n=0}^{\infty}n^{2}\left(\frac{(\nicefrac{{2}}{{m}})_{n}}{n!(mn+1)}\right)^{2}=+\infty.$

Figure 1.3: The graph of the c.d.f of the distribution of $\Re(Z_{\tau_{S}})$
where $S=\frac{2e^{\nicefrac{{-\pi i}}{{4}}}}{\sqrt{2}}f_{4}(\mathbb{D})$ is
the square of vertices $\\{\pm 1,\pm i\\}$. The presence of the two jumps at
$\pm 1$ is reflected by the atoms of the distribution.

In a very recent manuscript [2], the two authors investigated the geometric
map transforming $U$ to $G$. They proposed to call it Brownian symmetrization
(inspired by Steiner symmetrization) denoted by $\mathfrak{B}$. They proved
these two transformations are different and included some extremal problems.
Taking into consideration their geometric perspective, theorem 5 reads as :
Brownian symmetrization of a bounded domain does not increase the Skorokhod
energy. In other words

$\Lambda(\mathfrak{B}(U))\leq\Lambda(U).$

##  2 Tools and proofs

In this section, we provide all required ingredients of the results stated in
the previous section.

######  Definition 6.

The Hilbert transform of a $2\pi$\- periodic function $f$ is defined by

$H\\{f\\}(x):=PV\left\\{\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x-t)\cot(\frac{t}{2})dt\right\\}=\lim_{\eta\rightarrow
0}\frac{1}{2\pi}\int_{\eta\leq|t|\leq\pi}f(x-t)\cot(\frac{t}{2})dt$

where $PV$ denotes the Cauchy principal value, which is required here as the
trigonometric function $t\longmapsto\cot(\cdot)$ has poles at $k\pi$ with
$k\in\mathbb{Z}$. The Hilbert transform is a bounded operator on
$L_{2\pi}^{p}$ with $p>1$. More precisely we have

######  Theorem 7.

[6] If $f$ is in $L_{2\pi}^{p}$ then $H\\{f\\}$ does exit almost every where
for $p>1$. Furthermore

$\|H\\{f\\}\|_{p}\leq\lambda_{p}\|f\|_{p}$ (2.1)

for some positive constant $\lambda_{p}$.

The case $p=1$ is irregular and unfortunately $H$ becomes unbounded. For all
proofs and properties of the Hilbert transform, we refer the reader to [6,
11]. When $p=2$, which is related to our frame, Hilbert transform becomes an
isometry. That is

$\|H\\{f\\}\|_{2}=\|f\|_{2}.$ (2.2)

The equality (2.2) can be easily understood if one knows that
$H\\{\cos\\}=\sin$ and $H\\{\sin\\}=-\cos$.

######  Definition 8.

Let $f$ be an analytic function acting on the unit disc and $p>0$. The
$p^{th}$\- Hardy norm of $f$ is defined by

$H_{p}(f):=\sup_{0\leq r<1}\left\\{\frac{1}{2\pi}\int_{0}^{2\pi}|f(re^{\theta
i})|^{p}d\theta\right\\}^{\frac{1}{p}}=\sup_{0\leq r<1}N_{r}\\{f\\}.$ (2.3)

The Hardy norm is well defined since the quantity $N_{r}(f)$ is non decreasing
in terms of $r$ [18]. A crucial result about Hardy norms is that, if
$H_{p}(f)$ is finite then $f$ has a radial extension to the boundary
$f^{\dagger}(\xi):=\lim_{r\rightarrow 1}f(r\xi)$ exists a.e for all
$\xi\in\partial\mathbb{D}$. Moreover

$\|f^{\dagger}\|_{p}=H_{p}(f).$

For the sake of brevity, we keep writing $f(e^{ti})$ instead of
$f^{\dagger}(e^{ti})$. In [5], the author gave a very fancy result revealing
the tight connections between the theory of planar Brownian motion and the
theory of analytic functions. That is, he showed the following equivalences

$H_{p}(f)<\infty\Longleftrightarrow\mathbf{E}(\tau^{\nicefrac{{p}}{{2}}})<\infty\Longleftrightarrow\mathbf{E}(\sup_{0\leq
t\leq\tau}|Z_{t}|^{p})\Rightarrow\mathbf{E}(|Z_{t}|^{p})$

where $f$ is univalent and $\tau$ is the exit time of a planar Brownian motion
from $f(\mathbb{D})$. The last implication turns to be an equivalence provided
that $\mathbf{E}(\ln(1+\tau))<\infty$. This requirement has been shown to be
redundant when $f$ is univalent [1]. The following theorem concerns the range
of $p$ for which the $p^{th}$ Hardy norm of a univalent function is finite. In
particular, it gives sense to the statement in lemma 2.

######  Theorem 9.

[7] If $f$ is univalent then $H_{p}(f)$ is finite for all $p<\frac{1}{2}$.

Theorem 9 has triggered the question about the supremum of $p$ such that
$H_{p}(f)$ is finite. Such a number is called the Hardy number [9]. As the
Hilbert transform maps $\cos$ to sin and sin to cos, then it is not difficult
to see that

$f(e^{ti})=\Re\\{f((e^{ti})\\}+H\\{\Re\\{f(e^{ti})\\}\\}i.$

In other words, the imaginary part of $f(e^{ti})$ is simply the Hilbert
transform of its real part. By (2.2), we get

$H_{2}(f)^{2}=\|f\|_{2}^{2}=2\|\Re\\{f\\}\|_{2}^{2}.$

In particular

$\mathfrak{E}(f)=\frac{1}{4}H_{2}(f^{\prime})^{2}.$

###  2.1 Proof of proposition 4

The planar version of the isoperimetric inequality states that

$\text{area}(U)\leq{\textstyle\frac{1}{4\pi}}\ell(U)^{2}$

where $\text{area}(U)$ and $\ell(U)$ are the area of $U$ and its perimeter
[15, 17]. That is, it is enough to show that finite Skorokhod energy implies
finite perimeter. Let $f(z)=\sum_{n=1}^{+\infty}c_{n}z^{n}$ be a univalent
function mapping $\mathbb{D}$ onto $U$. The perimeter of $U$ is

$\ell(U)=2\pi H_{1}(f^{\prime}).$

which covers non smooth boundaries as well. the Cauchy–Bunyakovsky–Schwarz
inequality yields

$\displaystyle\ell(U)^{2}$ $\displaystyle\leq
2\pi\int_{-\pi}^{\pi}|f^{\prime}(e^{ti})|^{2}dt$ $\displaystyle\leq
4\pi^{2}H_{2}(f^{\prime})^{2}$
$\displaystyle\overset{\eqref{eq:isometry}}{=}8\pi^{2}\|\Re\\{f^{\prime}(e^{ti})\\}\|_{2}^{2}$
$\displaystyle=16\pi^{2}\Lambda(U).$

The result follows.

Note that an inequality between the area and the Skorokhod energy can be
obtained without involving the perimeter. In fact, a standard identity for the
area of $U$ is given by

$\text{area}(U)=\pi\sum_{n=1}^{+\infty}n|c_{n}|^{2}.$ (2.4)

In particular we get

$\text{area}(U)\leq
4\pi{\textstyle\frac{1}{4}}\sum_{n=1}^{+\infty}n^{2}|c_{n}|^{2}=4\pi\Lambda(U).$

The three quantities $\text{area}(U),\ell(U),\Lambda(U)$ are related by the
inequality

$\text{area}(U)\leq{\textstyle\frac{1}{4\pi}}\ell(U)^{2}\leq 4\pi\Lambda(U).$

The constants ${\textstyle\frac{1}{4\pi}}$ and $4\pi$ are optimal as for
$U=\mathbb{D}$ as all the terms are equal to $\pi$.

###  2.2 Proof of theorem 5

The proof of theorem 5 is based on the symmetric decreasing rearrangement. For
a survey about the subject as well as the main properties, we refer the reader
to [10]. Let $f$ be a non negative function and let
$\xi_{t}^{f}=\frac{|f^{-1}(t,\infty)|}{2}$. When $f$ is clear from the
context, we shall omit the superscript. It is obvious that the map
$t\mapsto\xi_{t}$ is non increasing.

######  Definition 10.

The symmetric decreasing rearrangement of $f$, denoted by $f^{*}$, is defined
by

$f^{*}(x)=\int_{0}^{+\infty}1_{(-\xi_{t},\xi_{t})}(x)dt$ (2.5)

The set $(-\xi_{t},\xi_{t})$ is the symmetrized of $f^{-1}(t,\infty)$. Formula
(2.5) shows that the function $f^{*}$ is even, i.e $f^{*}(x)=f^{*}(|x|)$. As
$s\mapsto 1_{(-\xi_{s},\xi_{s})}(|x|)$ is non increasing then

$f^{*}(|x|)=\sup\\{s\mid|x|<\xi_{s}={\textstyle\frac{|f^{-1}(s,+\infty)|}{2}}\\}.$
(2.6)

The original definition of $f^{*}$ (2.5) applies to non negative functions,
however the identity (2.6) extends to functions with negative values. That is,
if $f$ is a function bounded below then we define the symmetric decreasing
rearrangement of $f$ by (2.6). In particular we find

$f^{*}=g^{*}+\inf(0,\inf f)$

with $g=f-\inf(0,\inf f)$. Hence, when $f$ hits the negative real line then it
symmetric decreasing rearrangement is obtained by a simple shift. A crucial
property of the symmetric decreasing rearrangement is the equimeasurability.
In other words, the level sets $\\{x\mid f(x)>t\\}$ and $\\{x\mid
f^{*}(x)>t\\}$ have the same Lebesgue measure. From a probabilistic
perspective, this property interprets as : A random variable $\vartheta$ and
its symmetric decreasing rearrangement $\vartheta^{*}$ share the same
probability law. Note that $\vartheta$ and $\vartheta^{*}$ are not coupled,
i.e they do not necessarily act on the same sample space, unless it is
symmetric of course. The next lemma characterizes a function in terms of its
variation and its probability distribution, i.e the measure of its level sets.

######  Lemma 11.

Let $\vartheta,\varrho:I\rightarrow\mathbb{R}$ be two non decreasing functions
with the same distribution. Then $\vartheta=\varrho$ a.e.

###### Proof.

Since $\vartheta$ and $\varrho$ have the same distribution then by
monotonicity we get

$\displaystyle|\vartheta^{-1}(s,\infty)|$
$\displaystyle=|\varrho^{-1}(s,\infty)|$
$\displaystyle=1-\inf\\{x\mid\vartheta(x)>s\\}$
$\displaystyle=1-\inf\\{x\mid\varrho(x)>s\\}.$

Thus $\inf\\{x\mid\vartheta(x)>s\\}=\inf\\{x\mid\varrho(x)>s\\}$. By lemma 11
in [3], we get $\vartheta=\varrho$ a.e. ∎

######  Remark 12.

Lemma 11 can be generalized in the following way : If $I=\biguplus I_{j}$
where $\vartheta,\varrho$ share the same distribution and have the same
monotonicity on each $I_{j}$, then $\vartheta=\varrho$ a.e.

Let $\mu$ be a bounded distribution measure of c.d.f $F$, which assumed to
have a zero average. The underlying quantile function of $\mu$ is defined for
every $u\in(0,1)$ by

$Q(u)=F^{-1}(u)=\inf\\{x\mid F(x)\geq u\\}=\sup\\{x\mid F(x)<u\\}\,\,\,(a.e).$
(2.7)

That is, $Q$ is the pseudo-inverse of $F$. In particular, $Q$ simplifies to
the standard inverse function when $F$ is continuous. The interchange between
$\inf$ and $\sup$ in (2.7) is justified by the fact that the set of
discontinuity of $F$ is at most countable. The quantile function is very handy
to generate random variables. That is, if $U$ is uniformly distributed on
$(0,1)$ then $Q(U)$ samples as $\mu$. The following proposition yields a clue
relation between $Q$ and its symmetric decreasing rearrangement $Q^{*}$.

######  Proposition 13.

We have

$Q^{*}(x)=Q(1-2|x|)\,\,\,\,(a.e)$ (2.8)

where $x\in(-\frac{1}{2},\frac{1}{2})$.

###### Proof.

We have

$\displaystyle Q^{*}(|x|)$
$\displaystyle=\sup\\{s\mid|x|<{\textstyle\frac{|Q^{-1}(s,\infty)|}{2}}\\}$
$\displaystyle=\sup\\{s\mid|x|<{\textstyle\frac{1-F(s)}{2}}\\}$
$\displaystyle=\sup\\{s\mid F(s)<1-2|x|\\}$ $\displaystyle=F^{-1}(1-2|x|)$
$\displaystyle=Q(1-2|x|).$

∎

The identity (2.8) is consistent with the definition as $Q^{*}$ is defined on
$(-\frac{1}{2},\frac{1}{2})$ which is the symmetrized of $(0,1)$. The
following corollary follows immediately from proposition 13.

######  Corollary 14.

If $\psi:\theta\in(-\pi,\pi)\mapsto\mathbb{R}$ is a random variable sampling
as $\mu$ then for almost every $\theta\in(-\pi,\pi)$

$\psi^{*}(\theta)=Q(1-{\textstyle\frac{|\theta|}{\pi}}).$

In [8], the author showed the following result, cornerstone of his solution
for the planar Skorokhod embedding problem.

######  Theorem 15.

Let $(a_{n})_{n}$ be the sequence of Fourier coefficients of
$\theta\in(-\pi,\pi)\mapsto Q(\frac{|\theta|}{\pi})$ then

$\Psi(z)=\sum_{n=1}^{\infty}a_{n}z^{n}$

is univalent on the unit disc.

The Fourier expansion of $Q(\frac{|\theta|}{\pi})$ is simply
$\sum_{n=1}^{\infty}a_{n}\cos(n\theta)$ (no $\sin$ terms), and Gross smart
idea was to generate the power series out of the Fourier expansion. Now,
remark that the Fourier coefficients of $Q(\frac{|\theta|}{\pi})$ and
$Q^{*}(\theta)=Q(1-\frac{|\theta|}{\pi})$ differ by a $(-1)^{n}$ factor. That
is, if we denote by $\mathscr{S}(z)$ the power series generated by the Fourier
coefficients of $Q^{*}$ then we get

$\mathscr{S}(z)=\Psi(-z)$

and hence $\mathscr{S}$ is univalent, and it obviously maps the unit disc to
$\Psi(\mathbb{D})$, the same $\mu$-domains obtained by Gross technique. The
last ingredient of the proof is the famous Pólya-Szegö inequality.

######  Theorem 16.

[10] If $f$ is in $W^{1,p}$ then $f^{*}\in W^{1,p}$ and we have

$\|\nabla f^{*}\|_{p}\leq\|\nabla f\|_{p}.$

Pólya-Szegö inequality states that symmetric decreasing rearrangement does not
boost the norm of the gradient. In particular, the kinetic energy of a
function does not increase under symmetrization. The norm of the gradient (up
to to some factor) appears frequently in areas related to potential theory
like Electrostatics [16].

Our main result follows now. Let $\Phi$ be any univalent function mapping
$\mathbb{D}$ onto $U$ and fixing the origin. We have

$\Lambda(U)=\mathfrak{E}(\Re\\{\Phi(e^{\theta i})\\})$

and

$\Lambda(G)=\mathfrak{E}(\Re\\{\Psi(e^{\theta i})\\}).$

The conformal invariance principle implies that $\Phi(e^{\theta
i}),\Re(Z_{\tau_{U}}),\Psi(e^{\theta i}),\Re(Z_{\tau_{G}})$ sample as
$\mu_{U}.$ By corollary 14, the symmetric decreasing rearrangement of
$\Re\\{\Phi(e^{\theta i})\\}$ is $Q(1-\frac{|\theta|}{\pi})$ where $Q$ is the
quantile function of $\mu_{U}$. In particular, the univalent function
$\mathscr{S}$ maps the unit disc to Gross domain, solution of the Skorokhod
problem. By Pólya-Szegö inequality inequality we get

${\textstyle\frac{1}{2}}\|{\Re\\{\mathscr{S}(e^{ti})\\}^{*}}^{\prime}\|_{2}^{2}={\textstyle\frac{1}{2}}\|Q^{\prime}\|_{2}^{2}=\Lambda(G)\leq{\textstyle\frac{1}{2}}\|\Re\\{\Phi(e^{ti})\\}\|_{2}=\Lambda(U)$

which proves our claim.

##  3 Comments and further directions

As per Gross technique, the Skorokhod energy of any discrete distribution
$\mu=\sum_{j=1}^{n}p_{j}\delta_{x_{j}}$is infinite as its corresponding Gross
domain is unbounded. More precisely, it is a vertical strip $\\{\min
x_{j}<\Re(z)<\max x_{j}\\}$ minus some vertical half rays at the remaining
$x_{j}$’s. On the other hand, if $\mu$ is continuous of p.d.f $\rho$, then the
Skorokhod energy of $\mu$ simplifies to

$\Lambda(\mu)=\frac{1}{2}\int_{0}^{1}Q^{\prime}(x)^{2}dx=\frac{1}{2}\int_{0}^{1}\frac{1}{\rho(Q(x))^{2}}dx.$

where $\rho$ is the p.d.f of $\mu$.

Theorem 5 can be seen as follows (with the same notations): If
$\Phi=\sum_{n=1}^{+\infty}b_{n}z^{n}$ and
$\Psi=\sum_{n=1}^{+\infty}a_{n}z^{n}$ then

$\sum_{n=1}^{+\infty}n^{2}|a_{n}|^{2}\leq\sum_{n=1}^{+\infty}n^{2}|b_{n}|^{2}.$

Theorem 5 persuades us to believe that the Brownian symmetrization does not
increase neither the area nor the perimeter [2]. Using the area formula 2.4,
it implies in particular that

$\sum_{n=1}^{+\infty}n|a_{n}|^{2}\leq\sum_{n=1}^{+\infty}n|b_{n}|^{2}.$

However, the effect of Gross construction on the principal Dirichlet
eigenvalue is still far from being decided. In this scope, we think the
problem of finding optimal $\mu$-domains with respect to a certain property,
like the area, the perimeter etc is worth studying. In other words, find the
solution of

$\inf_{U:\mu\text{-domain}}\mathcal{J}(U)$ (3.1)

where $\mathcal{J}$ is some functional ($\mu$ being a bounded distribution).
This would definitely open new avenues of research that foster more the
connections between geometry, PDE’s and planar Brownian motion. So far and as
mentioned in the introduction, one problem is treated in this direction. It is
when $\mathcal{J}$ outputs the principal Dirichlet eigenvalue, where the
solution of 3.1 is given in [3].

## References

  * [1] M. Boudabra.  On some problems related to exit times of planar Brownian motion.  PhD thesis, Monash University, 2021.
  * [2] M. Boudabra and K. Hamza.  Brownian symmetrization of planar domains, 2024.
  * [3] M. Boudabra and G. Markowsky.  A new solution to the conformal Skorokhod embedding problem and applications to the Dirichlet eigenvalue problem.  Journal of Mathematical Analysis and Applications, 491(2):124351, 2020.
  * [4] M. Boudabra and G. Markowsky.  Remarks on Gross’ technique for obtaining a conformal Skorohod embedding of planar Brownian motion.  Electronic Communications in Probability, 2020.
  * [5] D. Burkholder.  Exit times of Brownian motion, harmonic majorization, and Hardy spaces.  Advances in Mathematics, 26(2):182–205, 1977.
  * [6] P. Butzer and R. Nessel.  Hilbert transforms of periodic functions.  In Fourier Analysis and Approximation, pages 334–354. Springer, 1971.
  * [7] P. L Duren.  Univalent functions, volume 259.  Springer Science & Business Media, 2001.
  * [8] R. Gross.  A conformal Skorokhod embedding.  Electronic Communications in Probability, 2019.
  * [9] L.J Hansen.  Hardy classes and ranges of functions.  The Michigan Mathematical Journal, 17(3):235–248, 1970.
  * [10] S. Kesavan.  Symmetrization and applications, volume 3.  World Scientific, 2006.
  * [11] F. King.  Hilbert transforms.  Cambridge University Press Cambridge, 2009.
  * [12] P. Mariano and H. Panzo.  Conformal Skorokhod embeddings and related extremal problems.  Electronic Communication in Probability, 25, 2020.
  * [13] G. Markowsky.  On the distribution of planar Brownian motion at stopping times.  Annales Academiæ Scientiarum Fennicæ Mathematica, 2018.
  * [14] J. Obłój.  The Skorokhod embedding problem and its offspring.  Probability Surveys, pages 321 – 392, 2004.
  * [15] R. Osserman.  The isoperimetric inequality.  Bulletin of the American Mathematical Society, 84(6):1182–1238, 1978.
  * [16] L. E. Payne.  Isoperimetric inequalities and their applications.  SIAM review, 9(3):453–488, 1967.
  * [17] G. A Philippin and G. Porru.  Isoperimetric inequalities and overdetermined problems for the Saint-Venant equation.  New Zealand J. Math, 25:217–227, 1996.
  * [18] W. Rudin.  Real and complex analysis (3r ed.).  McGraw-Hill Education, 2001.


# Code generation for productive portable scalable finite element simulation
in Firedrake.

Jack D. Betteridge, Patrick E. Farrell and David A. Ham J. D. Betteridge and
D. A. Ham are at Imperial College London. P. E. Farrell is at the University
of Oxford.

###### Abstract

Creating scalable, high performance PDE-based simulations requires a suitable
combination of discretizations, differential operators, preconditioners and
solvers. The required combination changes with the application and with the
available hardware, yet software development time is a severely limited
resource for most scientists and engineers. Here we demonstrate that
generating simulation code from a high-level Python interface provides an
effective mechanism for creating high performance simulations from very few
lines of user code. We demonstrate that moving from one supercomputer to
another can require significant algorithmic changes to achieve scalable
performance, but that the code generation approach enables these algorithmic
changes to be achieved with minimal development effort.

## Introduction

The simulation of continuum mechanics problems through the numerical solution
of partial differential equations (PDEs) is a classical high performance
computing task. The three-dimensional multiscale nature of science and
engineering challenges in fields such as fluid mechanics, solid mechanics, and
electromagnetics means that the number of variables and operations required
can become very large very quickly. The ubiquitous nature of these systems
means that many scientists and engineers need to solve problems in this class.
However, the specific equations to be solved, the discretizations and solver
strategies employed, and the coupling of these to other processes, all vary
greatly from case to case. It is also seldom the case that solving a PDE is
the whole problem: scientists and engineers need to quantify uncertainty in
solutions, compute their sensitivity to inputs, assimilate observed data, or
optimize parameters.

The combination of all of these factors means that bespoke simulations are the
norm: most simulations are put together for a particular purpose by one or a
small group of scientists or engineers, and run a limited number of times
before the problem specification changes in a way which demands alterations to
the equations being solved or the algorithms being used to solve them. To add
to this complexity, these users will rarely be able to design or procure their
own supercomputer. Rather, they will need to make effective use of whichever
institutional or national facilities they can access, making the required
algorithmic changes as they move from one machine to another. Even the
emergence of cloud computing resources offers little comfort, since the ever-
changing hardware landscape means that the performance characteristics of
processors such as memory bandwidth, core count, and vector length evolve
rapidly and may vary between providers, or even between the classes of cloud
node suitable for simulations of different sizes.

Making high performance parallel simulation viable for this long tail of
simulation science and engineering is a different proposition from the large-
scale operational simulations carried out in fields such as weather
forecasting and seismic inversion, in which the same computational problem is
solved an enormous number of times for different data inputs. In those
operational fields, the cost of laborious implementation and optimization work
can be amortized over countless simulation runs so that the cost per
simulation is not dominated by development. Conversely, for the typical
simulation scientist, the requirement is for high performance for as close to
zero programmer effort as can be achieved. Further, the cost of changing any
aspect of the simulation specification or algorithm must be minimal, lest the
code become fossilized and unable to adapt to the next change in task or
available hardware.

Here we demonstrate the utility of the Firedrake automated finite element
system in addressing this need. Firedrake enables users to write high-level
Python code in a symbolic language which reflects the mathematics of the
problem. It integrates tightly with the Portable Extensible Toolkit for
Scientific Computation (PETSc) [1] to enable runtime programmable
preconditioners and solvers that reach back to the discretization as required.
We will demonstrate the capabilities of Firedrake in providing productive
performance portability by taking a fluids simulation tuned for ARCHER, the
previous UK national supercomputer, and porting it to the ARM-based Isambard
and ARCHER2, the AMD-based new UK national supercomputer. We find that
algorithmic choices made for one platform are suboptimal on other platforms,
but that working from a high-level mathematical specification of the problem
and generating the implementation automatically enables the necessary
algorithmic changes to be made at minimal user cost.

## Firedrake

The Firedrake automated finite element system [2]
(https://firedrakeproject.org) is a Python package which generates numerical
solutions to PDEs from a very high level mathematical specification provided
by the user. The PDEs to be solved by Firedrake are specified in the Unified
Form Language (UFL) [3], a specialized computer algebra language also employed
by FEniCS [4] and DUNE [5]. Firedrake is distinguished from those projects by
its pure Python implementation, with a greater emphasis on code generation to
deliver high performance, and by its tight integration with the linear and
nonlinear solver capabilities of PETSc.

The operation of Firedrake is best illustrated with an example. We consider
the Poisson equation in three dimensions, which is simple to write and which
also acts as a model for the fluid dynamics problem which will be our target
application. The direct correspondence between the mathematics of the problem
and the user code is illustrated by highlighting the terms which which
correspond to lines of code in Listing 1, which is a minimal functional
Firedrake script which solves the Poisson equation in 3D.

We consider Poisson’s equation in a unit cube domain,
$\displaystyle\Omega=[0,1]^{3}\subset\mathbb{R}^{3}$ on a $\displaystyle
32\times 32\times 32$ tetrahedral mesh which is further refined twice to
create a hierarchy of 3 nested meshes which we will use in a multigrid solver.
We use a degree 3 continuous (CG) finite element space  We solve the Poisson
equation subhect to homogeneous (zero) Dirichlet boundary conditions (BCs):

$\left\\{\begin{aligned} -\nabla^{2}u&=f&&\text{on }\Omega,\\\ u&=0&&\text{on
}\partial\Omega,\end{aligned}\right.$ (1)

where the right hand side $\displaystyle f$ is given by:

$\displaystyle\begin{array}[]{l}f\\!\left(\begin{smallmatrix}\tiny x\\\ y\\\
z\end{smallmatrix}\right)=-\frac{\pi^{2}}{2}\\\ \hskip
35.00005pt\times\left(2\cos(\pi x)-\cos\left(\frac{\pi x}{2}\right)\right.\\\
\hskip 45.00006pt\left.-2(a^{2}+b^{2})\sin(\pi x)\tan\left(\frac{\pi
x}{4}\right)\right)\\\ \hskip 35.00005pt\times\sin(a\pi y)\sin(b\pi
z).\end{array}$ (2)

This has analytic solution

$\begin{array}[]{l}u\\!\left(\begin{smallmatrix}\tiny x\\\ y\\\
z\end{smallmatrix}\right)=\sin(\pi x)\tan\left(\frac{\pi x}{4}\right)\\\
\hskip 45.00006pt\times\sin(a\pi y)\sin(b\pi z),\end{array}$ (3)

which we use to calculate the error in our numerical solution.

Multiplying by a test function $\displaystyle v$, from the same finite element
space, and integrating by parts yields the prototypical weak form: find
$\displaystyle u\in V$ such that

$\hskip 40.00006pt\leavevmode\hbox to116.99pt{\vbox
to13.29pt{\pgfpicture\makeatletter\hbox{\hskip 58.49666pt\lower-4.14444pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{0.85,0.7,0.85}\pgfsys@color@rgb@fill{0.85}{0.7}{0.85}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} {\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{0.85,0.7,0.85}\pgfsys@color@rgb@fill{0.85}{0.7}{0.85}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{}\pgfsys@moveto{56.29666pt}{8.94444pt}\pgfsys@lineto{-56.29666pt}{8.94444pt}\pgfsys@curveto{-57.40125pt}{8.94444pt}{-58.29666pt}{8.04903pt}{-58.29666pt}{6.94444pt}\pgfsys@lineto{-58.29666pt}{-1.94444pt}\pgfsys@curveto{-58.29666pt}{-3.04903pt}{-57.40125pt}{-3.94444pt}{-56.29666pt}{-3.94444pt}\pgfsys@lineto{56.29666pt}{-3.94444pt}\pgfsys@curveto{57.40125pt}{-3.94444pt}{58.29666pt}{-3.04903pt}{58.29666pt}{-1.94444pt}\pgfsys@lineto{58.29666pt}{6.94444pt}\pgfsys@curveto{58.29666pt}{8.04903pt}{57.40125pt}{8.94444pt}{56.29666pt}{8.94444pt}\pgfsys@closepath\pgfsys@moveto{-58.29666pt}{-3.94444pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-56.29666pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\int_{\Omega}\nabla u\cdot\nabla v\ dx=\int_{\Omega}fv\
dx$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\qquad\forall v\in V.$
(4)

The solve call at line LABEL:line:solve creates a PETSc solver object which
solves the discretized system (a linear system, in this case). PETSc needs
callback functions which assemble the left hand side matrix and residual of
the weak form. Firedrake’s compiler creates optimized C representations of
these operations from the UFL provided, and passes the compiled C callbacks to
PETSc. The PETSc solver is fully programmable using the supplied dictionary of
parameters, and this can include further delegation of preconditioner stages
back to Firedrake. For example, line LABEL:line:pcpatch of Listing 2 delegates
the preconditioning of the multigrid smoother back to the PatchPC
preconditioner presented in [6]. Optimal preconditioners for PDE problems
usually depend on discretization-specific information. The tight coupling of
Firedrake and PETSc enables this mathematical fact to be easily translated
into executable code: PETSc controls the solving process and calls back to
Firedrake whenever discretization-specific code is required. The Firedrake
preconditioners can themselves contain (preconditioned) linear solves for
subsystems and will call back to PETSc to execute these.

A notable absence from listing 1 is any explicit parallel code: there are no
calls to MPI or threading directives. In fact, none are needed. The script
presented will run in parallel without modification, simply by running it
under MPI. The user code specifies the mathematical problem to be solved,
while Firedrake decides how to execute it on the parallel architecture.

⬇

1from firedrake import *

2

3N = 32 #*\label{line:meshresolution}*

4

5hierarchy = MeshHierarchy(mesh, 2) #*\label{line:meshhierarchy}*

6mesh = hierarchy[-1]

7

8 #*\label{line:degree}*

9

10u = TrialFunction(V)

11v = TestFunction(V)

12

13

14

15

16x, y, z = SpatialCoordinate(mesh)

17a = Constant(1)

18b = Constant(2)

19

20

21

22

23

24

25

26

27

28solver_parameters = {...} #*\label{line:solverparam}*

29u = Function(V)

30solve( , u, bcs=bcs, solver_parameters=solver_parameters)
#*\label{line:solve}*

Listing 1: Simple Firedrake script for solving the Poisson equation in 3D,
Listing 2 and Listing 3 give example solver parameters that can be used in
line LABEL:line:solverparam.

⬇

1solver_parameters = {

2 "ksp_type": "preonly",

3 "pc_type": "mg",

4 "pc_mg_log": None,

5 "pc_mg_type": "full",

6 "mg_levels": {

7 "ksp_type": "chebyshev",

8 "ksp_max_it": 2,

9 "ksp_norm_type":"unpreconditioned",

10 "ksp_convergence_test": "skip",

11 "pc_type": "python",

12 "pc_python_type": #*\color{blue}{"firedrake.PatchPC"}\label{line:pcpatch}*,

13 "patch_pc_patch": {

14 "construct_type": "star",

15 "construct_dim": 0

16 }

17 },

18 "mg_coarse_pc_type": "lu"

19}

Listing 2: PETSc solver options for the Poisson problem. These options use a
full multigrid method smoothed with 2 Chebyshev iterations at each level. The
preconditioner of the smoother is delegated back to Firedrake’s patch
preconditioner, and the coarse problem is solved with a direct (LU) solver.

## The performance portability challenge

Our objective will be to make effective use of different supercomputers, with
important architecture differences, to produce a numerical solution to the
steady incompressible Navier--Stokes equations, core equations of fluid
dynamics. The scenario we will simulate is the three-dimensional lid-driven
cavity (LDC) at Reynolds number 1000, a solution of which is illustrated in
Figure 1. This and similar scenarios have previously been solved in Firedrake
on up to one third of ARCHER, the previous UK national supercomputer [7, 8].
However, that machine has now been decommissioned so we find ourselves in the
familiar position of needing to move to a new machine.

Figure 1: Visualization of the velocity field in a Reynolds number 1000 lid-
driven cavity Navier--Stokes simulation.

The machines available to us are Isambard, an ARM-based supercomputer which
forms part of the UK tier-2 supercomputer offering, and ARCHER2, the new UK
national supercomputer which is based on AMD x86_64 processors. At the time of
writing, only a small portion of the latter machine has been commissioned.
Table I shows the characteristics of each compute node of these three
machines. In comparison with ARCHER, an Isambard node has nearly three times
the available memory bandwidth, and three times the core count, but each core
is significantly less powerful, so there is about 50% more bandwidth per
floating point operation (FLOP). Conversely, an ARCHER2 node has twice the
core count of an Isambard node, and these cores are much more powerful. The
consequence is that there is about three times less bandwidth per FLOP on
ARCHER2 than on Isambard. This has important consequences for how to solve our
equations, as we shall see.

Machine | CPUs | Peak DP FLOPs | Mem.BW | Architecture
---|---|---|---|---
ARCHER | 24 | 0.56TFLOPs | 119GiB/s | Intel x86_64
Isambard | 64 | 1.08TFLOPs | 318GiB/s | Aarch64
ARCHER2 | 128 | 4.61TFLOPs | 381GiB/s | AMD x86_64
TABLE I: HPC facilities used to run these demonstrations. All values are
reported _per node_. ARCHER is the previous generation of HPC; Isambard has
very high memory bandwidth (Mem. BW); ARCHER2 has very high FLOP throughput.

## Solving Poisson on Isambard

Solving a steady Navier--Stokes problem at Reynolds number 1000 is a somewhat
involved process. In order to ensure convergence of the nonlinear Newton
iteration, a sequence of problems at increasing Reynolds number must be
solved, and each of those problems is a coupled set of equations for velocity
and pressure, which must be discretized using a stable pair of finite
elements. The discretized nonlinear problem is solved using Newton’s method,
which results in a sequence of linear systems to solve. Each linear system is
solved by breaking it up into separate subproblems for the velocity and
pressure. In the augmented Lagrangian approach, a term is added to the
equations that renders the pressure subproblem straightforward, at the cost of
making the velocity solve more difficult. This velocity solve is conducted
using a multigrid algorithm. Before executing this rather complicated solver,
it is wise to consider the solution of the Poisson system with multigrid
first. By investigating the performance of multigrid for Poisson, we can
understand the consequences of the architecture in the simplest possible
situation, and then use this knowledge in the full Navier--Stokes system.

We therefore adapt the script for the Poisson system given in Listings 1 and 2
as our initial test problem. We will want to solve on a fine mesh in order to
correctly resolve all features of the solution, but this creates a very large,
globally coupled system on which iterative solves will converge very slowly.
Geometric multigrid algorithms overcome this by restricting the problem
residual (misfit) to ever-coarser meshes, and using iterative smoothers at
each level to reduce local error over ever-larger scales. On the coarsest mesh
the problem is small enough that a direct solver can be employed (such as a
Cholesky factorization of the linear system), which would be too expensive on
the fine mesh. The coarse corrections are transferred to the finer meshes,
again smoothing the local error, until we reach the finest mesh. This is
repeated a handful of times until the fully-resolved solution is recovered.

Figure 2: Strong scaling of the Poisson problem on Isambard with 3 multigrid
levels. The figure on the left shows the total runtime, as well as the
breakdown into various component processes. The right panel shows a breakdown
of the (dominant) multigrid preconditioner time into coarse solve and the two
multigrid smoothing laters. The numbers on the lowest line of each graph show
the DOF count per processor for the fine mesh. The coarse solver scales very
poorly, and an alternative approach is required. Figure 3: Solving the test
Poisson problem with 4 multigrid levels. The coarse solve no longer dominates
on 64 cores, but it scales disastrously. Figure 4: Solving the test Poisson
problem with 4 multigrid levels and telescoping the coarse grid achieves far
better strong scaling.

Figure 2 shows the time taken by the Poisson problem as it is strongly scaled
on Isambard, from one to four 64-processor nodes. In a strong scaling problem,
the total problem size (about 57 million DOFs in this case) is held constant
while the number of processors employed is increased. Eventually, the non-
parallelizable overheads in the simulation code and the similarly limited
network bandwidth or latency will dominate the calculation, with the
consequence that the code stops speeding up as the number of processors
increases. A simple model for this effect is Amdahl’s law, which models
execution time as:

$t=s+p/n$ (5)

where $\displaystyle s$ is the time taken by scalar, nonparallelizable work
$\displaystyle p$ is the time taken by parellelizable work, and $\displaystyle
n$ is the number of processors. For an optimal multigrid solver,
$\displaystyle p$ is expected to be linear in the total DOFs in the problem,
so the appropriate unit of strong scaling is the number of DOFs per processor.
This quantity will be used to measure strong scaling throughout this paper.

If we turn our attention first to the left plot of Figure 2, the three
subcomponents of the solve shown each demonstrate their own performance story.
The red line shows the amount of time taken by sparse matrix-vector
multiplications in the multigrid stages in PETSc. These are, unsurprisingly,
very fast and scale exceptionally well. The green line shows the cost of
assembling the matrices and residuals used in the various solver stages. This
is the stage of the process at which the Python layer, and the execution of
code generated in C come together. However, the execution of Python code is
duplicated across all processors, so were there a significant performance
overhead associated with the use of Python, it would show up as a strong
scaling problem here. In fact, the strong scaling of this component of the
solve time is excellent.

The same observation cannot, regrettably, be made about scaling of the
multigrid solution time. This dominates the solve time, which is expected, but
as configured completely fails to scale. The right plot in Figure 2 shows the
decomposition of this time into two multigrid smoothers on different meshes,
and the coarse solve. The smoothers show very good scaling with the finer
level 2 taking more time and the coarser level 1 scaling less well. It should
be remembered that on 512 processes, multigrid level 1 is down to around 14000
DOFs per processor (=111458/8).

The elephant in the room, however, is that the coarse solve dominates and does
not scale at all. Our first response to this is to note that the idea of
multigrid is to coarsen until the coarse solve is small and inexpensive, so we
should coarsen further. In Firedrake this can be achieved by changing
precisely two numbers in Listing 1. On line LABEL:line:meshresolution, 32 is
changed to 16 to halve the resolution of the coarse problem, and on line
LABEL:line:meshhierarchy, 2 is changed to 3 to introduce an additional
refinement (and hence multigrid level). Figure 3 shows the strong scaling
behaviour of the modified system. The significant change, other than the
additional multigrid level, is that the cost of the coarse solve on 64
processors has dropped from about 70 seconds to about 4 seconds, which is not
unexpected since the problem will have 8 times fewer degrees of freedom.
However, the scaling behaviour is still awful, and by 512 processors the
coarse solve takes 20 seconds and once more dominates the solve time.

This situation can be rectified using another algorithmic change. The problem
here is not that our coarse solve is too large, but rather that it
parallelizes poorly. This can be remedied by _telescoping_ the coarse solve:
duplicating the problem on each node and solving them independently [9]. This
can be achieved by changing the PETSc solver options. The red line in Listing
2 is replaced by the green lines in Listing 3. The results can be seen in
figure 4, which shows very good strong scaling of the whole solver out to
111458 DOFs per core. The tight integration between PETSc and Firedrake
enables sophisticated and effective solver strategies to be applied without
costly recoding of solver algorithms.

⬇

1solver_parameters = {

2 "ksp_type": "preonly",

3 "pc_type": "mg",

4 "pc_mg_log": None,

5 "pc_mg_type": "full",

6 "mg_levels": {

7 "ksp_type": "chebyshev",

8 "ksp_max_it": 2,

9 "ksp_norm_type":"unpreconditioned",

10 "ksp_convergence_test": "skip",

11 "pc_type": "python",

12 "pc_python_type": #*\color{blue}"firedrake.PatchPC"*,

13 "patch_pc_patch": {

14 "construct_type": "star",

15 "construct_dim": 0

16 }

17 },

18+ "mg_coarse": {

19+ "pc_type": "python",

20+ "pc_python_type": "firedrake.AssembledPC",

21+ "assembled": {

22+ "mat_type": "aij",

23+ "pc_type": "telescope",

24+ "pc_telescope_reduction_factor":

25+ args.telescope_factor,

26+ "pc_telescope_subcomm_type":

27+ "contiguous",

28+ "telescope_pc_type": "lu"

29+ }

30+ }

31}

Listing 3: Telescoped full multigrid PETSc solver options. This causes the
coarse problem to be solved in duplicate on each node. args.telescope_factor
is a variable that is set to yield one coarse solve per node.

## Solving Navier–Stokes on Isambard

Discretization | FE Pair | Expected | Stabilization | Refinement | Smoother
---|---|---|---|---|---
| velocity -- pressure | Convergence |  |  |
Extended Bernardi-Raugel | $\displaystyle[\mathbb{P}_{1}\oplus\mathbb{B}_{3}^{\text{F}}]^{3}-\mathbb{P}_{0}$ | 1--2 | SUPG | uniform |  additive star patches with GMRES
Scott--Vogelius | $\displaystyle[\mathbb{P}_{3}]^{2}-\mathbb{P}_{2}^{\text{disc}}$ | 4 | interior penalty | barycentric |  additive _macro_ star patches with GMRES
TABLE II: The two stable finite element (FE) pairings (for the velocity and
pressure fields) for the Navier--Stokes equations used in our simulations. The
Scott--Vogelius pair gives faster convergence at the cost of more DOFs per
element and a more expensive multigrid smoother. Expected convergence refers
to the order of convergence of the velocity in the $\displaystyle
L^{2}(\Omega)$ norm.

Having fixed the scaling of the Poisson solver, we now turn our attention to
the Navier--Stokes equations. We follow the Reynolds-robust solver approach
given in [7] and [8]. We therefore only sketch the solver algorithm here;
however, the full code executed to conduct the experiments in this paper can
be found in [10].

Table II summarizes the two discretizations we consider for Navier--Stokes.
Initially we solve using a low-order extended Bernardi-Raugel (EBR) pairing
($\displaystyle[\mathbb{P}_{1}\oplus\mathbb{B}_{3}^{F}]^{3}$\--$\displaystyle\mathbb{P}_{0}$),
which is a cheap and stable pairing with few degrees of freedom (DOFs) per
cell. An alternative is to use the Scott--Vogelius (SV) pairing
$\displaystyle[\mathbb{P}_{3}]^{2}$\--$\displaystyle\mathbb{P}_{2}^{\text{disc}}$.
This converges at higher order, and exactly enforces the incompressiblity
constraint, which is crucial at high Reynolds numbers [11]. Table II also
highlights that using an alternative discretization requires changing many
other components, including the stabilization used, the mesh hierarchy
employed, and the multigrid relaxation on each level.

Figure 5: A sequence of strong scaling results for the lid-driven cavity on
Isambard. In each case, reasonably good strong scaling is observed down to
tens of thousands of DOFs per processor. The absence of ideal scaling lines
for the highest Reynolds numbers on a single node is caused by the absence of
a single processor run for these cases: it is difficult to justify occupying a
whole node for many hours to solve a problem in serial that would take minutes
to run in parallel.

With the knowledge gained from solving Poisson we were able to adjust the
PETSc solver options for the Navier--Stokes LDC application and successfully
utilize more than a third of Isambard. Strong scaling results for three
different problem sizes (5, 40 and 320 million DOFs) with the EBR
discretization are shown in Figure 5. We solve the problem up to Reynolds
number 1000 using numerical continuation; the scaling for intermediate
solutions are also plotted, but the top green line with star markers
corresponds to the most numerically challenging solution to obtain. The figure
shows consistent strong scaling up to somewhere between $\displaystyle 40\
000$ and $\displaystyle 80\ 000$ DOFs per core on all three problem sizes.

### Solving Navier-Stokes on ARCHER2

At the time of writing the ARCHER2 HPC facility is not fully available and the
existing portion of the machine is still in the evaluation phase. As a result,
we are only currently able to report single node performance on this machine.
This still allows us to perform an ‘‘apples to apples’’ comparison of Isambard
and ARCHER2 single node performance, which is shown for the LDC Navier--Stokes
problem in Figure 6.

Figure 6: Single node strong scaling results: Isambard low order (EBR, left),
ARCHER2 low order (EBR, middle), and ARCHER2 high order (SV, bottom)
discretizations

In the left and middle plot, identical problems (size, discretization, and
solver options) are solved on Isambard and ARCHER2 respectively, and we see
two distinct behaviors. On Isambard, in the left plot, we observe good strong
scaling performance. In particular the top green line with star markers
(corresponding to the highest Reynolds number) scales well. At lower Reynolds
number the scaling is not as good, but this is a less challenging computation
and hence a less interesting case. The J-shaped runtime curve is indicative of
Amdahl’s law in action: as the number of processors increases the non-
paralellizable workload smoothly takes an increasing proportion of runtime. On
ARCHER2, in the middle graph, there is an abrupt halt to scaling and a
distinct ‘‘hockey stick’’ shape. Beyond approximately 32 cores the simulation
stops strong scaling, this suggest that we have run out of another resource;
in this case, there is not enough memory bandwidth to continue strong scaling.
This means that our seemingly quite good strong scaling curves are hiding our
failure to effectively utilize the abundance of FLOPs on an ARCHER2 node.

When a finite element simulation is limited by memory bandwidth, one way to
improve the strong scaling is to increase the degree of the finite element
used in the discretization, which increases the arithmetic intensity of the
simulation, pushing us back towards being compute bound. For the Poisson
benchmark problem it is straightforward to change polynomial degree, for
example by changing line LABEL:line:degree of Listing 1 from 3 to 4 (or
higher). For the Navier--Stokes equations, the situation is different. The
elements used for velocity and pressure must satisfy a compatibility condition
for stability and cannot be chosen independently. To obtain a more accurate
solution we change from the EBR finite element pair to the Scott--Vogelius
pair, and make the necessary modifications to the stabilization, mesh
hierarchy, and multigrid relaxation, as outlined in Table II.

Figure 6 shows the different scaling behaviour for the Scott--Vogelius
($\displaystyle[\mathbb{P}_{3}]^{3}$\--$\displaystyle\mathbb{P}_{2}^{\text{disc}}$)
discretization in the right graph. The problem size has to be reduced to fit
on a single node. One might initially observe that, whilst the strong scaling
has improved, the total number of DOFs for the problem has significantly
decreased (5 million down to 1 million), but the runtime has not changed.
However, Figure 7 demonstrates the utility of this new discretization: the
more expensive per DOF algorithm obtains a significantly more accurate
solution.

Figure 7: Norm of the difference between the computed solution and the
analytic solution for the pressure field at the different Reynolds numbers
used in the continuation steps for the LDC problem.

## CONCLUSION

Making effective use of whichever supercomputing resources are available to a
scientist or engineer at any given time demands that simulation algorithms,
and not just their implementations, are adapted to the available hardware. The
changes presented here in porting to Isambard and ARCHER2 span the
discretization used, choice and configuration of preconditioner, stabilization
scheme, and refinement algorithm. These changes were required for just one
application in fluid dynamics to be effectively ported between three machines
with relatively conventional CPU-based architectures. The immense diversity of
continuum mechanics simulation challenges across science and engineering and
their ever-changing nature means that the space of algorithmic combinations
that will be required by different users is essentially unbounded. Further,
the amount of developer time available to achieve acceptable performance for a
given application is typically severely limited.

Here we have demonstrated that, using Firedrake, it is possible to achieve
this combination of short, expressive code that can be easily modified, and
scalable parallel performance. The key to this is the tight integration of
code generation for the discretized operators with the composable,
programmable solvers and preconditioners provided by PETSc. Using Python as
the user-level language facilitates the expressivity of the user code, while
code generation and PETSc’s callback mechanisms ensure that the expensive
inner loops are optimized, compiled code. A further advantage of this
composable approach is that developments created for one application become
immediately available to other applications without the need for costly re-
implementation.

Firedrake has already been applied in many applications resulting in hundreds
of papers, but there are always limits to its capabilities. Among the many
directions in which it could be extended, support for adaptive mesh
refinement, local polynomial order refinement, and support for GPUs are
developments currently under consideration that would be particularly useful
in high performance computing contexts.

## ACKNOWLEDGMENTS

This work used the Isambard UK National Tier-2 HPC Service
(http://gw4.ac.uk/isambard/) operated by GW4 and the UK Met Office, and funded
by the Engineering and Physical Sciences Research Council (EPSRC)
(EP/P020224/1). This work used the ARCHER2 UK National Supercomputing Service
(https://www.archer2.ac.uk). This work was supported by the United Kingdom
Research and Innovation (UKRI) ExCALIBUR program (EPSRC grant EP/V001493/1).

## References

  * [1] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, and H. Zhang, ‘‘PETSc users manual,’’ Argonne National Laboratory, Tech. Rep. ANL-95/11 - Revision 3.15, 2021. [Online]. Available: https://www.mcs.anl.gov/petsc
  * [2] 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_ , vol. 43, no. 3, pp. 1--27, 2016.
  * [3] M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells, ‘‘Unified form language: A domain-specific language for weak formulations of partial differential equations,’’ _ACM Transactions on Mathematical Software_ , vol. 40, no. 2, pp. 9:1--9:37, 2014.
  * [4] A. Logg, K.-A. Mardal, and G. Wells, _Automated solution of differential equations by the finite element method: The FEniCS book_. Springer Science & Business Media, 2012, vol. 84.
  * [5] P. Bastian, M. Blatt, A. Dedner, N.-A. Dreier, C. Engwer, R. Fritze, C. Gräser, C. Grüninger, D. Kempf, R. Klöfkorn, M. Ohlberger, and O. Sander, ‘‘The dune framework: Basic concepts and recent developments,’’ _Computers & Mathematics with Applications_, vol. 81, pp. 75--112, 2021, development and Application of Open-source Software for Problems with Numerical PDEs.
  * [6] 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_ , 2021, in press.
  * [7] 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,’’ _SMAI Journal of Computational Mathematics_ , vol. 7, pp. 75--96, 2021.
  * [8] 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_ , vol. 41, no. 5, pp. A3073--A3096, 2019.
  * [9] D. A. May, P. Sanan, K. Rupp, M. G. Knepley, and B. F. Smith, ‘‘Extreme-scale multigrid components within PETSc,’’ in _Proceedings of the Platform for Advanced Scientific Computing Conference_ , 2016.
  * [10] ‘‘Software used in ’Code generation for productive portable scalable finite element simulation in Firedrake’,’’ Apr 2021. [Online]. Available: https://doi.org/10.5281/zenodo.4694370
  * [11] 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_ , vol. 59, no. 3, pp. 492--544, 2017.

|  Jack D. Betteridge, is a research software engineer at Imperial College
London. His research interests lie in code generation applications for
numerically solving PDEs on HPC. His doctorate is from the University of Bath.
Contact him at<EMAIL_ADDRESS>
---|---
|  Patrick E. Farrell, is an associate professor in the numerical analysis
group at the University of Oxford. His research interests are in fast
preconditioners for PDEs, and in bifurcation analysis of nonlinear problems.
He was jointly awarded the 2015 Wilkinson Prize for Numerical Software for the
development of the dolfin-adjoint automated inverse simulation system, a 2015
Leslie Fox prize in numerical analysis, and the 2021 Charles Broyden prize in
optimization. He holds a doctorate from Imperial College London. Contact him
at<EMAIL_ADDRESS>
---|---
|  David A. Ham, is a reader in computational mathematics at Imperial College
London. His research interests are in the development of high level
abstractions and code generation technology for PDEs. He was jointly awarded
the 2015 Wilkinson Prize for Numerical Software. He is chief-executive editor
of the journal Geoscientific Model Development. He holds a doctorate from
Delft University of Technology. Contact him at<EMAIL_ADDRESS>
---|---