# Effective codescent morphisms of $n$-quasigroups and $n$-loops

Dali Zangurashvili

###### Abstract.

Effective codescent morphisms of $n$-quasigroups and of $n$-loops are
characterized. To this end, it is proved that, for any $n\geq 1$, every
codescent morphism of $n$-quasigroups (resp. $n$-loops) is effective. This
statement generalizes our earlier results on qusigroups and loops. Moreover,
it is shown that the elements of the amalgamated free products of
$n$-quasigroups (resp. $n$-loops) have unique normal forms, and that the
varieties of $n$-quasigroups and $n$-loops satisfy the strong amalgamation
property. The latter two statements generalize the corresponding old results
on quasigroups and loops by Evans.

Key words and phrases: effective codescent morphism; $n$-quasigroup; $n$-loop;
normal form of an element of the amalgamated free product; strong amalgamation
property.

2020 Mathematics Subject Classification: 18E50, 18C20, 20N15, 20N05, 08B25,
68Q42.

##  1\. Introduction

The notion of an effective (co)descent morphism in a category is one of the
main notions of the Grothendieck descent theory. The problem of characterizing
effective codescent morphisms in varieties of universal algebras was posed by
Janelidze. This problem was studied in a number of works [9], [6], [13], [14],
[15], [10], [16], [11]. In [15], we proved that every codescent morphism of
quasigroups (resp. loops) is effective. In the present paper, we generalize
this result to the case of $n$-quasigroups and $n$-loops, for arbitrary $n\geq
1$. This generalization together with the criterion for a monomorphism of a
regular category with pushouts and the strong amalgamation property to be a
codescent morphism, found in [13], gives the characterization of effective
codescent morphisms of $n$-quasigroups.

We also show that, for any $n\geq 1$, the elements of the amalgamated free
products of $n$-quasigroups (resp. $n$-loops) have unique normal forms and
that the varieties of $n$-quasigroups and $n$-loops satisfy the strong
amalgamation property. These statements generalize the corresponding old
results on quasigroups and loops by Evans.

The author gratefully acknowledges the financial support from Shota Rustaveli
National Science Foundation of Georgia (FR-22-4923).

##  2\. The convergent representations of the varieties of n-quasigroups and
$n$-loops

We follow the notation used in the treatise [2] by Belousov. Namely, for
natural numbers $i$ and $j$, the symbol $x_{i}^{j}$ denotes the sequence
$x_{i},x_{i+1},...,x_{j}$ if $i<j$; it denotes the symbol $x_{i}$ if $i=j$,
and denotes the empty sequence if $i>j$. The symbol $\stackrel{{\scriptstyle
m}}{{e}}$ denotes the sequence $e,e,...,e$ ($m$ times) if $m\geq 1$, and
denotes the empty sequence if $m=0$.

Let $n$ be a natural number. Recall [2] that an $n$-quasigroup is defined as a
set $Q$ equipped with an $n$-ary operation $f$ such that, for arbitrary
$a_{1},a_{2},...,a_{n}\in Q$ and $i$ $(1\leq i\leq n)$, there is a unique
element $b\in Q$ with

$f(a_{1}^{i-1},b,a_{i+1}^{n})=a_{i}.$ (2.1)

An $n$-quasigroup $Q$ is called an $n$-loop [2] if there is an element $e\in
Q$ such that, for any $i$ $(1\leq i\leq n)$, the identity

$f(\stackrel{{\scriptstyle i-1}}{{e}},x,\stackrel{{\scriptstyle n-i}}{{e}})=x$
(2.2)

is satisfied in $Q$.

One can easily verify that the category of $n$-quasigroups (with morphisms
being mappings preserving $f$) is the variety of universal algebras where the
signature consists of the $n$-ary operation symbols $f,g_{i}$, while the
identities are

$f(x_{1}^{i-1},g_{i}(x_{1}^{n}),x_{i+1}^{n})=x_{i}.$ (2.3)

and

$g_{i}(x_{1}^{i-1},f(x_{1}^{n}),x_{i+1}^{n})=x_{i}$ (2.4)

$(1\leq i\leq n)$. Similarly, the category of $n$-loops is the variety of
universal algebras where the signature, in addition to the $n$-ary operation
symbols $f,g_{i}$, contains one nullary operation $e$, while the identities
are (2.2)-(2.4) $(1\leq i\leq n)$.

Obviously, $2$-quasigroups are precisely quasigroups, and $2$-loops are
precisely loops. Observe also that $1$-quasigroups are precisely sets equipped
with permutations, while $1$-loops are just sets.

Let $\mathcal{C}$ be a category with pushouts, and $p:B\rightarrow E$ be its
morphism. The morphism $p$ induces the change-of-cobase functor

$p_{*}:B/\mathcal{C}\rightarrow E/\mathcal{C}$

between the coslice categories; it sends a morphism $\varphi:B\rightarrow C$
to the pushout of $\varphi$ along $p$. As is well-known, this functor has a
right adjoint. Recall that $p$ is called a codescent morphism (resp. an
effective codescent morphism) if the functor $p_{*}$ is precomonadic (resp.
comonadic) [5].

In [13], we gave the necessary and sufficient condition for a morphism of a
variety of universal algebras with the strong amalgamation property to be a
codescent morphism. In [15], we gave the sufficient condition for any
codescent morphism of a variety of universal algebras with the amalgamation
property to be effective. Applying this sufficient condition, in the joint
paper with Samsonadze [10], we related the problem of characterizing effective
codescent morphisms in varieties of universal algebras to the notion of
confluency that is one of the central notions in term rewriting theory. Below
we recall some definitions from this theory (for more details, we refer the
reader to the treatise [1] by Baader and Nipkow).

Let $\mathfrak{F}$ be a signature, i.e., a set of operation symbols equipped
with arities. Let $X$ be a countable set with $\mathfrak{F}\cap
X=\varnothing$. We assume that the reader is familiar with the notions of a
$\mathfrak{F}$-term (or simply a term) over $X$, the size of a term, a
position in a term, and the subterm of a term $t$ at a position $p$. We denote
the latter subterm by $t|_{p}$, and denote the term obtained from $t$ by
replacing the subterm $t|_{p}$ by a term $t^{\prime}$ by the symbol
$t[t^{\prime}]_{p}$.

As usual, we identify positions in a term with nodes in the corresponding
tree, and use the well-known rule of enumerating them with strings of natural
numbers. For instance, the tree and the enumeration of nodes for the term
$t=f(x_{1},x_{2},g_{3}(x_{1},x_{2},x_{3},x_{4}),x_{4})$ is the following:

$\textstyle{f\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{1\,\,\,\,x_{1}}$$\textstyle{2\,\,\,\,x_{2}}$$\textstyle{g_{3}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{4\,\,\,\,x_{4}}$$\textstyle{31\,\,\,\,x_{1}}$$\textstyle{32\,\,\,\,x_{2}}$$\textstyle{33\,\,\,\,x_{3}}$$\textstyle{34\,\,\,\,x_{4}}$

If, for instance, $p=(3)$, then $t|_{p}=g_{3}(x_{1},x_{2},x_{3},x_{4})$. The
node $f$ corresponds to the root position.

Let $T(\mathfrak{F},X)$ be the set of all terms over $X$, and $\sigma$ be a
substitution (i.e., a mapping $\sigma:X\rightarrow T(\mathfrak{F},X)$ such
that the set $\\{x\mid\sigma(x)\neq x\\}$ is finite). Since
$T(\mathfrak{F},X)$ has the structure of an $\mathfrak{F}$-algebra and is free
over the set $X$, there is a unique $\mathfrak{F}$-homomorphism
$\widehat{\sigma}:T(\mathfrak{F},X)\rightarrow T(\mathfrak{F},X)$ such that
its restriction on $X$ is $\sigma$.

An oriented identity is a formal expression of the form $l=r$, where $l$ and
$r$ are terms (the order of the terms matters). In this paper, when no
confusion might arise, ‘identity’ means an oriented identity.

We denote the set of variables that occur in at least one of terms
$t_{1},t_{2},...,t_{n}$ by $Var(t_{1},t_{2},...,t_{n})$.

Let $\Sigma$ be a set of identities. A pair $(\mathfrak{F},\Sigma)$ is called
a term rewriting system if $l$ is not a variable and $Var(r)\subseteq Var(l)$,
for any identity $l=r$ from $\Sigma$.

Let $(\mathfrak{F},\Sigma)$ be a term rewriting system. One introduces the
following binary relation $\rightarrow$ on the set $T(\mathfrak{F},X)$ of
terms: $t\rightarrow t^{\prime}$ if the condition (C) below is satisfied.

(C) there exists an identity $l=r$ from $\Sigma$, a substitution $\sigma$ and
a position $p$ of $t$ such that $t|_{p}=\widehat{\sigma}(l)$ and
$t^{\prime}=t[\widehat{\sigma}(r)]_{p}$. Roughly speaking, $t\rightarrow
t^{\prime}$ means that the term $t^{\prime}$ can be obtained from $t$ by
replacing some subterm, and this replacement is compliant with some identity
from $\Sigma$. At that, we are permitted to use identities only in one
direction – from left to right.

We use the symbol $\xrightarrow{*}$ for the reflexive transitive closure of
$\rightarrow$.

A term rewriting system $(\mathfrak{F},\Sigma)$ (or simply $\Sigma$) is called
terminating if there is no infinite sequence of terms

$t_{1}\rightarrow t_{2}\rightarrow...$

A term rewriting system $(\mathfrak{F},\Sigma)$ (or simply $\Sigma$) is called
confluent if, for any terms $t,t_{1},t_{2}$ with
$t_{1}\xleftarrow{*}t\xrightarrow{*}t_{2}$, there is a term $t^{\prime}$ such
that $t_{1}\xrightarrow{*}t^{\prime}\xleftarrow{*}t_{2}$ (in that case, the
terms $t_{1}$ and $t_{2}$ are called joinable):

$\textstyle{t\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{*}$$\scriptstyle{*}$$\textstyle{t_{1}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{*}$$\textstyle{t_{2}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{*}$$\textstyle{t^{\prime}}$

Further, recall that a unifier of terms $t$ and $t^{\prime}$ is defined as a
substitution $\sigma$ such that
$\widehat{\sigma}(t)=\widehat{\sigma}(t^{\prime})$. A unifier $\sigma$ is
called most general if any unifier of these terms can be obtained by composing
$\widehat{\sigma}$ with some substitution. If two terms have at least one
unifier, then they have also the most general unifier; it is unique up to
renaming (i.e. an injective substitution $\varrho$ with $\varrho(X)\subseteq
X$) [1].

Let $\iota_{1}:l_{1}=r_{1}$ and $\iota_{2}:l_{2}=r_{2}$ be identities from
$\Sigma$. Let us rename their variables so that

$Var(l_{1},r_{1})\cap Var(l_{2},r_{2})=\varnothing.$

Let $p$ be a position of $l_{1}$ such that $l_{1}|_{p}$ is not a variable and
the terms $l_{1}|_{p}$ and $l_{2}$ unify (in that case, one says that the term
$l_{1}$ overlaps the term $l_{2}$ at the position $p$). Let $\sigma$ be the
most general unifier of $l_{1}|_{p}$ and $l_{2}$. The pair of terms

$(\widehat{\sigma}(r_{1}),\widehat{\sigma}(l_{1})[\widehat{\sigma}(r_{2})]_{p})$

is called the critical pair determined by the identities $\iota_{1}$ and
$\iota_{2}$ (at the position $p$ of $l_{1}$).

Observe that the term $\widehat{\sigma}(l_{1})$ is in the relation
$\rightarrow$ with both terms $\widehat{\sigma}(r_{1})$ and
$\widehat{\sigma}(l_{1})[\widehat{\sigma}(r_{2})]_{p}$ participating in the
critical pair:

$\textstyle{\widehat{\sigma}(l_{1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\widehat{\sigma}(r_{1})}$$\textstyle{\widehat{\sigma}(l_{1})[\widehat{\sigma}(r_{2})]_{p}}$
(2.5)

Moreover, the latter two terms are obviously joinable if $p$ is a root
position of $\iota_{1}$ and $\iota_{2}$ is a renamed copy of $\iota_{1}$.

######  Theorem 2.1.

[1] Let $(\mathfrak{F},\Sigma)$ be terminating. Then it is confluent if and
only if all critical pairs arisen from the identities of $\Sigma$ (including
the critical pairs determined by identities and their renamed copies) are
joinable.

######  Lemma 2.2.

The system of identities (2.3)-(2.4) is confluent if and only if $n=1$.

###### Proof.

First observe that for any $i$ and $j$ ($1\leq i,j\leq n$), the left-hand term
$l_{1}$ of identity (2.4) overlaps the left-hand term $l_{2}$ of the identity

$f(y_{1}^{j-1},g_{j}(y_{1}^{n}),y_{j+1}^{n})=y_{j}$ (2.6)

at the position $p=(i)$. The most general unifier $\sigma$ of $l_{1}|_{p}$ and
$l_{2}$ is given as follows: $x_{j}\mapsto g_{j}(y_{1},y_{2},...,y_{n})$ and
$x_{k}\mapsto y_{k}$, for $k\neq j$. Assume now that $1\leq i<j\leq n$. Then
diagram (2.5) takes the form

---
$\textstyle{g_{i}(y_{1}^{i-1},f(y_{1}^{j-1},g_{j}(y_{1}^{n}),y_{j+1}^{n}),y_{i+1}^{j-1},g_{j}(y_{1}^{n}),y_{j+1}^{n})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\,\,\,\,\,\,\,\,y_{i}}$$\textstyle{\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,g_{i}(y_{1}^{i-1},y_{j},y_{i+1}^{j-1},g_{j}(y_{1}^{n}),y_{j+1}^{n})}$

The terms in the obtained critical pair obviously are not joinable with
respect to the system (2.3)-(2.4) of identities. Theorem 2.1 implies that this
system is not confluent. The ”if” part of the claim is immediate. ∎

The arguments given in the proof of Lemma 2.2. imply that, in the variety of
$n$-quasigroups with $n\geq 2$, the identity

$g_{i}(x_{1}^{i-1},x_{j},x_{i+1}^{j-1},g_{j}(x_{1}^{n}),x_{j+1}^{n})=x_{i}$
(2.7)

is satisfied if $i<j$, and the identity

$g_{i}(x_{1}^{j-1},g_{j}(x_{1}^{n}),x_{j+1}^{i-1},x_{j},x_{i+1}^{n})=x_{i}$
(2.8)

is satisfied if $i>j$.

######  Lemma 2.3.

Let $n\geq 2$. The system of identities (2.3)-(2.4) with $1\leq i,j\leq n$,
(2.7) with $1\leq i<j\leq n$ and (2.8) with $1\leq j<i\leq n$ is a confluent
representation of the variety of $n$-quasigroups.

###### Proof.

To account all critical pairs, note that the left-hand terms $t_{1}$ and
$t_{2}$ of identities (2.4) and (2.7) do not unify, for any $i<j$. Indeed, if
$\sigma$ were their unifier, then the size $k_{1}$ of the subterm $t|_{(i)}$
of the term $t=\widehat{\sigma}(t_{1})$ would be greater than the size $k_{2}$
of the subterm $t|_{(j)}$. On the other hand, $k_{1}<k_{2}$ since
$t=\widehat{\sigma}(t_{2})$, and we arrive to the contradiction.

Consider now, for example, the case where $\iota_{1}$ is identity $(2.3)$,
$p=(i)$, and $\iota_{2}$ is the renamed copy of $(2.7)$ (where $x_{i}$’s are
replaced by $y_{i}$’s; $i<j$). The most general unifier of $l_{1}|_{p}$ and
$l_{2}$ sends $x_{i}$ to $y_{j}$, $x_{j}$ to $g_{j}(y_{1},...,y_{n})$, and
sends all other $x_{k}$’s to $y_{k}$’s. It is easy to see that diagram (2.5)
takes the form

---
$\textstyle{f(y_{1}^{i-1},g_{i}(y_{1}^{i-1},y_{j},y_{i+1}^{j-1},g_{j}(y_{1}^{n}),y_{j+1}^{n}),y_{i+1}^{j-1},g_{j}(y_{1}^{n}),y_{j+1}^{n})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\,\,\,\,\,\,\,\,\,\,\,\,\,y_{j}}$$\textstyle{\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,f(y_{1}^{j-1},g_{j}(y_{1}^{n}),y_{j+1}^{n})}$

We obtain the critical pair that obviously is joinable due to identity (2.3).
Similarly, one can verify that all other critical pairs are joinable.

∎

We are going now to deal with the variety of $n$-loops ($n\geq 1$).

######  Lemma 2.4.

(i) Let $n$ be an arbitrary natural number. The system of identities
(2.2)-(2.4) $(1\leq i\leq n)$ is not confluent.

(ii) Let $n>1$. The union of the system of identities given in Lemma 2.3 with
the system of identities (2.2) ($1\leq i\leq n$) is not confluent.

###### Proof.

Let $\iota_{1}$ be identity (2.2), $\iota_{2}$ be identity (2.3), and $p$ be
the root position of $l_{1}$. The most general unifier of $l_{1}$ and $l_{2}$
sends $x_{j}$ to $e$, for $j\neq i$, and sends $x$ to
$g_{i}(e^{i-1},y_{i},e^{n-i})$. Then diagram (2.5) takes the form

$\textstyle{f(\stackrel{{\scriptstyle i-1}}{{e}},g_{i}(\stackrel{{\scriptstyle
i-1}}{{e}},x_{i},\stackrel{{\scriptstyle n-i}}{{e}}),\stackrel{{\scriptstyle
n-i}}{{e}})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{g_{i}(\stackrel{{\scriptstyle
i-1}}{{e}},x_{i},\stackrel{{\scriptstyle n-i}}{{e}})}$$\textstyle{x_{i}}$

We see that the terms in the critical pair that arises from this overlap is
not joinable with respect to the set (2.2)-(2.4). If $n>1$, these terms are
not joinable also with respect to the system of identities mentioned in the
claim (ii). ∎

Equating the terms in the critical pair arisen in the proof of Lemma 2.4, we
arrive to the identity

$g_{i}(\stackrel{{\scriptstyle i-1}}{{e}},x,\stackrel{{\scriptstyle
n-i}}{{e}})=x.$ (2.9)

Therefore, it is satisfied in the variety of $n$-loops, for any $n\geq 1$ and
any $i$ $(1\leq i\leq n$).

######  Lemma 2.5.

The system of identities (2-2)-(2.4), (2.9) $(1\leq i\leq n)$ is confluent if
and only if $n=1$. If $n>1$, then the union of of the system of identities
given in Lemma 2.4(ii) with the system of identities (2.9) ($1\leq i\leq n$)
is not confluent.

###### Proof.

For ”only if” part of the claim, let $1\leq i<j\leq n$. Let $\iota_{1}$ be
identity $(2.4)$, $p=(i)$, and $\iota_{2}$ be the identity

$f(\stackrel{{\scriptstyle j-1}}{{e}},y,\stackrel{{\scriptstyle
n-j}}{{e}})=y.$

The most general unifier of $l_{1}|_{p}$ and $l_{2}$ sends $x_{j}$ to $y$, and
$x_{k}$ to $e$, for $k\neq j$. Therefore, diagram (2.5) takes the form

$\textstyle{g_{i}(\stackrel{{\scriptstyle i-1}}{{e}},f(\stackrel{{\scriptstyle
j-1}}{{e}},y,\stackrel{{\scriptstyle n-j}}{{e}}),\stackrel{{\scriptstyle
j-i-1}}{{e}},y,\stackrel{{\scriptstyle
n-j}}{{e}})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{e}$$\textstyle{g_{i}(\stackrel{{\scriptstyle
i-1}}{{e}},y,\stackrel{{\scriptstyle j-i-1}}{{e}},y,\stackrel{{\scriptstyle
n-j}}{{e}})}$

The terms $e$ and $g_{i}(\stackrel{{\scriptstyle
i-1}}{{e}},y,\stackrel{{\scriptstyle j-i-1}}{{e}},y,\stackrel{{\scriptstyle
n-j}}{{e}})$ obviously are not joinable with respect to any system of
identities considered above. The ”if” part of the claim is trivial. ∎

In view of the proof of Lemma 2.5, consider the identity

$g_{i}(\stackrel{{\scriptstyle i-1}}{{e}},x,\stackrel{{\scriptstyle
j-i-1}}{{e}},x,\stackrel{{\scriptstyle n-j}}{{e}})=e,$ (2.10)

for any $i<j$, and the identity

$g_{i}(\stackrel{{\scriptstyle j-1}}{{e}},x,\stackrel{{\scriptstyle
i-j-1}}{{e}},x,\stackrel{{\scriptstyle n-i}}{{e}})=e,$ (2.11)

for any $j<i$. The arguments of the above-mentioned proof imply that these
identities are satisfied in the variety of $n$-loops, for any $n\geq 2$.
Similarly to Lemma 2.3, one can verify the following

######  Lemma 2.6.

Let $n\geq 2$. The system of identities (2.2)-(2.4) with $1\leq i\leq n$,
(2.7) with $1\leq i<j\leq n$, (2.8) with $1\leq j<i\leq n$, (2.9) with $1\leq
i\leq n$, (2.10) with $1\leq i<j\leq n$, and (2.11) with $1\leq j<i\leq n$ is
a confluent representation of the variety of $n$-loops.

##  3\. Effective codescent morphisms in the varieties of $n$-quasigroups and
$n$-loops

In [10], we considered the following conditions on the signature
$\mathfrak{F}$ and identities of a term rewriting system
$(\mathfrak{F},\Sigma)$:

(*) for any identity $l=r$ from $\Sigma$, no variable occurs in $r$ more often
than in $l$, and moreover, the size of $l$ is greater than the size of $r$;

(**) if the set $\mathfrak{F}_{0}$ of constants from $\mathfrak{F}$ is not
empty, then, for any non-trivial algebra $A$ from the variety determined by
$\mathfrak{F}$ and $\Sigma$, the mapping $\mathfrak{F}_{0}\rightarrow A$
sending a constant to its value in $A$ is injective;

(***) for any identity $l=r$ from $\Sigma$, any subterm $l^{\prime}$ of $l$
which is neither a variable nor a constant, we have $Var(l^{\prime})=Var(l)$.

Before continue, recall some definitions.

Let $\mathcal{V}$ be the variety of universal algebras determined by
$(\mathfrak{F},\Sigma)$. Let $I$ be a non-empty set, and
$(m_{i}:B\rightarrowtail A_{i})_{i\in I}$ be a family of
$\mathcal{V}$-algebras with an amalgamated subalgebra $B$. We assume that
$A_{i}\cap A_{j}=B$, for $i\neq j$ (and also that
$A_{i}\cap\mathfrak{F}=\emptyset$).

One can introduce the binary relation $\rightsquigarrow$ on the set of
$\mathfrak{F}$-terms over the set $\underset{i\in I}{\cup}A_{i}$ as follows:
$t\rightsquigarrow t^{\prime}$ if either the above-mentioned condition (C) is
satisfied (for $X=\underset{i\in I}{\cup}A_{i}$) or $t^{\prime}$ can be
obtained from $t$ by replacing some subterm – a subterm such that all its
variables (being elements of algebras) belong to one and the same algebra
$A_{i}$ – by the value of this subterm in $A_{i}$.

The condition (*) implies that any element $\alpha$ of the free product $A$ of
$(A_{i})_{i}\in I$ with the amalgamated subgroup $B$ can be written as an
irreducible term over the set $\underset{i\in I}{\cup}A_{i}$, i.e., a term
$\tau$ such that $\tau\rightsquigarrow\tau^{\prime}$ for no term
$\tau^{\prime}$ over the same set. The term $\tau$ is called a normal form of
the element $\alpha$ [10]. In general, an element $\alpha$ may have more than
one normal form.

Finally, recall that a variety is said to satisfy the amalgamation property
if, for any pushout

$\textstyle{A\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{f}$$\scriptstyle{g}$$\textstyle{B\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{g^{\prime}}$$\textstyle{C\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{f^{\prime}}$$\textstyle{D}$
(3.1)

with monomorphic $f$ and $g$, the homomorphisms $f^{\prime}$ and $g^{\prime}$
also are such. A variety is said to satisfy the strong amalgamation property
if it satisfies the amalgamation property and, moreover, for any pushout (3.1)
with monomorphic $f$ and $g$, one has $f^{\prime}(C)\cap
g^{\prime}(B)=f^{\prime}g(A)$ [8].

######  Theorem 3.1.

[10] Let a variety $\mathcal{V}$ of universal algebras be represented by a
confluent term rewriting system $(\mathfrak{F},\Sigma)$ that satisfies the
conditions (*)-(***). Then

(a) every codescent morphism of $\mathcal{V}$ is effective;

(b) the elements of amalgamated free products in $\mathcal{V}$ have unique
normal forms;

(c) the variety $\mathcal{V}$ satisfies the strong amalgamation property.

From Lemmas 2.2, 2.3, 2.5, 2.6 and Theorem 3.1 we obtain

######  Theorem 3.2.

Let $n$ be an arbitrary natural number.

(a) Every codescent morphism of $n$-quasigroups (resp. $n$-loops) is
effective;

(b) the elements of amalgamated free products in the variety of
$n$-quasigroups (resp. $n$-loops) have unique normal forms;

(c) the variety of $n$-quasigroups (resp. $n$-loops) satisfies the strong
amalgamation property.

Note that, for $n=2$, the claim (a) of Theorem 3.2 was first given in [15]
(see also [10]), while the claims (b) and (c) were first given in the paper
[3] by Evans (see also [7] and [15]). For $n=1$, the claim (a) is obvious
since the variety of $1$-quasigroups (being isomorphic to the category
$\mathcal{S}et^{\mathbb{Z}}$ of functors from the group $\mathbb{Z}$ of
integers, viewed as a category, to the category of sets) and the variety of
$1$-loops are topoi. The claim (c), for $1$-loops is obvious. For
$1$-quasigroups, the latter claim immediately follows from the main result of
[12].

Recall

######  Theorem 3.3.

[4] Let $\mathcal{C}$ be a category with pushouts and equalizers. A morphism p
is a codescent morphism if and only if it is a couniversal regular
monomorphism. i.e. a morphism whose any pushout is a regular monomorphism.

It is easy to observe that any monomorphism is regular in a variety with the
strong amalgamation property.

######  Theorem 3.4.

[13] Let a variety $\mathcal{V}$ of universal algebras satisfy the strong
amalgamation property. Then a monomorphism $A\rightarrowtail A^{\prime}$ of
$\mathcal{V}$ is a codescent morphism if and only if it satisfies the
congruence extension property, i.e., for any congruence $R$ on $A$, there is a
congruence $R^{\prime}$ on $A^{\prime}$ such that $R^{\prime}\cap(A\times
A)=R$.

Note that congruences mentioned in the congruence extension property in
Theorem 3.4 can be congruences with respect to any representation of the
variety (while normal forms mentioned in Theorem 3.1(c) are ones with respect
to the considered representation $(\mathfrak{F},\Sigma)$). Therefore, Theorem
3.2(a), Theorem 3.2(c) and Theorem 3.4 imply

######  Theorem 3.5.

Let $n$ be an arbitrary natural number. A monomorphism $Q\rightarrowtail
Q^{\prime}$ of $n$-quasigroups (resp. $n$-loops) is an effective codescent
morphism if and only if it satisfies the $\mathfrak{F}$-congruence extension
property, where $\mathfrak{F}=\\{f\\}\cup\\{g_{i}\mid 1\leq i\leq n\\}$ (and
all operation symbols are $n$-ary).

Note that, for $n=2$, the set of $\mathfrak{F}$-congruences on a quasigroup
$Q$ do not coincide with that of $\\{f\\}$-congruences on $Q$ [7]. The former
congruences are referred to as normal congruences in literature. Not any
monomorphism of quasigroups (resp. loops) satisfies the normal-congruence
extension property [7].

Let $n=1$. Then any monomorphism is an effective codescent morphism in both
the category of $1$-quasigroups and the category of $1$-loops since they are
topoi. Theorem 3.5 implies that, in these varieties, any monomorphism
satisfies the $\mathfrak{F}$-congruence extension property (this can be easily
seen immediately). Observe also that, not any $\\{f\\}$-congruence on a
$1$-quasigroup is an $\mathfrak{F}$-congruence (for an example, consider a
$1$-quasigroup $(\mathbb{Z},f)$ with $f(k)=k+1$, for any $k\in\mathbb{Z}$. Let
$\theta$ be the equivalence relation on $(\mathbb{Z},f)$ with the following
equivalence classes: $\\{k\\}$ ($k<1$), and $\mathbb{N}$. The relation
$\theta$ is obviously a congruence with respect to the signature $\\{f\\}$,
but is not a congruence with respect to the one
$\mathfrak{F}=\\{f,f^{-1}\\}$). Nevertheless, we have

######  Proposition 3.6.

Let $Q$ be a finite 1-quasigroup. Then any $\\{f\\}$-congruence on $Q$ is an
$\mathfrak{F}$-congruence.

###### Proof.

Let $\theta$ be an $\\{f\\}$-congruence on $Q$. Let $C_{a}$ denote the
$\theta$-class of an element $a\in Q$. One obviously has

$f(C_{a})\subseteq C_{f(a)}.$ (3.2)

Inequality (3.2) implies that we have the increasing chain of natural numbers:

$card(C_{a})\leq card(C_{f(a)})\leq card(C_{f^{2}(a)})\leq...$

Since $Q$ is finite, this chain stabilizes at some step $k\geq 0$ (assuming
that $f^{0}=id_{Q}$). Inequality (3.2) implies that, for any $m\geq k$, we
have

$f(C_{f^{m}(a)})=C_{f^{m+1}(a)}.$ (3.3)

Assume now that $k\geq 1$. Since $Q$ is finite and no distinct
$\theta$-classes have non-empty intersection, (3.3) implies that there is the
smallest $l\geq k$ with

$f(C_{f^{l}(a)})=C_{f^{l}(a)}.$

Since $l\geq 1$ and $f$ is a bijection, we arrive to the contradiction.
Therefore, $k=0$, and we have

$f(C_{a})=C_{f(a)}.$

∎

## References

  * [1] F. Baader, T. Nipkow, Term rewriting and all that, Cambridge University Press, 2006.
  * [2] V. D. Belousov, n-ary Quasigroups. Shtiintsa, Kishinev, 1972 (In Russian).
  * [3] T. Evans, On multiplicative systems defined by generators and relations. I. Normal form theorems. Proc. Camb. Philos. Soc. 47, 637–649 (1951).
  * [4] G. Janelidze, W. Tholen, How algebraic is the change-of-base functor? in: Lecture Notes in Math., 1488, Springer, Berlin, 1991, 174-186.
  * [5] G. Janelidze, W. Tholen, Facets of descent, I. Appl. Categ. Structures, 2(1994), 245-281.
  * [6] G. Janelidze, W. Tholen, Facets of descent, III. Appl. Categ. Structures, 12(5-6), 461-477(2004).
  * [7] J. Ježek, T. Kepka, Varieties of quasigroups determined by short strictly balanced identities, Czechoslovak Math. J. 29 (104)(1979), 84–96.
  * [8] E.W. Kiss, L. Márki, L. Pröhle, W. Tholen, Categorical algebraic properties. A compendium on amalgamation, congruence extension, epimorphisms, residual smallness and injectivity. Stud. Sci. Math. Hung., 18(1983), 79-141.
  * [9] B. Mesablishvili, Pure morphisms of commutative rings are effective descent morphisms for modules – a new proof, Theory Appl. Categ., 7(3)(2000), 38-42.
  * [10] G. Samsonadze, D. Zangurashvili, Effective codescent morphisms in the varieties determined by convergent term rewriting systems, Tbilisi Math. J. 9(1)(2016), 49-64.
  * [11] G. Samsonadze, D. Zangurashvili, Descent in the dual category of ternary rings, 2023 https://arxiv.org/abs/2308.04747
  * [12] D. Zangurashvili, On some categorical algebraic properties of the category of functors with the values in concrete categories, Bull. Acad. Sci. Georgian. SSR., 135(2)(1989), 2-4.
  * [13] D. Zangurashvili, The strong amalgamation property and (effective) codescent morphisms, Theory Appl. Categ. 11(20)(2003), 438-449.
  * [14] D. Zangurashvili, Effective codescent morphisms, amalgamations and factorization systems, J. Pure Appl. Algebra 209(1)(2007), 255-267.
  * [15] D. Zangurashvili, Effective codescent morphisms in some varieties of universal algebras, Appl. Categ. Structures, 22(2014), 241-252.
  * [16] D. Zangurashvili, Admissible Galois Structures on the categories dual to some varieties of universal algebras, Georgian Math. Journal, 28(4)(2021), 651-664.

Author’s address:

Dali Zangurashvili, A. Razmadze Mathematical Institute of Tbilisi State
University, 6 Alexidze Str., 0193, Georgia;

e-mail<EMAIL_ADDRESS>


# Deep neural operators can serve as accurate surrogates for shape
optimization: A case study for airfoils

*Khemraj Shukla, *Vivek Oommen, *Ahmad Peyvan
Brown University, Providence, RI, 02912

&*Michael Penwarden
University of Utah, Salt Lake City, UT, 84112
Luis Bravo, Anindya Ghoshal
Weapons and Materials Directorate, U.S. Army Research Laboratory
Aberdeen Proving Ground, MD, 21005

&Robert M. Kirby
University of Utah, Salt Lake City, UT, 84112

&George Em Karniadakis
Brown University, Providence, RI, 02912

&*Equal contribution

###### Abstract

Deep neural operators, such as DeepONets, have changed the paradigm in high-
dimensional nonlinear regression from function regression to (differential)
operator regression, paving the way for significant changes in computational
engineering applications. Here, we investigate the use of DeepONets to infer
flow fields around unseen airfoils with the aim of shape optimization, an
important design problem in aerodynamics that typically taxes computational
resources heavily. We present results which display little to no degradation
in prediction accuracy, while reducing the online optimization cost by orders
of magnitude. We consider NACA airfoils as a test case for our proposed
approach, as their shape can be easily defined by the four-digit
parametrization. We successfully optimize the constrained NACA four-digit
problem with respect to maximizing the lift-to-drag ratio and validate all
results by comparing them to a high-order CFD solver. We find that DeepONets
have low generalization error, making them ideal for generating solutions of
unseen shapes. Specifically, pressure, density, and velocity fields are
accurately inferred at a fraction of a second, hence enabling the use of
general objective functions beyond the maximization of the lift-to-drag ratio
considered in the current work.

_K_ eywords Neural operators $\cdot$ DeepONet $\cdot$ Airfoil shape
optimization $\cdot$ Navier-Stokes equations $\cdot$ Surrogate models

## Nomenclature

$m$ = | Maximum camber in percentage of the chord
---|---
$p$ = | Position of maximum camber in percentage of the chord
$t$ = | Maximum thickness of the airfoil in percentage of the chord
$\rho$ = | non-dimensional density
$u$ = | non-dimensional velocity of the fluid in x direction
$v$ = | non-dimensional velocity of the fluid in y direction
$p$ = | non-dimensional pressure
$T$ = | non-dimensional temperature
$\mathcal{G}$ = | output of the DeepONet model
$\theta_{b}$ = | set of all trainable weights and biases of the branch network
$\theta_{t}$ = | set of all trainable weights and biases of the trunk network
$N_{\phi}$ = | number of basis functions learned by the DeepONet
$\tau_{ij}$ = | Components of viscous stress tensors
$f$ = | generalized input function to the DeepONet
$\xi_{g}$ = | Geometric parameter
$c$ = | Chord length
$Ma$ = | Mach number
$\xi_{f}$ = | Flow parameter
$Re$ = | Reynolds number
$Pr$ = | Prandtl number
$\bm{\alpha}$ = | output of the branch network
$\bm{\phi}$ = | output of the trunk network
$\tau_{w}$ = | Wall Shear Stress
$L$ = | Lift
$D$ = | Drag
$\overrightarrow{n}$ = | Unit Normal
$\overrightarrow{t}$ = | Unit Tangent

##  1 Introduction

Neural networks that solve regression problems map input data to output data,
whereas neural operators map functions to functions. This recent paradigm
shift in perspective, starting with the original paper on the deep operator
network or DeepONet [1, 2], provides a new modeling capability that is very
useful in engineering – that is, the ability to replace very complex and
computational resource-taxing multiphysics systems with neural operators that
can provide functional outputs in real-time. Specifically, unlike other
physics-informed neural networks (PINNs) [3] that require optimization during
training and testing, a DeepONet does not require any optimization during
inference, hence it can be used in real-time forecasting, including design,
autonomy, control, etc. An architectural diagram of a DeepONet with the
commonly used nomenclature for its components is shown in Figure 1. DeepONets
can take a multi-fidelity or multi-modal input [4, 5, 6, 7, 8] in the branch
network and can use an independent network as the trunk, a network that
represents the output space, e.g. in space-time coordinates or in parametric
space in a continuous fashion. In some sense, DeepONets can be used as
surrogates in a similar fashion as reduced order models (ROMs) [9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19]. However, unlike ROMs, they are over-parametrized
which leads to both generalizability and robustness to noise that is not
possible with ROMs, see the recent work of [20].

Figure 1: Schematic of a DeepONet. $f$ represents a general input function
provided to the DeepONet. A discrete representation of $f$ with resolution
(=res) is provided as the input to the branch network. In this figure,
$y(\in\mathbb{R}^{d})$ represents the domain coordinates of the output
function.

In the present work, we investigate the possibility of using DeepONets for
representing functions over different solutions domains, something that has
not been explored before. In particular, we focus on aerodynamic design by
considering a classical problem of optimizing the shape of an airfoil at
subsonic flow conditions. Aerodynamic shape optimization (ASO) of airfoils
plays a vital role in the design of efficient modern commercial aircraft. The
aerodynamic geometric optimization process is usually applied to the airfoil
shape that forms the cross sections of the three-dimensional (3D) airfoil. The
geometric optimization is performed to reduce the drag force while increasing
the lift to enhance the aircraft’s fuel efficiency, hence reducing the
transportation cost. The optimization process often requires a numerical model
that predicts the flow field around a given airfoil geometry and computes the
lift and drag for desired flow conditions. Traditionally, the numerical models
are compressible flow numerical solvers, which are computationally intensive
to accurately realize the flow field around a complex airfoil. Surrogate
models can be introduced to circumvent the time-consuming part of the
optimization loop where the numerical solver calculates aerodynamic forces.

ASO traditionally uses two main paradigms, gradient-based and gradient-free
optimization approaches. The gradient-based approaches require the calculation
of the cost function derivative with respect to the design variables. When the
number of design variables exceeds a certain threshold, the gradient-based
optimization becomes infeasible due to its expensive computational cost [21].
The adjoint formulation is derived from either the Euler or Navier-Stokes
equations to make the gradient optimization independent of the design
variables [22, 23, 24, 25]. Solving the adjoint equations can be as time-
consuming as solving the governing equations (i.e., the rule-of-thumb is that
forward and adjoint solutions together is at least twice the cost of the
forward solver). Also, the optimization process can fall into local minima
leading to a non-optimized geometry [26]. Gradient-free approaches [27, 28]
can avoid local minima by employing the direct optimization approach that uses
many costly numerical simulations of the compressible flows. The numerous
simulations of the flow can help achieve the global minimum of the cost
function but at the expense of high computational costs. Surrogate models can
be deployed to realize the flow field with acceptable accuracy and an immense
speedup compared to the full CFD simulations. The surrogate model can then be
used in both a gradient-based of a gradient-free global optimization process.
The surrogate-based models are usually coupled with gradient-free
optimizations such as genetic algorithm and particle swarms optimization (PSO)
[29] methods. Krige [30] proposed the Kriging surrogate model that is employed
in aerospace design experiments [31, 27]. The Kriging surrogate model must be
trained both before and during the optimization process since the surrogate
model is usually inaccurate based solely on the initial training. High
generalization error of the Kriging surrogate model results in inaccurate
model prediction by the initial training stage. A deep operator network
(DeepONet), as proposed in this work, can alleviate this issue since it maps a
function to another function, which improves the generalization error
significantly [1].

The parameterization of the geometry drastically affects ASO’s computational
cost and accuracy. Parameterizing the geometry reduces the number of design
variables through which the optimization algorithm must search. Reducing the
number of design variables simplifies the optimization process and decreases
the sensitivity to noise. The parametric model must reproduce a wide range of
airfoil shapes and keep the number of design variables minimal. Various
geometric parametric models have been employed for ASO. Carpentieri _et al._
[23] employed orthogonal Chebyshev polynomials to construct the airfoil
curves. An orthogonal polynomial is used to cover the entire design space.
Lepine _et al._ [32] used Non-Uniform Rational Basis Spline (NURBS) to
parameterize a large class of airfoil shapes by only using $13$ control
points. Following the Lepine _et al._ [32] idea, Srinath and Mittal [25] also
employed NURBS for the parameterization. Other researchers have used B-Spline
[33], and Bezier [34] curves, Hicks and Henn’s functions [35] for airfoil
shape construction. Painchaud-Ouellet _et al._ [36] used NURBS for the shape
optimization of an airfoil within transonic regimes. They showed that using
NURBS ensures the regularity of the airfoil shape. The airfoil shape can also
be constructed using a deformation method [35]. This method adds a linear
combination of bumps to a baseline airfoil shape for parameterization [37,
38]. The Class function/shape function Transformation (CST) approach [39]
employs Bernstein polynomials [40] to parameterize airfoils and other
aerodynamic geometries. Other researches employed proper orthogonal
decomposition (POD) [39] to reduce the number of design variables. In the
current study, we employ the NACA 4-digits airfoil and NURBS parameterizations
for the airfoil shape construction.

With the significant advancement in computational power, Deep Neural Network
(DNN) tools have gained much attention for developing surrogate models [41,
42, 43]. The DNN approach can be readily trained for numerous input design
variables to predict the cost function of the optimization loop. Du _et al._
[44] trained a feed-forward DNN to receive airfoil shapes and predict drag and
lift coefficients. They also used RNN models for estimating the pressure
coefficient. The optimal airfoil design determined using the surrogate model
was compared with an airfoil design obtained with a CFD-based optimization
process [44]. Liao _et al._ [45] designed a surrogate model using a multi-
fidelity Convolutional Neural Network (CNN) with transfer learning. This
learning method transfers the information learned in a specific domain to a
similar field. The low-fidelity samples are taken as the source, and the high-
fidelity ones are assigned as targets. Tao and Sun [46] introduced a Deep
Belief Network (DBN) to be trained with low-fidelity data. The trained DBN was
later combined with high-fidelity data using regression to create a surrogate
model for shape optimization. Existing surrogate models for shape optimization
are all trained to predict lift, drag, or pressure coefficients. In contrast,
the flow field around the aerodynamic shape is not inferred. Here, we
construct a surrogate model that predicts the flow field around the airfoil
shape using a DeepONet. Predicting the flow field provides additional
information that can be used in the cost function of the optimization loop. We
aim to develop an aerodynamic shape optimization framework using a surrogate
model that can infer the flow field around the geometry. The surrogate model
is constructed using a DeepONet and is trained using high-fidelity CFD
simulations of airfoils in a subsonic flow regime. The surrogate model is then
implemented in two different optimization frameworks for shape optimization.
The novelties of this study include the following:

  * •

Generating a DeepONet-based surrogate model, which is an efficient and
inexpensive instantiation of the exorbitant CFD solver.

  * •

The construction of the surrogate model is invariant to input space, which can
be defined as low or high-dimensional parameterizations.

  * •

Prediction of high-dimensional flow field can be used for various cost
functions in the optimization loop.

  * •

Drag and lift coefficients are computed using the inferred high-dimensional
flow field, resulting in more accurate predictions.

  * •

Integration of the Dakota optimization framework with the DeepOnet surrogate.

The rest of the article is organized as follows. We begin by defining the
optimization process. We then present the data generation for training the
surrogate model using the open-source spectral/$hp$ element Nektar++ CFD
solver. The training procedure of the DeepONet-based surrogate model is
explained. Later, the optimization results using two different methods are
represented. Finally, we summarize our findings in the Conclusions section. In
the Appendix, we provide verification on the accuracy of the data generated by
repeating selected simulations using different codes. Additionally, we provide
validation of the Dakota optimizer by comparing against multiple approachs,
and find that all the approaches studied converge to the same solution.

##  2 Problem setup

To highlight the capabilities of DeepONets as function-to-function maps that
can be used within the airfoil shape optimization process, we start by
reviewing the traditional end-to-end shape optimization pipeline augmented
with DeepONet training. A schematic of the pipeline is shown in Figure 2.
Reviewing the figure from upper left to lower right, we start with an
experimental setup. This represents the determination of the feasible set from
which the parametric airfoils in training will be drawn, the aerodynamic
conditions, and any other engineering constraints related to the problem. We
choose NACA four-digit airfoils as our geometric representation, which
provides the upper and lower surface equations, given a random draw of
parameters. This representation is then used in three places: (1) to directly
mesh the flowfield around the airfoil for which we use Gmsh; (2) to use in
querying the surrogate model on the surface of the airfoil when predicting the
objective lift-to-drag ratio; and (3) as the branch input to the DeepONet
function-to-function map, which is pre-processed using NURBS to lower the
dimensionality.

In terms of creating the surrogate model, this can be viewed as the following
forward problem. One deciphers the geometric and flow parameters, which are
then used to create the inputs to a CFD solver: a geometric representation of
the airfoil and its corresponding mesh to be used for approximating the
flowfield and a flow parameter file. These are then input to a flow solver –
in our case, the CFD solver Nektar++. Results from this solver are used to
generate training data used to train our DeepONet surrogate. Finally, the
lower-right quadrant of the diagram denotes the shape optimization process
using the trained DeepONet surrogate. This process is recursive as the
optimizer, for which we use Dakota, queries the DeepOnet surrogate model with
new design parameters until the objective function is sufficiently minimized.
This process differs from existing approaches because the bulk of the
optimization process is done offline. Generating data using the CFD solver can
be expensive, but with the final trained DeepOnet, the online cost of
geometric optimization is orders of magnitude faster than other methods.
Furthermore, as long as the objective can be created by the DeepONet
flowfields trained over, a different objective function can be defined, not
only lift-to-drag, and an airfoil can be quickly optimized with respect to the
new objective without any additional cost.

For the purposes of our experiments, we have focused on using 2D compressible
Navier-Stokes fields at Reynolds number $Re=500$ and Mach number $Ma=0.5$ for
our training. These values have been chosen to allow us to focus on DeepONet’s
ability to capture variations in domains (instead of the compounding effects
of unsteadiness, etc.). Given the success of DeepONets under this experimental
setup, future work will extend this pipeline to more complex flows and
experimental conditions, such as varying the angle of attack, morphing
geometry, handling unsteady flow, or going into the high speed flow regimes.

Figure 2: Diagram of airfoil geometry optimization with a DeepONet surrogate
model. The blue ovals indicate the beginning and end states of the method. The
main intermediate steps are highlighted in colored boxes, whereas the
auxiliary steps are grey parallelograms.

##  3 Methodology

###  3.1 Data Generation

For each example in the dataset, we define a set of geometric parameters
$(\xi_{g})$ and flow parameters $(\xi_{f})$. The geometric parameters are then
converted into another representation, such as surface coordinates derived
from the NACA airfoil equations; this transformation is given by
$\Gamma(\xi_{g})$. These points are then used to mesh the flowfield domain
with Gmsh [47], which is then input into the flowfield simulation software
Nektar++ [48, 49]. We obtain the solutions fields from Nektar++ through post-
processing for density, x-velocity, y-velocity, and pressure. This is saved in
two sets, one in a subdomain around the airfoil for training the DeepONet, and
one at airfoil sensors on the surface for validation. We also save the
Nektar++ lift and drag forces for validation of the discrete integration of
the airfoil forces.

####  3.1.1 Geometry generation

NACA 4-digit airfoils provide an excellent testbed application for geometry
optimization using DeepONets since they can represent a wide range of shapes
from well-known and studied parameterized geometric equations. Our geometry
optimization framework could be easily extended to any parameterized geometry,
such as NACA 5-digits or beyond. Following [50], we define the parametric
equations for the surface of an airfoil with a chord length of one, as
follows:

$\displaystyle
y_{t}=\frac{t}{0.2}\left(a_{0}\sqrt{x}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+a_{x}x^{4}\right)$
(1) $\displaystyle
y_{c}=\begin{cases}\frac{m}{p^{2}}\left(2px-x^{2}\right)\quad&\text{if}\;x<p\\\
\frac{m}{\left(1-p\right)^{2}}\left(1-2p+2px-x^{2}\right)\quad&\text{if}\;x>p\\\
\end{cases}$ $\displaystyle\theta=\tan^{-1}\left(\frac{dy_{c}}{dx}\right)$ (2)
$\displaystyle x_{u}=x-y_{t}\sin(\theta),\quad y_{u}=y_{c}+y_{t}\cos(\theta)$
$\displaystyle x_{l}=x+y_{t}\sin(\theta),\quad y_{l}=y_{c}-y_{t}\cos(\theta)$
(3)

where $a_{0}=0.2969,a_{1}=-0.1260,a_{2}=-0.3516,a_{3}=0.2843,a_{4}=-0.1015$.
We can therefore define our geometry as a point cloud with coordinate sets
$\left(x_{u},y_{u},x_{l},y_{l}\right)$ parameterized by
$\xi_{g}=\left(t,p,m\right)$. A series of $x$ locations are found using cosine
spacing with $100$ points; this increases the geometric fidelity around the
leading and trailing edge, increasing the accuracy of the mesh and flowfield
simulation at these important locations. To simplify the problem and reduce
the likelihood of flow separation or turbulence, we constrain $\xi_{g}$. The
maximum thickness ($t$) is set to a constant $0.15$, and the domain of the
parametric space left by the position of maximum camber ($p$) and maximum
camber ($m$) is $p\times m\in[0.2,0.5]\times[0.0,0.09]$. Therefore, for one
geometric example in either the train or test set, we draw a $\xi_{g}$ tuple
where the variable parameters are drawn from a uniform distribution within
their domains. We perform this draw $50$ times and obtain the surface
coordinates from Equations 1 \- 3, splitting it into $40$ training and $10$
testing examples as seen in Fig. 3. The test/train split is an essential
aspect of deep learning; here, we choose a relatively sparse sampling
highlighting the ability of DeepONets to generalize well to unseen parameters.

Next, to lower the input dimensionality into the DeepONet branch, we fit the
airfoil surface with Non-Uniform Rational B-Splines (NURBS) with $30$ control
points using geomdl [51]. This reduces the input dimensionality from $200$
$(x,y)$ pairs to only $30$.

Figure 3: Train and test set geometries sampled over the $\xi_{g}$ domain.
Note that the numbers correspond to NACA airfoils and duplicate numbers are
due to rounding to the nearest integer for readability. The true NACA
parameters are drawn from a uniform distribution and are, therefore, real-
valued.

####  3.1.2 Mesh generation

The meshes are generated using Gmsh [47] for parametrized airfoils geometry
with a minimum characteristic length of $0.01$ at parametrized locations. The
airfoil flowfields are meshed using the $200$ surface points exactly from the
NACA equations, not the NURBS fit, which contains some inaccuracy. A spline
function is used to represent the 1D geometry of the boundaries of airfoils.
To resolve the flow at leading and trailing edges, meshes are refined by using
the splitting approach [52]. The NURBS low-dimensional representation is a
step to reduce overparameterization in the DeepONet, which is unnecessary for
generating high-fidelity training data. The mesh generation for all $50$
airfoils is automated using Gmsh’s Python API integrated with the geometry
generation in Python, so no manual operations are needed. Fig. 4 shows the
mesh of the entire simulated domain $\Omega_{S}$, which is then input into
Nektar++ along with the flow parameters to generate the DeepONet training
data. As seen in the figure, only a subset of the solved steady-state domain
$\Omega_{T}$ is used in training the DeepONet. This is because the DeepOnet is
a function-to-function map and does not strictly obey boundary conditions or
is affected by phenomena such as reflections due to the boundaries. It is
performing regression on the dataset, not solving the system of equations, and
therefore can be taken as a smaller domain, even without freestream
conditions. Since the objective is geometric optimization, this subdomain
simplifies the DeepONet training problem and cost of training.

Figure 4: Discretization of an airfoil with bounding domain. $\Omega_{\rm S}$
represents the entire domain over which the compressible Navier-Stokes
equations are solved. $\Omega_{\rm T}$ (inset image with red border)
represents the domain for which a DeepONet is trained. Inset image with yellow
boundary represents the subdomain $(\Omega_{\rm W})$ showing the mesh
refinement around the airfoil.

####  3.1.3 Flowfield simulation

We used a compressible flow solver implemented in Nektar++ to generate the
flow field data. The Nektar++ (www.nektar.info) is an open-source software
framework [48, 49] designed to support the development of high-performance,
scalable solvers for partial differential equations using the spectral/hp
element method. The 2D compressible flow solver uses the two-dimensional
compressible Navier-Stokes equations expressed as,

$\displaystyle\begin{aligned}
&\frac{\partial{\rho}}{\partial{t}}+\frac{\partial{\rho}{u}}{\partial{x}}+\frac{\partial{\rho}{v}}{\partial{y}}=0\\\
&\frac{\partial{\rho}{u}}{\partial{t}}+\frac{\partial{\rho}{u}^{2}+{p}}{\partial{x}}+\frac{\partial{\rho}{u}{v}}{\partial{y}}=\frac{1}{Re}\left(\frac{\partial{\tau}_{xx}}{\partial{x}}+\frac{\partial{\tau}_{yx}}{\partial{y}}\right)\\\
&\frac{\partial{\rho}{v}}{\partial{t}}+\frac{\partial{\rho}{u}{v}}{\partial{x}}+\frac{\partial{\rho}{v}^{2}+{p}}{\partial{y}}=\frac{1}{Re}\left(\frac{\partial{\tau}_{xy}}{\partial{x}}+\frac{\partial{\tau}_{yy}}{\partial{y}}\right)\\\
&\frac{\partial{E}}{\partial{t}}+\frac{\partial({E}+{p}){u}}{\partial{x}}+\frac{\partial({E}+{p}){v}}{\partial{y}}=\frac{1}{Re}\left[\frac{\partial\left({u}{{\tau}}_{xx}+{v}{{\tau}}_{xy}+{{\kappa}}\tfrac{\partial{T}}{\partial{x}}\right)}{\partial{x}}+\frac{\partial\left({u}{{\tau}}_{xy}+{v}{{\tau}}_{yy}+{{\kappa}}\tfrac{\partial{T}}{\partial{y}}\right)}{\partial{y}}\right],\end{aligned}$
(4)

where ${p}={R}{\rho}{T}$, ${R}=\frac{1}{\gamma Ma^{2}}$,
${k}=\frac{\gamma}{\gamma-1}\frac{{{\mu}}{R}}{Pr}$, with ${{\mu}}$ being the
non-dimensional viscosity and computed by using the Sutherland law as

$\displaystyle{{\mu}}=\frac{{T}^{3/2}}{Re}\frac{1+C/T_{\infty}}{{T}+C/T_{\infty}}.$
(5)

Expressions for $\tau_{xx}$, $\tau_{yy}$, and $\tau_{xy}$ are as follows

$\displaystyle\tau_{xx}$
$\displaystyle=2\mu\left(u_{x}-\frac{u_{x}+v_{y}}{3}\right),$
$\displaystyle\tau_{yy}$
$\displaystyle=2\mu\left(v_{y}-\frac{{u}_{x}+{v}_{y}}{3}\right),$
$\displaystyle\tau_{xy}$ $\displaystyle=\mu\left(u_{y}+{v}_{x}\right).$

We aim to simulate the flow past airfoils by solving the compressible Navier-
Stokes equations given by Equation (4) with free-stream parameters
$M_{\infty}=0.5$,
$Re_{L=1}=500,~{}u_{\infty}=1,~{}v_{\infty}=0,~{}T_{\infty}=1$, ${\rm AoA=0}$
and $Pr=0.72$. The flow domain is $[-3,11]\times[-3,3]$ and discretized by
conforming triangular elements. To solve (4), we use the discontinuous
Galerkin spectral element method (DGSEM) with basis functions spanned in 2D by
Legendre polynomials of the second degree. For advection and diffusion terms,
weak and interior penalty-based dG approach with Roe upwinding is used in
space. A diagonally implicit Runge–Kutta (DIRK) method is used as a time
integrator for advection and diffusion terms. The boundary conditions of
inflow, outflow, adiabatic wall at airfoil surface, and high-order boundary
conditions at the top and bottom are imposed weekly. For detailed descriptions
of solvers and methods, readers are advised to read the article by [53].

###  3.2 Surrogate model: DeepONet

####  3.2.1 Brief review of DeepONets

Neural operators are neural network models developed on the basis of the
universal operator approximation theorem [54]. The neural operators learn the
mapping between spaces of function and directly learn the underlying operator
from the available training data. DeepONets [1] and Fourier Neural Operators
(FNO) [55] are the two popular neural operators extensively used for solving a
wide spectrum of problems in diverse scientific areas. A DeepONet consists of
a branch network that encodes the input function and a trunk network that
learns a collection of basis functions. The DeepONet output is computed by
taking the inner product between the branch and trunk network outputs.

####  3.2.2 Training and Testing of DeepONets

We train four different DeepONet models to learn the pressure ($p$), density
($\rho$), and velocity ($u$, $v$) fields for a given airfoil geometry
($\xi_{g}$) from the training data. The trunk network learns a collection of
basis ($\bm{\phi}$) as functions of spatial coordinates, and the branch
network learns the corresponding coefficients ($\bm{\alpha}$) as a function of
the airfoil geometry. The DeepONet output is defined as

$\mathcal{G}^{q}(\xi_{g})(x,y)=\sum_{i=1}^{N_{\phi}}\alpha_{i}(\xi_{g};\theta_{b}^{q})\phi_{i}(x,y;\theta_{t}^{q})\quad
q\in\\{p,\rho,u,v\\}.$ (6)

In the case of airfoils, the geometry can be fed into the branch network by
directly providing the geometric parameter, $\xi_{g}$. However, $\xi_{g}$ need
not exist explicitly for a general arbitrary geometry. Under such a scenario,
the geometry is often represented using the NURBS control points. To
demonstrate the effectiveness of using a DeepONet in either of the situations,
we investigate parameter-DeepONet that directly takes $\xi_{g}$ as the branch
network input and NURBS-DeepONet that takes NURBS control points of the
airfoil geometry as the input to the branch network. The hyperparameters of
the DeepONet used in this study are provided in Table 2.

Table 2: Hyperparameters of NURBS-DeepONet and the Parameter DeepONet NURBS-DeepONet |  | Parameter-DeepONet
---|---|---
Branch network architecture: | [30,100,100,50] |  | Branch network architecture: | [2,100,100,50]
Branch network activation: | tanh |  | Branch network activation: | tanh
Trunk network architecture: | [2,100,100,50] |  | Trunk network architecture: | [2,100,100,50]
Trunk network activation: | tanh |  | Trunk network activation: | tanh
$N_{\phi}$: | 50 |  |  $N_{\phi}$: | 50
Optimizer: | Adam |  | Optimizer: | Adam
Learning rate: | 1.00E-04 |  | Learning rate: | 1.00E-04

####  3.2.3 Lift and Drag calculation

We present two ways to evaluate the objective function – in our case, lift-to-
drag. The map created by the trained DeepONets takes an input geometry and
spatial $(x,y)$ location in the output space and returns the value of the
field at that location. So it follows that to estimate the lift and drag, we
simply evaluate the discrete integral over the surface of the airfoil, for
which the $(x,y)$ points are easily generated for objective function queries
by Equations 1-3. The discrete integrals for lift and drag given by Equations
8 & 9 subject to the approximation of wall shear stress in Equation 7.

$\displaystyle\tau_{w}=\mu\frac{dU}{d\overrightarrow{n}}$ (7) $\displaystyle
L=\int dF_{x}=\sum
p\overrightarrow{n_{x}}dS+\sum\tau_{w}\overrightarrow{t_{x}}dS$ (8)
$\displaystyle D=\int dF_{y}=\sum
p\overrightarrow{n_{y}}dS+\sum\tau_{w}\overrightarrow{t_{y}}dS$ (9)

For the first term in the aerodynamics forces, the pressure is directly
obtained from one of the DeepONet predictions. In the second term, the
viscosity $\mu(p,\rho)$ is obtained by Equation 5 as a function of the
pressure and density DeepONets. Finally, the change in speed of the flow over
the airfoil surface $\frac{dU}{d\overrightarrow{n}}$, is obtained in two ways:

  1. 1.

Finite-difference (A):

$\displaystyle\frac{dU}{d\overrightarrow{n}}=\frac{-U_{2}+4U_{1}-3U_{0}}{2h}$
(10)

  2. 2.

Automatic-differentiation (B):

$\displaystyle\frac{dU}{d\overrightarrow{n}}=(v_{y}-u_{x})sin\theta cos\theta-
v_{x}sin^{2}\theta+u_{y}cos^{2}\theta$ (11)

where $h=0.001$ and $\theta$ is the angle between the x-y axis and each
segment’s normal-tangental axis. Approach (A) is a second-order forward finite
difference approximation obtained by sampling the x-velocity and y-velocity
DeepONets at the appropriate locations defined by the surface normal and
spacing $h$. Approach (B) utilizes the now well-known development in automatic
differentiation [56], primarily utilized in physics-informed machine learning
for approximating partial derivatives to obtain the PDE residual. Here, since
the DeepONet directly takes in the spatial $(x,y)$ coordinates and outputs the
velocity components $(u,v)$, the computational graph is complete and the
partials $(u_{x},u_{y},v_{x},v_{y})$ can be estimated with this method. The
required sampling for each method is shown in Fig. 5. While (A) takes three
times the amount of point evaluations, (B) requires the gradients to be
computed, so the cost of each can be viewed as similar. However, the
flexibility of using automatic differentiation in this way may allow for more
complex objective functions in the future given the right mapping and
subsequent computational graph.

Figure 5: Illustration of the discrete integral for lift and drag. The points
in red indicate the surrogate model samples used in the construction of the
approximation to $\frac{dU}{d\overrightarrow{n}}$ with a finite difference
approximation of the gradient or the automatic differentiation approximation
using the direct network gradients.

##  4 Results

###  4.1 Results from DeepONet model

The training and testing relative $L^{2}$ error of the NURBS and parameter-
based DeepONets for all four different fields are reported in table 3. We
observe that the NURBS and parameter-based DeepONets predict fields with
similar accuracy. NURBS-DeepONet has marginally better predictions on the
pressure field, while parameter-DeepONet generates marginally better
predictions for velocity and density fields. The main takeaway is that either
representation is sufficient for the geometry optimization of this experiment.
However, we must consider that in the future more complex geometries may be
used, particularly in the sense of local morphing. Therefore, the NURBS
representation will likely be necessary as a direct parameter mapping may miss
local nuances. The predicted fields and the absolute pointwise error by the
best DeepONet models for the flowfield parameters are shown in Fig. 7. It can
be seen that the global prediction is, in general, accurate; the error is
primarily localized to the airfoil’s leading edge. In the future, adaptive
weighting schemes will be used to improve DeepONet training, particularly at
the points of difficulty, such as the surface and leading edge. The
corresponding relative $L^{2}$ error of the fields over the entire dataset is
shown in Fig. 6. We observe minimal generalization error and that DeepONets
are globally accurate.

Table 3: Relative $L^{2}$ errors of the state variables trained DeepONet models | NURBS-DeepONet | Parameter-DeepONet
---|---|---
| Train rel. $L^{2}$ Error | Test rel. $L^{2}$ Error | Train rel. $L^{2}$ Error | Test rel. $L^{2}$ Error
$\mathcal{G}^{p}$ | 4.68e-03 | 6.05e-03 | 5.23e-03 | 6.85e-03
$\mathcal{G}^{u}$ | 4.97e-03 | 6.21e-03 | 4.12e-03 | 5.38e-03
$\mathcal{G}^{v}$ | 3.73e-03 | 4.60e-03 | 3.31e-03 | 4.25e-03
$\mathcal{G}^{\rho}$ | 4.57e-03 | 5.89e-03 | 4.00e-03 | 5.18e-03
Figure 6: Error Scatter Plots. The relative $L^{2}$ error corresponding to
each of the train and test samples for pressure, velocity, and density fields,
with respect to the DeepONet predictions, are shown in this figure. Figure 7:
DeepONet Predictions. The pressure, density, and velocity fields around the
test set airfoil NACA 7315 predicted by the DeepONet, and the corresponding
pointwise absolute errors are also provided.

Regarding geometry optimization, we must concern ourselves not only with the
global flowfield accuracy, but also with the accuracy on the surface of the
airfoil in particular. Fig. 8 shows the corresponding surface prediction plots
for the same airfoil presented in Fig. 7 as a function of the x-direction over
the airfoil. We observe good accuracy in the pressure and density fields,
which will provide very accurate predictions of the lift and drag force
components due to pressure as well as the viscosity $\mu(p,\rho)$, which is a
function of these fields per Equation 5. The surface’s $x$ and $y$ velocity
fields do not agree because the DeepONet does not strictly obey a no-slip
condition. However, aside from the leading edge, the predictions are close to
zero. Furthermore, the fields are not directly related to the objective lift-
to-drag but indirectly related through the estimate of the change in flow
speed over the surface $\frac{dU}{d\overrightarrow{n}}$ obtained by approaches
(A) and (B) in Equation 10 and 11. Therefore, the inaccuracy does not
significantly affect the overall objective prediction. This is corroborated by
Fig. 9, which shows the sorted lift-to-drag ratio for the entire dataset. The
results of both numerical integration approaches (A) and (B) are very accurate
to the stored lift-to-drag results from the Nektar++ data generation step. We
can also see in the error plot that it is quite uniform, and there are no
discernible biases in the geometric parameter space, indicating we have
learned the entire space well enough for the final optimized result to be
accurate.

Figure 8: Plot of the flowfield variables on the surface of the test set
airfoil NACA $\mathbf{7315}$. All plots display accurate predictions on the
surface, which are then used to compute the lift and drag forces. The no-slip
condition is not directly enforced by the DeepONet, which results in the
velocity plot difference. However, it can be seen by the y-scale that the
prediction is close to zero, aside from the leading edge, and does not greatly
affect the overall lift and drag computation. Figure 9: Plot of the computed
lift-to-drag objective for the entire dataset sorted by the Nektar++ reference
values. As seen in both plots, the approximation to the high-fidelity CFD
solution is very accurate and consistent throughout the entire parametric
domain. Particularly, we note that the testing set performs comparably to the
training set, meaning there is little to no generalization error, which is
necessary when inferring unseen queried geometries during optimization.

###  4.2 Shape Optimization Results

The objective of shape optimization is to maximize the lift-to-drag ratio over
a feasible region of parameters, which are $m$ and $p$ for this case. Equation
4.2 gives this objective in the form of a minimization problem, as is standard
for most optimizers that perform gradient-based or gradient-free optimization.
Therefore, the definition of a constrained optimization problem for airfoil is
expressed as m,p-f(m, p) m_min ≤m ≤m_max p_min ≤p ≤p_max,

where $f(m,p)$ represents ratio of lift to drag and $[m_{min},m_{max}]$ and
$[p_{min},m_{max}]$ are bounds for feasible search region.

One of the present study’s goals is to optimize the shape for any arbitrary
geometry. Therefore, we integrated the DeepONet-based surrogate model with
Dakota, which is a multilevel parallel object-oriented framework for design
optimization, parameter estimation, uncertainty quantification, and
sensitivity analysis [57]. Dakota is freely available and offers a very
efficient and scalable implementation. We integrated the DeepONet with Dakota
in a modular approach as shown in Algorithm 1, where $\mathcal{D}$ is an
algorithm chosen from a set of optimizers provided by Dakota and DeepONet-
based model $\Phi$ is passed as an argument to $\mathcal{D}$. For example, to
achieve the solution of Equation (4.2), we use an efficient global algorithm
(EGO), which is a derivative-free approach that uses a Gaussian process model
for the optimization of the expected improvement function and is based on the
NCSU Direct algorithm [58]. The reason behind choosing this method is to avoid
tuning various hyperparameters. To use the algorithm to solve the problem in
Equation (4.2), we set a seed, which is to be used for Latin Hypercube
Sampling (LHS) to generate the initial set of points for constructing the
initial Gaussian process. To gain efficiency, we used batch-sequential
parallelization offered by Dakota on an eight-core CPU (2.3 GHz Intel core
i9).

Algorithm 1 Integration of DeepONet-based surrogate model with Dakota

$m$, $p$, $x$, $y$: maximum camber, the position of maximum camber, spatial
coordinates for flow-field prediction

Trained DeepONet: $\mathcal{G}(m,p,x,y)$

$\mathcal{D}$ : $\mathcal{D}\in{\text{Algorithms in Dakota}}$

$\Psi(u,v,\rho,p,\xi_{g})$: Function producing the lift $L$ and drag $D$

$N$: Number of objective function evaluations

$m_{\min},m_{\max},p_{min},p_{\max}$: Bounds of feasible region

$n\leftarrow 0$, $m_{\text{opt}}\leftarrow m_{0}$, $p_{\text{opt}}\leftarrow
p_{0}$ $\triangleright$ Initialize

while $n~{}!=~{}N$ do

$m_{\text{opt}},p_{\text{opt}}\leftarrow\mathcal{D}(-f(m_{\text{opt}},p_{\text{opt}}),m_{\min},m_{\max},p_{min},p_{\max},\mathcal{G},\Psi)$
$\triangleright$ Optimization process for parameters

$n\leftarrow n+1$

end while

The optimization landscapes for approaches (A) and (B) are shown in Fig. 10
obtained by brute force evaluation of the respective objectives. Also plotted
are the locations of the dataset in the parameter space
$(p,m)\in[0.2,0.5]\times[0.0,0.09]$, displaying the sparse sampling used to
obtain accurate optimization results. We also observe that the landscape for
this experimental setup is simple and convex. In future work, more complex
conditions like local morphing will make the landscape more complex, likely
nonconvex, requiring sophisticated optimization methods. While not needed
here, we still utilize state-of-the-art optimization in our framework with
Dakota. Additionally, we can see that the local minimum given the constrained
parameter bounds, set to the dataset sampling bounds, is at the edge. This
implies that the global minimum lies outside of our trained bounds, which is
not known a priori when generating data. In future work, we hope to evaluate
DeepONets efficacy in extrapolating outside the trained bounds, potentially
using transfer learning.

The optimization process finds the minimizer of function in $15$ evaluations,
and the optimum value which maximizes $L/D$ is $(m_{*},p_{*})=(0.2,0.067)$. To
achieve consistency in the optimization process, we ran the EGO algorithm 15
times with different seeds and we observed the same optimal point.
Furthermore, we validated the optimization results with a comparison to other
approaches in Appendix 5.2. The results obtained from all the approaches are
in excellent agreement and reported in detail in Table 6 along with their wall
clock times. Finally, the most significant contribution of the paper is shown
in Table 4. As we can see, the integration of DeepONet into an airfoil
geometry optimization framework has lowered the online cost of new objective
evaluations by $32,000+$ times. This makes it entirely possible to have almost
real-time optimization results, costing a few minutes instead of days.
Furthermore, the trained models can be put on any hardware, such as a standard
laptop, and real-time accurate flowfields can be predicted in seconds, meaning
the geometry optimization is not hardware dependent at test time. We have
demonstrated that integrating DeepONets into a geometry optimization pipeline
suffers little in accuracy and provides the tradeoff of obtaining and training
on an offline dataset with almost instantaneous optimization results when used
online compared to a traditional CFD method.

Figure 10: Visualization of the lift/drag landscape obtained from brute force sampling of $p$ and $m$ using a $10\times 10$ grid. The train and test sets are also plotted to show the sparse dataset used by the DeepONet (A) Landscape obtained using finite difference approximation. (B) Landscape obtained using automatic differentiation approximation. Model Type | Relative Cost of Single Objective Function Evaluation
---|---
Baseline CFD (Nektar++ placeholder) | $32,253$
DeepONet (A) | $1.34$
DeepONet (B) | $1$
Table 4: Relative cost of single objective function evaluation during geometry
optimization A Flowfield mapping with finite-difference approximation. B
Flowfield mapping with automatic-differentiation approximation.

To validate the parameters of optimized airfoil $(p=0.2,0.067)$, we compare
the streamline plots in Fig. 11 constructed using the flowfields $(u,v)$
obtained from trained DeepONet and Nektar++. In general, the streamline plots
show an excellent agreement and therefore validate the workflow of shape
optimization presented in this work. DeepONet is able to successfully detect
and predict the circulation near the trailing edge of the airfoil. However, a
closer look reveals non-physical streamlines originating from the surface of
the airfoil. This is due to the error in the velocity fields near the airfoil
surface, predicted by the DeepONet. The results can be further improved by
giving more importance to the region near the airfoil surface (via proper
weighting) during training of the DeepONet.

Figure 11: Optimized airfoil.  The streamline plot represents the flow past
the optimized airfoil with $(p,m)=(0.2,0.067)$ in $\Omega_{W}$. The velocity
fields simulated in Nektar++ and predicted by the DeepONets are shown here.

##  5 Conclusions

We have successfully integrated DeepONets as a surrogate model into the shape
optimization framework for airfoils. Having summarized prior work in this
field, we empirically demonstrate the efficacy of DeepONets in terms of
retaining sufficient flowfield accuracy used in evaluating the objective
function of lift-to-drag, as well as the significant computational speed up as
a replacement for a traditional CFD solver during online geometry
optimization. We have provided thorough validation of the results presented as
well as extensive experimentation such as two approaches (A) and (B) when
approximating the wall shear stress and two forms of DeepONet inputs (NURBS
and $\xi_{g}$) to ensure a robust pipeline. Importantly, DeepONets exhibit
almost no generalization error over the dataset, so it follows that the
resulting optimized geometry ($p=0.2,m=0.067$) is accurate and achieved
$32,253$ speed-up compared to the CFD baseline. The framework is general and
can address more complex problems with multiple inputs, e.g. different Mach
numbers and different angles of attack that can be inputs to either the branch
or the trunk networks. Hence, with relatively small modifications, such a
framework can handle optimization in the high speed flow regimes that exhibit
flow unsteadiness, shocks, non-equilibrium chemistry, and even morphing
geometry. Furthermore, the approaches presented are flexible due to the
integration of machine learning in the form of function-to-function maps using
DeepONet. Therefore, improvements such as the introduction of multi-fidelity
training and physics-informed machine learning can be leveraged to reduce the
cost of data generation. We also successfully show the application of
automatic differentiation, which performs comparably to the traditional
approach of finite differences in the wall shear stress calculation. Finally,
we hope to utilize transfer learning and uncertainty quantification using the
recently developed library NeuralUQ [59] to extrapolate outside of the trained
geometric parameter domain to find global optima with confidence.

## Acknowledgments

The research reported in this document is performed in connection with
cooperative agreement contract/instrument W911NF-22-2-0047 with the U.S. Army
Research Laboratory. L.B. and A.G. were supported by the US Army Research
Laboratory 6.1 basic research program in vehicle power and propulsion
sciences. The views and conclusions contained in this document are those of
the authors and should not be interpreted as representing the official
policies or positions, either expressed or implied, of the U.S. Army Research
Laboratory or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for Government purposes notwithstanding any
copyright notation herein.

## Appendix

###  5.1 Nektar++ cross-verification

To ensure a steady-state solution, we selected the flow parameters as $M=0.5$
and $Re=500$. Calculating the steady-state solution requires careful
consideration. To ensure the steady state solution, for each NACA profile, we
recorded the value of conservative variables at six different locations in the
wake. We then examined the conservative variables’ time history to ensure the
steady state is reached and the solution update has stopped. Fig. 12 shows the
location of the history points, where the transient flow could last longer
than other spatial locations. As a sample, we plot the time history of
variables at point 5 (Fig. 13), which experiences the highest flow
fluctuations in time. According to Fig. 13, the solution reached a stationary
state when the conservative variables approached constant values.

Figure 12: History points locations in the wake for NACA airfoils. The conservative variable values are stored in time at these locations to monitor the steady-state solution. |
---|---
(a) Density | (b) x-momentum
|
(c) y-momentum | (d) Total energy
Figure 13: Time history of (a) $\rho$, (b) $\rho u$, (c) $\rho v$, and (d) $E$
at point 5 (see Fig. 12) for the NACA7315 airfoil. The profiles show that the
simulation has reached a steady-state solution where the flow variables reach
constant states.

We also validated the simulation setup of NACA airfoils in NekTar++ with the
results obtained by an in-house code based on the discontinuous spectral
element method of Kopriva [60, 61]. We selected the NACA0020 airfoil for
cross-verification. The steady-state solutions of the flow with $M=0.5$ and
$Re=500$ are computed, and the results are shown in Fig. 14. The flow field
primitive variables computed by both solvers agree and show the validity of
the NekTar++ simulation setup, including mesh and simulation parameters.
According to Fig. 14, the number of elements employed for the NekTar++
simulations is sufficient for an accurate prediction. After validating the
flow field, we performed an extra flow simulation around the NACA4402 airfoil.
The drag and lift coefficients of this airfoil are reported by Kunz [62] but
for an incompressible flow at $Re=1000$.

We employed the automatic mesh generation setup used for the training set to
create the mesh and used the same simulation setup as the training set. We
computed the drag and lift coefficients and compared them with the literature.
Table. 5 compares the drag and lift coefficients computed by Nektar++ with
values reported by Kunz [62] for a similar but not exactly the same set up.

|
---|---
(a) x-Velocity Nektar++ | (b) x-Velocity DSEM
|
(c) y-Velocity Nektar++ | (d) y-Velocity DSEM
|
(e) Density Nektar++ | (f) Density DSEM
|
(g) Pressure Nektar++ | (h) Pressure DSEM
Figure 14: Cross verification of steady state flow around NACA0020 obtained by the DSEM code and Nektar++. The flow parameters are set as $Re=500$ and $M=0.5$. | $C_{d}$ | $C_{l}$ | Error
---|---|---|---
Nektar++ | 0.1050 | 0.1700 | $1.7\%$
Kunz [62] | 0.1032 | 0.1852 | $8.2\%$
Table 5: Comparison between the drag and lift coefficients computed by
Nektar++ with Kunz [62]. For the reference solution, the setup is different,
i.e., the Reynolds number is $1000$, and the flow is assumed to be
incompressible.

###  5.2 Geometry optimization validation

To validate the main geometry optimization findings in the manuscript using
Dakota, we also evaluate the objective function using brute force and SciPy’s
[63] dual-annealing method. Brute force is evaluated using a $10\times 10$
grid on the geometric parameter space $(p,m)\in[0.2,0.5]\times[0.0,0.09]$, and
therefore requires 100 evaluations of the objective. The dual-annealing method
is set to have a maximum amount of 50 evaluations but likely could be set to
fewer. As seen in Table 6, all methods discover the same optimal set of
parameters and $\frac{L}{D}$, regardless of optimizer or partial approximation
given by (A) and (B).

Model Type | Optimized Parameters (p, m) | Optimized $\frac{L}{D}$ | Computation Cost (sec)
---|---|---|---
Brute Force (A) | (0.2, 0.070) | 0.272 | 332.63
Brute Force (B) | (0.2, 0.070) | 0.281 | 264.12
SciPy dual-annealing (A) | (0.2, 0.064) | 0.269 | 196.09
SciPy dual-annealing (B) | (0.2, 0.064) | 0.281 | 165.57
Dakota (A) | (0.2, 0.067) | 0.269 | 157.71
Dakota (B) | (0.2, 0.063) | 0.282 | 105.00
Table 6: Geometry optimization results for different DeepONet methodologies.
(A) Flowfield mapping with finite-difference approximation. (B) Flowfield
mapping with automatic-differentiation approximation. The cost of these
methods were obtained on a eight-core CPU (2.3 GHz Intel core i9) with 16 GB
2667 MHz DDR4 of MacBook Pro 2019.

## References

  * [1] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature Machine Intelligence, 3(3):218–229, 2021.
  * [2] Lu Lu, Pengzhan Jin, and George Em Karniadakis.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  arXiv preprint arXiv:1910.03193, 2019.
  * [3] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics, 378:686–707, 2019.
  * [4] Subhayan De, Malik Hassanaly, Matthew Reynolds, Ryan N. King, and Alireza Doostan.  Bi-fidelity modeling of uncertain and partially unknown systems using deeponets. Preprint at https://arxiv.org/abs/2204.00997, 2022\.
  * [5] Amanda A. Howard, Mauro Perego, George E. Karniadakis, and Panos Stinis.  Multifidelity deep operator networks. Preprint at https://arxiv.org/abs/2204.09157, 2022\.
  * [6] Lu Lu, Raphaël Pestourie, Steven G. Johnson, and Giuseppe Romano.  Multifidelity deep neural operators for efficient learning of partial differential equations with application to fast inverse design of nanoscale heat transport.  Phys. Rev. Research, 4:023210, 2022.
  * [7] Pengzhan Jin, Shuai Meng, and Lu Lu.  Mionet: Learning multiple-input operators via tensor product.  arXiv preprint arXiv:2202.06137, 2022.
  * [8] Min Zhu, Handi Zhang, Anran Jiao, George Em Karniadakis, and Lu Lu.  Reliable extrapolation of deep neural operators informed by physics or sparse observations.  arXiv preprint arXiv:2212.06347, 2022.
  * [9] Jan S Hesthaven and Stefano Ubbiali.  Non-intrusive reduced order modeling of nonlinear problems using neural networks.  Journal of Computational Physics, 363:55–78, 2018.
  * [10] Jan S Hesthaven, Gianluigi Rozza, Benjamin Stamm, et al.  Certified reduced basis methods for parametrized partial differential equations, volume 590.  Springer, 2016.
  * [11] Peter Benner, Mario Ohlberger, Anthony Patera, Gianluigi Rozza, and Karsten Urban.  Model reduction of parametrized systems.  Springer, 2017.
  * [12] Matthew O Williams, Ioannis G Kevrekidis, and Clarence W Rowley.  A data–driven approximation of the koopman operator: Extending dynamic mode decomposition.  Journal of Nonlinear Science, 25(6):1307–1346, 2015.
  * [13] Eliodoro Chiavazzo, Charles W Gear, Carmeline J Dsilva, Neta Rabin, and Ioannis G Kevrekidis.  Reduced models in chemical kinetics via nonlinear data-mining.  Processes, 2(1):112–140, 2014.
  * [14] Chad Lieberman, Karen Willcox, and Omar Ghattas.  Parameter and state model reduction for large-scale statistical inverse problems.  SIAM Journal on Scientific Computing, 32(5):2523–2542, 2010.
  * [15] Tan Bui-Thanh, Karen Willcox, and Omar Ghattas.  Model reduction for large-scale systems with high-dimensional parametric input space.  SIAM Journal on Scientific Computing, 30(6):3270–3288, 2008.
  * [16] Peter Benner, Serkan Gugercin, and Karen Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM review, 57(4):483–531, 2015.
  * [17] David Amsallem, Matthew Zahr, Youngsoo Choi, and Charbel Farhat.  Design optimization using hyper-reduced-order models.  Structural and Multidisciplinary Optimization, 51(4):919–940, 2015\.
  * [18] Kevin Carlberg and Charbel Farhat.  A compact proper orthogonal decomposition basis for optimization-oriented reduced-order models.  In 12th AIAA/ISSMO multidisciplinary analysis and optimization conference, page 5964, 2008.
  * [19] Youngsoo Choi, Gabriele Boncoraglio, Spenser Anderson, David Amsallem, and Charbel Farhat.  Gradient-based constrained optimization using a database of linear reduced-order models.  Journal of Computational Physics, 423:109787, 2020.
  * [20] Katiana Kontolati, Somdatta Goswami, Michael D Shields, and George Em Karniadakis.  On the influence of over-parameterization in manifold based surrogates and deep neural operators.  arXiv preprint arXiv:2203.05071, 2022.
  * [21] Yin Yu, Zhoujie Lyu, Zelu Xu, and Joaquim RRA Martins.  On the influence of optimization algorithm and initial design on wing aerodynamic shape optimization.  Aerospace Science and Technology, 75:183–199, 2018.
  * [22] James Reuther, Antony Jameson, James Farmer, Luigi Martinelli, and David Saunders.  Aerodynamic shape optimization of complex aircraft configurations via an adjoint formulation.  In 34th aerospace sciences meeting and exhibit, page 94, 1996.
  * [23] Giampietro Carpentieri, Barry Koren, and Michel JL van Tooren.  Adjoint-based aerodynamic shape optimization on unstructured meshes.  Journal of Computational Physics, 224(1):267–287, 2007.
  * [24] Siva Nadarajah and Antony Jameson.  Studies of the continuous and discrete adjoint approaches to viscous automatic aerodynamic shape optimization.  In 15th AIAA computational fluid dynamics conference, page 2530, 2001.
  * [25] DN Srinath and Sanjay Mittal.  An adjoint method for shape optimization in unsteady viscous flows.  Journal of Computational Physics, 229(6):1994–2008, 2010.
  * [26] Oleg Chernukhin and David W Zingg.  Multimodality and global optimization in aerodynamic design.  AIAA journal, 51(6):1342–1354, 2013.
  * [27] Jichao Li, Jinsheng Cai, and Kun Qu.  Surrogate-based aerodynamic shape optimization with the active subspace method.  Structural and Multidisciplinary Optimization, 59(2):403–419, 2019\.
  * [28] Xiaojing Wu, Zijun Zuo, and Long Ma.  Aerodynamic data-driven surrogate-assisted teaching-learning-based optimization (tlbo) framework for constrained transonic airfoil and wing shape designs.  Aerospace, 9(10):610, 2022.
  * [29] Russell Eberhart and James Kennedy.  A new optimizer using particle swarm theory.  In MHS’95. Proceedings of the sixth international symposium on micro machine and human science, pages 39–43. Ieee, 1995.
  * [30] Daniel G Krige.  A statistical approach to some basic mine valuation problems on the witwatersrand.  Journal of the Southern African Institute of Mining and Metallurgy, 52(6):119–139, 1951.
  * [31] J Liu, W-P Song, Z-H Han, and Y Zhang.  Efficient aerodynamic shape optimization of transonic wings using a parallel infilling strategy and surrogate models.  Structural and Multidisciplinary Optimization, 55(3):925–943, 2017\.
  * [32] Jerome Lepine, Francois Guibault, Jean-Yves Trepanier, and Francois Pepin.  Optimized nonuniform rational b-spline geometrical representation for aerodynamic design of wings.  AIAA journal, 39(11):2033–2041, 2001.
  * [33] Kun Wang, Shengjiao Yu, Zheng Wang, Renzhong Feng, and Tiegang Liu.  Adjoint-based airfoil optimization with adaptive isogeometric discontinuous galerkin method.  Computer Methods in Applied Mechanics and Engineering, 344:602–625, 2019.
  * [34] Dimitrios I Papadimitriou and Costas Papadimitriou.  Aerodynamic shape optimization for minimum robust drag and lift reliability constraint.  Aerospace Science and Technology, 55:24–33, 2016.
  * [35] Raymond M Hicks and Preston A Henne.  Wing design by numerical optimization.  Journal of Aircraft, 15(7):407–412, 1978.
  * [36] Simon Painchaud-Ouellet, Christophe Tribes, Jean-Yves Trépanier, and Dominique Pelletier.  Airfoil shape optimization using a nonuniform rational b-splines parametrization under thickness constraint.  AIAA journal, 44(10):2170–2178, 2006.
  * [37] Guodong Chen and Krzysztof Fidkowski.  Airfoil shape optimization using output-based adapted meshes.  In 23rd AIAA Computational Fluid Dynamics Conference, page 3102, 2017.
  * [38] Xiaolong He, Jichao Li, Charles A Mader, Anil Yildirim, and Joaquim RRA Martins.  Robust aerodynamic shape optimization—from a circle to an airfoil.  Aerospace Science and Technology, 87:48–61, 2019.
  * [39] Xiaojing Wu, Weiwei Zhang, Xuhao Peng, and Ziyi Wang.  Benchmark aerodynamic shape optimization with the pod-based cst airfoil parametric method.  Aerospace Science and Technology, 84:632–640, 2019.
  * [40] Md Tausif Akram and Man-Hoe Kim.  Cfd analysis and shape optimization of airfoils using class shape transformation and genetic algorithm—part i.  Applied Sciences, 11(9):3791, 2021.
  * [41] Xinshuai Zhang, Fangfang Xie, Tingwei Ji, Zaoxu Zhu, and Yao Zheng.  Multi-fidelity deep neural network surrogate model for aerodynamic shape optimization.  Computer Methods in Applied Mechanics and Engineering, 373:113485, 2021.
  * [42] SUN Zhiwei, WANG Chen, Yu Zheng, BAI Junqiang, LI Zheng, XIA Qiang, and FU Qiujun.  Non-intrusive reduced-order model for predicting transonic flow with varying geometries.  Chinese Journal of Aeronautics, 33(2):508–519, 2020.
  * [43] S Ashwin Renganathan, Romit Maulik, and Jai Ahuja.  Enhanced data efficiency using deep neural networks and gaussian processes for aerodynamic design optimization.  Aerospace Science and Technology, 111:106522, 2021.
  * [44] Xiaosong Du, Ping He, and Joaquim RRA Martins.  Rapid airfoil design optimization via neural networks-based parameterization and surrogate modeling.  Aerospace Science and Technology, 113:106701, 2021.
  * [45] Peng Liao, Wei Song, Peng Du, and Hang Zhao.  Multi-fidelity convolutional neural network surrogate model for aerodynamic optimization based on transfer learning.  Physics of Fluids, 33(12):127121, 2021.
  * [46] Jun Tao and Gang Sun.  Application of deep learning based multi-fidelity surrogate model to robust aerodynamic design optimization.  Aerospace Science and Technology, 92:722–737, 2019.
  * [47] Christophe Geuzaine and Jean-Francois Remacle.  Gmsh.
  * [48] C.D. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo, D. De Grazia, S. Yakovlev, J.-E. Lombard, D. Ekelschot, B. Jordi, H. Xu, Y. Mohamied, C. Eskilsson, B. Nelson, P. Vos, C. Biotto, R.M. Kirby, and S.J. Sherwin.  Nektar++: An open-source spectral/hp element framework.  Computer Physics Communications, 192:205–219, 2015.
  * [49] David Moxey, Chris D. Cantwell, Yan Bao, Andrea Cassinelli, Giacomo Castiglioni, Sehun Chun, Emilia Juda, Ehsan Kazemi, Kilian Lackhove, Julian Marcon, Gianmarco Mengaldo, Douglas Serson, Michael Turner, Hui Xu, Joaquim Peiró, Robert M. Kirby, and Spencer J. Sherwin.  Nektar++: Enhancing the capability and application of high-fidelity spectral/hp element methods.  Computer Physics Communications, 249:107110, 2020.
  * [50] Eastman N Jacobs, Kenneth E Ward, and Robert M Pinkerton.  The characteristics of 78 related airfoil sections from tests in the variable-density wind tunnel.  National Advisory Committee for Aeronautics, 1933.
  * [51] Onur Rauf Bingol and Adarsh Krishnamurthy.  NURBS-Python: An open-source object-oriented NURBS modeling framework in Python.  SoftwareX, 9:85–94, 2019.
  * [52] de Berg Mark, Cheong Otfried, van Kreveld Marc, and Overmars Mark.  Computational geometry algorithms and applications.  Spinger, 2008.
  * [53] Gianmarco Mengaldo, Daniele De Grazia, Freddie Witherden, Antony Farrington, Peter Vincent, Spencer Sherwin, and Joaquim Peiro.  A guide to the implementation of boundary conditions in compact high-order methods for compressible aerodynamics.  In 7th AIAA Theoretical Fluid Mechanics Conference, page 2923, 2014\.
  * [54] Tianping Chen and Hong Chen.  Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems.  IEEE Transactions on Neural Networks, 6(4):911–917, 1995.
  * [55] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Graph kernel network for partial differential equations, 2020.
  * [56] Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  Journal of Marchine Learning Research, 18:1–43, 2018.
  * [57] Brian Adams, William Bohnhoff, Keith Dalbey, Mohamed Ebeida, John Eddy, Michael Eldred, Russell Hooper, Patricia Hough, Kenneth Hu, John Jakeman, Mohammad Khalil, Kathryn Maupin, Jason A. Monschke, Elliott Ridgway, Ahmad Rushdi, Daniel Seidl, John Stephens, Laura Painton Swiler, Anh Tran, and Justin Winokur.  Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis: Version 6.16 user’s manual.
  * [58] D Finkel and Carl Tim Kelley.  Convergence analysis of the direct algorithm.  Technical report, North Carolina State University. Center for Research in Scientific Computation, 2004.
  * [59] Zongren Zou, Xuhui Meng, Apostolos F Psaros, and George Em Karniadakis.  Neuraluq: A comprehensive library for uncertainty quantification in neural differential equations and operators.  arXiv preprint arXiv:2208.11866, 2022.
  * [60] David A. Kopriva and John H. Kolias.  A conservative staggered-grid chebyshev multidomain method for compressible flows.  Journal of Computational Physics, 125(1):244–261, 1996.
  * [61] Ahmad Peyvan, Jonathan Komperda, Dongru Li, Zia Ghiasi, and Farzad Mashayek.  Flux reconstruction using jacobi correction functions in discontinuous spectral element method.  Journal of Computational Physics, 435:110261, 2021.
  * [62] Peter Josef Kunz.  Aerodynamics and design for ultra-low Reynolds number flight.  Stanford University, 2003.
  * [63] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.