11affiliationtext: Department of Mathematics, KTH Royal Institute of
Technology, Stockholm, Sweden22affiliationtext: Department of Mathematics and
Statistics, McMaster University, Hamilton, Ontario, Canada

# Adjoint-Based Enforcement of State Constraints in PDE Optimization Problems

Pritpal Matharu Email address for correspondence<EMAIL_ADDRESS>Bartosz
Protas

###### Abstract

This study demonstrates how the adjoint-based framework traditionally used to
compute gradients in PDE optimization problems can be extended to handle
general constraints on the state variables. This is accomplished by
constructing a projection of the gradient of the objective functional onto a
subspace tangent to the manifold defined by the constraint. This projection is
realized by solving an adjoint problem defined in terms of the same adjoint
operator as used in the system employed to determine the gradient, but with a
different forcing. We focus on the “optimize-then-discretize” paradigm in the
infinite-dimensional setting where the required regularity of both the
gradient and of the projection is ensured. The proposed approach is
illustrated with two examples: a simple test problem describing optimization
of heat transfer in one direction and a more involved problem where an optimal
closure is found for a turbulent flow described by the Navier-Stokes system in
two dimensions, both considered subject to different state constraints. The
accuracy of the gradients and projections computed by solving suitable adjoint
systems is carefully verified and the presented computational results show
that the solutions of the optimization problems obtained with the proposed
approach satisfy the state constraints with a good accuracy, although not
exactly.

Keywords: PDE Optimization; Adjoint Analysis; State Constraints; Heat
Transfer; Turbulence

##  1 Introduction

Numerous problems in modern science and engineering are cast in terms of
optimization of systems described by partial differential equations (PDEs),
often subject to complex constraints imposed on the control (decision)
variables. A standard framework for handling such constrained PDE optimization
problems is based on Lagrange multipliers [Luenberger, 1969; Tröltzsch, 2010]
where the objective functional is considered a function of both the state and
the control variable, the two related by suitable constraints. However, this
formalism can lead to challenging computational problems since the system of
the first-order optimality conditions has the form of a saddle-type problem
which is often difficult to solve numerically. As an alternative, one can
consider formulations based on the reduced objective functional depending on
the control variable only (i.e., where the state variable is not a separate
independent variable, but is considered a function of the control variable),
such that the optimization problem can be solved using some discrete form of
the gradient flow in which the constraints are built into the definition of
the space in which the gradient flow is constructed. An advantage of this
approach is that a saddle-type (min-max) problem is avoided and one can take
advantage of computational methods of unconstrained optimization. When dealing
with infinite-dimensional PDE optimization problems a key element of this
approach is determination of the gradient of the reduced objective functional
with respect to the control variable which can be conveniently computed by
solving a suitably-defined adjoint system [Lions, 1968; Gunzburger, 2003].
These methods are already quite well developed and there exists a large body
of literature on this topic, with applications in various areas [Bukshtynov
and Protas, 2013; Vergnault and Sagaut, 2014; Costanzo et al., 2022; Sirignano
et al., 2023]. However, methods based on discrete gradient flows can be
challenging to apply in the presence of complicated constraints, especially
when they involve both state and control variables as is the case in some
applications. In the present study we demonstrate how such constraints can be
approximately enforced using the solution of a suitable adjoint problem,
similar but distinct from the adjoint system used to determine the gradient of
the reduced objective functional. We focus on the “optimize-then-discretize”
approach where all elements of the algorithm are derived in the continuous
setting and only then discretized for the purpose of the numerical solution
[Gunzburger, 2003]. To the best of our knowledge, this is the first such
generalization of adjoint-based techniques in the solution of PDE optimization
problems.

To illustrate the ideas described above, we consider the state
$u\in\mathcal{X}$ and the control variable ${\varphi}\in\mathcal{Y}$ with
$\mathcal{X}$ and $\mathcal{Y}$ denoting the appropriate functional spaces, to
be specified below. The objective functional then is
$j\;:\;\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}$ and the PDE
constraint is given in terms of some function
${S}\;:\;\mathcal{X}\times\mathcal{Y}\rightarrow\mathcal{X}^{*}$ representing
the weak form of the PDE with $\mathcal{X}^{*}$ the dual space of
$\mathcal{X}$. In addition, we will also assume that the state and the control
variables are subject to $m\in\mathbb{N}^{+}$ scalar constraints given by the
function $c\;:\;\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}^{m}$ assumed
sufficiently smooth. Our constrained PDE optimization problem can then be
stated as

$\textrm{(A)}\hfill\qquad\qquad\qquad\begin{aligned}
&\min_{(u,{\varphi})\in\mathcal{X}\times\mathcal{Y}}j(u,{\varphi})\\\
&\textrm{subject to}:\quad\left\\{\begin{aligned} {S}(u,{\varphi})&=0,\\\
c(u,{\varphi})&=0\end{aligned}\right..\end{aligned}$

Introducing the Lagrange multipliers $\lambda\in\mathcal{X}$,
$\mu\in\mathbb{R}^{m}$ and the augmented objective functional, i.e., the
Lagrangian,
$\mathcal{L}(u,{\varphi},\lambda,\mu):=j(u,{\varphi})+\left\langle\lambda,{S}(u,{\varphi})\right\rangle_{\mathcal{X}\times\mathcal{X}^{*}}+\langle\mu,c(u,{\varphi})\rangle_{\mathbb{R}^{m}}$,
where $\langle\cdot,\cdot\rangle_{\mathcal{X}\times\mathcal{X}^{*}}$ is the
duality pairing between the spaces $\mathcal{X}$ and $\mathcal{X}^{*}$
[Luenberger, 1969], $\langle\cdot,\cdot\rangle_{\mathbb{R}^{m}}$ the inner
product in $\mathbb{R}^{m}$, and “:=” means “equal to by definition”, the
constrained problem (A) can be recast in the corresponding unconstrained form
as [Luenberger, 1969]

$\textrm{(B)}\hfill\qquad\qquad\qquad\min_{(u,{\varphi})\in\mathcal{X}\times\mathcal{Y}}\quad\max_{(\lambda,\mu)\in\mathcal{X}^{*}\times\mathbb{R}^{m}}\mathcal{L}(u,{\varphi},\lambda,\mu).$

Critical points of this problem are characterized by the appropriate first-
order optimality conditions. Besides handling the saddle-type (min-max) nature
of these critical points, an extra complication is the need to also
approximate the optimal Lagrange multipliers $\widecheck{\lambda}$ and
$\widecheck{\mu}$ in addition to the optimal state and control variables
$\widecheck{u}$ and $\widecheck{\phi}$ (hereafter, the symbol
$\widecheck{\cdot}$ will denote the optimal value of a variable).

Defining the reduced objective functional
$\mathcal{J}({\varphi}):=j(u({\varphi}),{\varphi})$, which assumes that for
each ${\varphi}$ the map $u=u({\varphi})$ is well defined in terms of the
solution of the PDE problem ${S}(u({\varphi}),{\varphi})=0$, Problem (A) can
be recast such that minimization is performed with respect to ${\varphi}$
only, i.e.,

$\textrm{(C)}\hfill\qquad\qquad\qquad\begin{aligned}
&\min_{{\varphi}\in\mathcal{Y}}\mathcal{J}({\varphi})\\\ &\textrm{subject
to}:\quad c(u({\varphi}),{\varphi})=0\end{aligned},$

where the constraint can be interpreted as defining a codimension-$m$ manifold
${\mathcal{M}}:=\left\\{{\varphi}\in\mathcal{Y}:\quad
c(u({\varphi}),{\varphi})=0\right\\}$. Thus, since the optimizer
${\widecheck{\varphi}}$ must be sought on this constraint manifold,
${\widecheck{\varphi}}\in{\mathcal{M}}$, problem (C) defines a Riemannian
optimization problem [Absil et al., 2008].

A local minimizer of Problem (C) can then be approximated with a discrete
projected gradient flow as
${\widecheck{\varphi}}=\lim_{n\rightarrow\infty}{\varphi}^{(n)}$, where

$\displaystyle{{\varphi}}^{(n+1)}$
$\displaystyle={{\varphi}}^{(n)}-{\tau_{n}}\,{P_{\mathcal{TM}_{{\varphi}^{(n)}}}\left(\nabla_{\varphi}\mathcal{J}({\varphi}^{(n)})\right)},\quad\quad
n=0,1,\dots,$ (1a) $\displaystyle{{\varphi}}^{(0)}$
$\displaystyle={{\varphi}}_{0},$ (1b)

in which
${P}_{\mathcal{TM}_{{\varphi}}}\>:\;\mathcal{Y}\rightarrow\mathcal{TM}_{{\varphi}}$
is the operator representing the orthogonal projection onto the subspace
$\mathcal{TM}_{{\varphi}}\in\mathcal{Y}$ tangent to the constraint manifold
$\mathcal{M}$ at ${\varphi}\in\mathcal{M}$, $\tau_{n}$ is the step size along
the projected gradient, $\nabla_{{\varphi}}\mathcal{J}$ is the gradient of the
reduced functional $\mathcal{J}(\varphi)$ with respect to ${\varphi}$, and
${\varphi}_{0}$ the initial guess. As illustrated in Figure 1, in iterations
(1) the constraint is satisfied approximately only, with an error
$\mathcal{O}(\tau_{n}^{2})$ at each iteration. This error is eliminated if one
can introduce the retraction operator
${\mathcal{R}_{\mathcal{M}}}\;:\>\mathcal{TM}\rightarrow\mathcal{M}$, where
$\mathcal{TM}$ is the tangent bundle, such that the right-hand side (RHS) in
(1a) is replaced with
${\mathcal{R}_{\mathcal{M}}\left({\varphi}^{(n)}-\tau_{n}\,P_{\mathcal{TM}_{{\varphi}^{(n)}}}\left(\nabla_{\varphi}\mathcal{J}({\varphi}^{(n)})\right)\right)}$
yielding the Riemannian gradient approach [Absil et al., 2008]. However, for
general manifolds $\mathcal{M}$ the retraction operator can be very difficult
to implement and as a result one often needs to resort to using (1). While the
gradient ${\nabla_{{\varphi}}\mathcal{J}({\varphi})}$ of the (reduced)
objective functional can be conveniently computed by solving an adjoint
problem [Gunzburger, 2003], the main contribution of the present study is to
demonstrate that an analogous computation can in fact be also used to
determine the projection operator ${P}_{\mathcal{TM}_{{\varphi}}}$. This
problem is nuanced by the fact that both the gradient
${\nabla_{{\varphi}}\mathcal{J}({\varphi})}$ and the projector
${P}_{\mathcal{TM}_{{\varphi}}}$ need to be determined with respect to a
particular topology, typically induced by the norm in the space $\mathcal{Y}$.
This contribution thus expands the scope of applications of adjoint analysis
in PDE constrained optimization.

We introduce and illustrate the new approach by way of two examples: the first
one is rather academic and concerns a simple control problem involving a one-
dimensional (1D) heat equation; the second is more involved and represents an
extension of a recently formulated problem concerning finding optimal
turbulence closures in two-dimensional (2D) Navier-Stokes flows [Matharu and
Protas, 2022]. While the first example is intended to serve as a simple
illustration only of the proposed approach, the second one showcases an
application to a nontrivial problem with a nonstandard structure. The
structure of the paper is as follows: in the next section we introduce the two
PDE optimization problems whereas the proposed approach is presented in
Section 3 with a particular emphasis on how the projection onto the subspace
tangent to the constraint manifold can be constructed based on the solution of
an adjoint system; computational results are then presented in Section 5 with
a summary and conclusions deferred to Section 6.

$\nabla_{{\varphi}}\mathcal{J}$$\mathcal{N}_{\varphi}$${\mathcal{TM}_{\varphi}}$${P_{\mathcal{TM}_{\varphi}}}$${P_{\mathcal{TM}_{\varphi}}}\left(\nabla_{{\varphi}}\mathcal{J}\right)$$\mathcal{R}_{\mathcal{M}}\left({P_{\mathcal{TM}_{\varphi}}}\left(\nabla_{{\varphi}}\mathcal{J}\right)\right)$$\mathcal{M}$$\varphi$$\varphi^{\text{new}}$
Figure 1: Schematic representation of the projection of the gradient
${\nabla_{\varphi}\mathcal{J}}$ onto the subspace $\mathcal{TM}_{{\varphi}}$
tangent to a general manifold $\mathcal{M}$ at ${\varphi}$. The projection
operator ${P_{\mathcal{TM}_{\varphi}}}$ uses the element
$\mathcal{N}_{\varphi}$ normal to the tangent subspace
$\mathcal{TM}_{{\varphi}}$, cf. 8ab. The retraction operator
$\mathcal{R}_{\mathcal{M}}$ then maps elements of the tangent subspace
$\mathcal{TM}_{{\varphi}}$ back to the manifold $\mathcal{M}$, cf. 8ac, to a
new element ${\varphi}^{\text{new}}$. The objects shown in the figure are to
be interpreted as elements of a function space, typically a Sobolev space
[Adams and Fournier, 2003], and can be functions of time as in Problems 1 and
2, or state as in Problem 3.

##  2 Test Problems

In this section we introduce two constrained PDE optimization problems that
will serve as our test cases. While the first problem is classical and is used
primarily to illustrate our approach, the second one is more challenging due
to its non-standard structure. With a slight abuse of notation motivated by
the desire to highlight the analogies between them, in the two test problems
we reuse the symbols denoting key elements of the formulation of the
optimization problems, with the exception of $\phi$ and $\varphi$ which
represent the control variable in Section 2.1 and Section 2.2, respectively.

###  2.1 Control of Heat Conduction in 1D

We consider heat conduction in a finite rod $\Omega:=[a,b]\subset\mathbb{R}$,
where $-\infty<a<b<\infty$, in which the time-dependent heat flux $\phi(t)$,
$t\in[0,T]$ with some $T>0$, is applied at the left boundary $x=a$ such that
the resulting temperature at the right boundary $x=b$ should match some
prescribed history $\bar{u}_{b}(t)$, $t\in[0,T]$, and at the same time the
average “energy” of the system should be given by $E_{0}>0$. The problem is
thus described by the system

$\displaystyle\frac{\partial u}{\partial t}-\Delta u$ $\displaystyle=0,$
$\displaystyle\qquad(t,x)$ $\displaystyle\in(0,T]\times\Omega,$ (2a)
$\displaystyle\frac{\partial u}{\partial x}\Big{|}_{x=a}$
$\displaystyle=\phi(t),$ $\displaystyle\quad t$ $\displaystyle\in(0,T],$ (2b)
$\displaystyle\frac{\partial u}{\partial x}\Big{|}_{x=b}$ $\displaystyle=0,$
$\displaystyle\quad t$ $\displaystyle\in(0,T],$ (2c) $\displaystyle u(t=0)$
$\displaystyle=u_{0},$ $\displaystyle\quad x$ $\displaystyle\in\Omega,$ (2d)

where $u_{0}\;:\;\Omega\rightarrow\mathbb{R}$ is the initial condition and the
solution $u=u(t,x;\phi)$ is assumed to depend on the flux $\phi$ as the
control variable. The energy $E$ is defined as

$E(t;\phi):=\frac{1}{2}\int_{\Omega}u(t,x;\phi)^{2}\,dx,$ (3)

whereas $[f]_{T}:=(1/T)\int_{0}^{T}f(t)\,dt$ will denote the time average of
some function $f\;:\;[0,T]\rightarrow\mathbb{R}$. Next, we introduce the
Sobolev space $H^{p}(\Lambda)$ with $p\in\mathbb{N}$ of functions defined on
some domain $\Lambda$ with $p$ square-integrable distributional derivatives
and endowed with the inner product

$\forall\ z_{1},z_{2}\in H^{p}(\Lambda)\qquad\left\langle
z_{1},z_{2}\right\rangle_{H^{p}(\Lambda)}{:=}\sum_{q=0}^{p}\ell_{q}^{2q}\left\langle\frac{d^{q}z_{1}}{dx^{q}},\frac{d^{q}z_{2}}{dx^{q}}\right\rangle_{L^{2}(\Lambda)},$
(4)

where $\langle
z_{1},z_{2}\rangle_{L^{2}(\Lambda)}:=\int_{\Lambda}z_{1}z_{2}\,dx$, whereas
$\ell_{q}\in\mathbb{R}^{+}$, $q=0,\dots,p$, are “length-scale” parameters (in
the multidimensional case, i.e., when $\Lambda\subseteq\mathbb{R}^{n}$, $n>1$,
the definition (4) admits a natural generalization [Adams and Fournier,
2003]). While for different values of $0<\ell_{q}<\infty$ the inner products
in (4) are equivalent (in the precise sense of norm equivalence), in Section 5
we will show that the flexibility represented by these parameters
significantly improves the convergence of the our iterative minimization
algorithm (1).

It is instructive to consider two formulations: one in which the flux is
simply assumed to be a square-integrable function of time together with its
first derivative, i.e., $\phi\in\mathcal{S}:=H^{1}(0,T)$, and another one
where it is additionally assumed to belong to a quadratic manifold,
$\phi\in\mathcal{M}\subset\mathcal{S}$, which is defined as

${\mathcal{M}}:=\left\\{\phi\in{\mathcal{S}}:\quad\left[E(\cdot;\phi)\right]_{T}=\frac{1}{2T}\int_{0}^{T}\int_{\Omega}u(t,x;\phi)^{2}\,dxdt=E_{0}\right\\}.$
(5)

Our (reduced) objective functional
${\mathcal{J}}\;:\;{\mathcal{M}}\rightarrow\mathbb{R}$ is a measure of the
error between the actual temperature at the right boundary $u(t,b;\phi)$ and
the prescribed profile $\bar{u}_{b}(t)$

$\displaystyle{\mathcal{J}}(\phi)=\frac{1}{2}\int_{0}^{T}[u(\phi)|_{b}-\bar{u}_{b}]^{2}\,dt.$
(6)

We note that the quadratic constraint defining the manifold ${\mathcal{M}}$ in
(5) simplifies when the governing system (2) is subject to the homogeneous
initial condition $u_{0}(x)=0$, $\forall\ x\in\Omega$ ($u_{0}\equiv 0$). In
such case the quadratic constraint becomes homogeneous, i.e.,
$\left[E(\cdot;{\beta}\phi)\right]_{T}={\beta}^{2}\left[E(\cdot;\phi)\right]_{T}$
for some $\beta\in\mathbb{R}$, which reduces enforcement of this constraint to
a simple rescaling. Otherwise, enforcement of the constraint is more involved
and necessitates the use of a retraction operator. For comparison, we thus
consider two versions of each of the optimization problems depending on
whether or not the initial condition $u_{0}$ is homogeneous, i.e.,

######  Problem 1.

Given system 2 with $u_{0}\equiv 0$ and objective functional 6, find

$\displaystyle=\underset{\phi\in{\mathcal{S}}}{\operatorname*{arg\,min}}\,{\mathcal{J}}(\phi),$
(Problem 1.A)
$\displaystyle=\underset{\phi\in{\mathcal{M}}}{\operatorname*{arg\,min}}\,{\mathcal{J}}(\phi).$
(Problem 1.B)

######  Problem 2.

Given system 2 with $u_{0}\not\equiv 0$ and objective functional 6, find

$\displaystyle=\underset{\phi\in{\mathcal{S}}}{\operatorname*{arg\,min}}\,{\mathcal{J}}(\phi),$
(Problem 2.A)
$\displaystyle=\underset{\phi\in{\mathcal{M}}}{\operatorname*{arg\,min}}\,{\mathcal{J}}(\phi).$
(Problem 2.B)

We add that since in Problem 1.A and Problem 2.A minimization is performed
over the entire linear subspace $\mathcal{S}$, we will consider these problems
as “unconstrained”. In contrast, Problem 1.B and Problem 2.B, where
optimization is carried out over the manifold $\mathcal{M}\subset\mathcal{S}$,
will be regarded as “constrained”.

###  2.2 Optimal Eddy Viscosity in Closure Models for 2D Turbulent Flows

Turbulent flows governed by the incompressible Navier-Stokes system remain
very challenging to compute when the Reynolds number is high due to
difficulties resolving motions occurring at small spatial and temporal scales.
A solution often adopted in practice relies on solving a filtered version of
the governing equations such that the number of resolved degrees of freedom is
reduced, an approach referred to as the Large-Eddy Simulation (LES) [Lesieur,
1993; Pope, 2000; Davidson, 2015]. However, the LES system is not closed as it
involves terms explicitly depending on the unresolved degrees of freedom and
closing this system in terms of the resolved degrees of freedom leads to the
celebrated turbulence closure problem. While a significant body of literature
has been devoted to deriving turbulence models adapted to different flows,
most of these approaches have been empirical in nature. Recent advances in
machine learning have enabled the development of data-based techniques for
deducing closure models [Kutz, 2017; Gamahara and Hattori, 2017; Jimenez,
2018; Duraisamy et al., 2019; Duraisamy, 2021; Pawar and San, 2021;
Waschkowski et al., 2022]. As an alternative to these approaches, a data-based
technique relying on calculus of variations and methods of PDE optimization
was recently proposed in [Matharu and Protas, 2020, 2022]. It allows one to
find an optimal, in a mathematically precise sense, functional form of the
eddy viscosity appearing in a commonly used family of turbulence closure
models. The problem considered here arises as an extension of that approach.

We consider the flow of viscous incompressible fluids on a two-dimensional
(2D) periodic domain $\Omega=[0,2\pi]^{2}$ and over a time interval $[0,T]$
for some $T>0$. Assuming uniform fluid density $\rho=1$, the motion is
governed by the Navier-Stokes system, which written in vorticity-
streamfunction form, is

$\displaystyle\partial_{t}{w}+\boldsymbol{\nabla}^{\perp}\psi\cdot\boldsymbol{\nabla}{w}$
$\displaystyle=\nu_{N}\boldsymbol{\Delta}{w}-\alpha{w}+f_{\omega}$
$\displaystyle\text{in}\quad{(0,T]\times\Omega},$ (7a)
$\displaystyle\boldsymbol{\Delta}\psi$ $\displaystyle=-{w}$
$\displaystyle\text{in}\quad{(0,T]\times\Omega},$ (7b) $\displaystyle{w}(t=0)$
$\displaystyle={w_{0}}$ $\displaystyle\text{in}\quad\Omega,$ (7c)

where $w:=-\boldsymbol{\nabla}^{\perp}\cdot\boldsymbol{u}$, with
$\boldsymbol{\nabla}^{\perp}:=[\partial_{x_{2}},-\partial_{x_{1}}]^{T}$ and
$\boldsymbol{u}$ the velocity field, is the vorticity component perpendicular
to the plane of motion, $\psi$ is the streamfunction, $\nu_{N}$ is the
coefficient of the kinematic viscosity (for simplicity of notation, we reserve
the symbol $\nu$ for the eddy viscosity), and ${w_{0}}$ is the initial
condition. In (7a) the term proportional to $\alpha>0$ represents large-scale
dissipation (Ekman friction), whereas $f_{\omega}$ is a time-independent band-
limited forcing acting on low wavenumbers. The parameters of these two terms
are chosen such that the flow is in a statistical equilibrium with a well-
developed enstrophy cascade and a rudimentary inverse energy cascade [Lesieur,
1993].

Using $\widetilde{(\cdot)}$ to denote a suitable low-pass filter, we can write
the LES version of 7 as (the reader is referred to [Matharu and Protas, 2022]
for derivation details)

$\displaystyle\partial_{t}\widetilde{\omega}+$
$\displaystyle=\boldsymbol{\nabla}\cdot{\stackon[-6pt]{$\left(\left[\nu_{N}+\nu(s)\right]\boldsymbol{\nabla}\widetilde{\omega}\right)$}{\vstretch{1.3}{\hstretch{3.0}{\widetilde{\phantom{\;\;\;\;\;\;\;\;\;}}}}}}-\alpha\widetilde{\omega}+f_{\omega}$
$\displaystyle\text{in}\quad{(0,T]\times\Omega},$ (8a)
$\displaystyle\boldsymbol{\Delta}\widetilde{\psi}$
$\displaystyle=-\widetilde{\omega}$
$\displaystyle\text{in}\quad{(0,T]\times\Omega},$ (8b)
$\displaystyle\widetilde{\omega}(t=0)$
$\displaystyle=\widetilde{\omega}_{0}{:=}\widetilde{w}_{0}$
$\displaystyle\text{in}\quad\Omega,$ (8c)

where $\widetilde{\omega}$ is the LES vorticity, and the initial condition
$\widetilde{\omega}_{0}$ is given as the filtered initial condition (7c) and
the forcing term is unaffected by the filter as it acts on the low wavenumbers
only. The LES equation (8a) features a Smagorinsky-type closure model with a
state-dependent eddy viscosity expressed as

$\nu(s)=\left[{\eta^{3}}\,\sqrt{s}+\nu_{0}\right]\varphi\left(\frac{s}{s_{\text{max}}}\right)\qquad\text{with}\quad
s:=|\boldsymbol{\nabla}\widetilde{\omega}|^{2}\in\mathcal{I}:=[0,s_{\text{max}}],$
(8i)

where $\nu_{0}>0$, $\eta=2\pi/k_{c}$ is the width of the LES filter with
$k_{c}$ the largest resolved wavenumber, $s_{\text{max}}>0$ is a sufficiently
large number to be specified later and
$\varphi\;:\;[0,1]\rightarrow\mathbb{R}$ a non-dimensional function. The form
of equation (8a) suggests that $\nu=\nu(s)$, and hence also
$\varphi=\varphi(s/s_{\text{max}})$, must be at least piecewise $C^{1}$
functions on $\mathcal{I}$ and $[0,1]$, respectively. However, as will become
evident in Section 3.3, our solution approach imposes some additional
regularity requirements, namely, $\nu=\nu(s)$ needs to be piecewise $C^{3}$ on
$\mathcal{I}$ with the first and third derivatives vanishing at
$s=0,s_{\text{max}}$. Since gradient-based solution approaches to PDE-
constrained optimization problems are preferably formulated in Hilbert spaces
[Protas et al., 2004], we shall look for an optimal function $\varphi$
parametrizing the eddy viscosity as an element of the following linear space
which is a subspace of the Sobolev space $H^{2}(\mathcal{I})$, cf (4),

$\mathcal{S}:={\left\\{\varphi\in{C^{3}([0,1])}\>:\>\frac{d}{d\xi}\varphi(\xi)=\frac{d^{3}}{d\xi^{3}}\varphi(\xi)=0\
\text{at}\ \xi=0,1\right\\}}.$ (8j)

We will seek an optimal form of the function $\varphi$ such that solutions of
the LES system (8) with eddy viscosity (8i) best match the corresponding
solutions of the original Navier-Stokes system (7) in a sense to be specified
below. We add that when $\nu_{0}=0$ and $\varphi(\xi)=1$, $\xi\in[0,1]$, then
the eddy viscosity (8i) reduces to the Leith model which is a counterpart of
the Smagorinsky model in 2D flows [Leith, 1968, 1971, 1996].

While in [Matharu and Protas, 2022] we considered a pointwise match, both in
space and in time, between the LES and the Navier-Stokes flows, here we
determine an optimal eddy viscosity $\widecheck{\nu}(s)$ that will allow the
LES flow to match the original Navier-Stokes flow in a certain average sense.
Important quantities characterizing 2D flows are the enstrophy and
palinstrophy defined as follows

$\displaystyle\mathcal{E}(t)$
$\displaystyle:=\frac{1}{2}\int_{\Omega}\widetilde{w}(t,\mathbf{x})^{2}\,d\mathbf{x},$
(8ka) $\displaystyle\mathcal{P}(t)$
$\displaystyle:=\frac{1}{2}\int_{\Omega}\left|\boldsymbol{\nabla}\widetilde{w}(t,\mathbf{x})\right|^{2}\,d\mathbf{x},$
(8kb)

where $\widetilde{w}$ is the filtered solution of the Navier-Stokes system
(7), whereas $\widetilde{\mathcal{E}}(t;\varphi)$ and
$\widetilde{\mathcal{P}}(t;\varphi)$ will denote the corresponding quantities
defined in terms of the solution $\widetilde{\omega}$ of the LES system (8)
and the notation emphasizes their dependence on the function $\varphi$
parameterizing the eddy viscosity, cf. ansatz (8i).

The optimal eddy viscosity $\widecheck{\nu}(s)$ will be chosen so as to yield
a good match between the palinstrophy in the LES flow (with the given eddy
viscosity) and in the original Navier-Stokes flow, so that we shall consider
the following error functional

$\mathcal{J}({{\varphi}}):=\frac{1}{2\,D}\int_{0}^{T}\left[\mathcal{P}(t)-\widetilde{\mathcal{P}}(t;{\varphi})\right]^{2}\,dt,$
(8l)

where $D>0$ serves as a normalization factor. In addition, we will also
require the LES flow with the optimal eddy viscosity to have a time-averaged
enstrophy equal to the enstrophy in the original Navier-Stokes flow, i.e.,

$\left[\widetilde{\mathcal{E}}(\cdot;{\varphi})\right]_{T}=\Big{[}\mathcal{E}(\cdot)\Big{]}_{T}=:\mathcal{E}_{0}.$
(8m)

This condition thus defines the following nonlinear manifold in the space
$\mathcal{S}$, cf. (8j), of nondimensional functions parameterizing the eddy
viscosity

${\mathcal{M}}:=\left\\{{\varphi}\in\mathcal{S}\>:{\left[\widetilde{\mathcal{E}}(\cdot;{\varphi})\right]_{T}}=\mathcal{E}_{0},\right\\}.$
(8n)

The optimal eddy viscosity $\widecheck{\nu}=\widecheck{\nu}(s)$ can then be
obtained from the solution of the following optimization problems via ansatz
(8i) with the function $\varphi$ replaced with $\widecheck{\varphi}$, where as
in the previous section we consider the unconstrained and constrained
formulation with and without the nonlinear constraint

######  Problem 3.

For the system 8 and objective functional 8l, find

$\displaystyle\widecheck{{\varphi}}:=$
$\displaystyle\underset{{\varphi}\in{\mathcal{S}}}{\operatorname*{arg\,min}}\,\mathcal{J}({\varphi}),$
(Problem 3.A) $\displaystyle\widecheck{{\varphi}}:=$
$\displaystyle\underset{{\varphi}\in{\mathcal{M}}}{\operatorname*{arg\,min}}\,\mathcal{J}({\varphi}).$
(Problem 3.B)

We emphasize that Problem 3 has a fundamentally different structure than
Problems 1 and 2, since the control variable $\varphi(s/s_{\text{max}})$ is a
function of the dependent (state) variable $s$, cf. (8i), rather than the
independent variables ($t$ or $x$) as is typical in PDE-constrained
optimization problems. In other words, solving Problem 3 can be interpreted as
finding an optimal form of the nonlinearity in the closure model. As described
in the next section, this aspect will have significant ramifications for how
Problem 3 is solved.

##  3 Solution Approach

In this section we first briefly present an adjoint-based approach to compute
the gradient $\nabla_{\phi}{\mathcal{J}}(\phi)$ of the (reduced) objective
functional in Problem 1 and 2, which is quite standard, and then describe how
an analogous approach can be used to define the projection operator
$P_{\mathcal{TM}_{\phi}}$ realizing the orthogonal projection onto the
subspace tangent to the constraint manifold ${\mathcal{M}}$, cf. (1a) and
Figure 1. Finally, we provide some details how these approaches can be adapted
to perform analogous tasks in the solution of Problem 3. Throughout these
derivations particular attention will be paid to ensuring the required
regularity of the obtained solutions which will be done using the framework
established in [Protas et al., 2004].

###  3.1 Evaluation of the Cost Functional Gradients

In order to determine the gradient $\nabla_{\phi}{\mathcal{J}}$ of the
objective functional (6), we begin by computing its Gâteaux (directional)
differential in the direction of some arbitrary perturbation $\phi^{\prime}\in
H^{1}(0,T)$ of the heat flux

${\mathcal{J}}^{\prime}(\phi;\phi^{\prime})=\frac{d}{d\epsilon}\,{\mathcal{J}}\left(\phi+\epsilon\phi^{\prime}\right)\Big{|}_{\epsilon=0}=\int_{0}^{T}[u(\phi)|_{b}-\bar{u}_{b}]\,u^{\prime}(x,t;\phi,\phi^{\prime})\,dt,$
(8o)

where $u^{\prime}(t,x;\phi,\phi^{\prime})$ satisfies the system obtained as a
perturbation of the governing system (2)

$\displaystyle{\mathcal{K}}u^{\prime}:=$ $\displaystyle\frac{\partial
u^{\prime}}{\partial t}-\Delta u^{\prime}$ $\displaystyle=0$
$\displaystyle\qquad(t,x)$ $\displaystyle\in(0,T]\times\Omega,$ (8pa)
$\displaystyle\frac{\partial u^{\prime}}{\partial x}\Big{|}_{x=a}$
$\displaystyle=\phi^{\prime}(t)$ $\displaystyle\quad t$
$\displaystyle\in(0,T],$ (8pb) $\displaystyle\frac{\partial
u^{\prime}}{\partial x}\Big{|}_{x=b}$ $\displaystyle=0$ $\displaystyle\quad t$
$\displaystyle\in(0,T],$ (8pc) $\displaystyle u^{\prime}(t=0)$
$\displaystyle=0$ $\displaystyle\quad x$ $\displaystyle\in\Omega.$ (8pd)

In order to extract the gradient $\nabla_{\phi}{\mathcal{J}}$ from the Gâteaux
differential (8o), we use the fact that the differential is a bounded linear
functional on both $L^{2}(0,T)$ and $H^{1}(0,T)$ when viewed as a function of
$\phi^{\prime}$, and invoke the Riesz representation theorem to obtain
[Berger, 1977]

${\mathcal{J}}^{\prime}({\phi};\phi^{\prime})=\left\langle{\nabla_{\phi}^{L^{2}}}{\mathcal{J}},\phi^{\prime}\right\rangle_{L^{2}{(0,T)}}=\Big{\langle}{\nabla_{\phi}}{\mathcal{J}},\phi^{\prime}\Big{\rangle}_{H^{1}{(0,T)}},$
(8q)

where the gradients ${\nabla_{\phi}^{L^{2}}}{\mathcal{J}}\in L^{2}(0,T)$ and
${\nabla_{\phi}}{\mathcal{J}}\in H^{1}(0,T)$ are the Riesz representers in the
two spaces. While in 1 we require the gradient in the space $H^{1}(0,T)$, it
is convenient to first obtain the gradient with respect to the $L^{2}$
topology. Since the expression on the RHS of (8o) is not consistent with the
Riesz form (8q) as the perturbation $\phi^{\prime}$ does not appear explicitly
in it, but is instead hidden in the boundary condition (8pb), we introduce the
adjoint field $u^{*}\;:\;[0,T]\times\Omega\rightarrow\mathbb{R}$ and define
the following duality-pairing relation

$\displaystyle\left({\mathcal{K}}u^{\prime},u^{*}\right):$
$\displaystyle=\int_{0}^{T}\int_{\Omega}{(\mathcal{K}u^{\prime})\,u^{*}}\,dx\,dt$
(8r)
$\displaystyle=\int_{0}^{T}\int_{\Omega}{u^{\prime}\,(\mathcal{K}^{*}u^{*}})\,dx\,dt-\overbrace{\int_{0}^{T}[u(\phi)|_{b}-\bar{u}_{b}]\,u^{\prime}(b,t;\phi,\phi^{\prime})\,dt}^{{{\mathcal{J}}}^{\prime}({{\phi}};{\phi}^{\prime})}-\int_{0}^{T}u^{*}\big{|}_{x=a}\,\phi^{\prime}(t)\,dt=0,$

where integration by parts was performed with respect to both space and time
and the field $u^{*}$ solves adjoint system

$\displaystyle{\mathcal{K}^{*}u^{*}}:=-$ $\displaystyle\frac{\partial
u^{*}}{\partial t}-\Delta u^{*}$ $\displaystyle=0$ $\displaystyle\qquad(t,x)$
$\displaystyle\in(0,T]\times\Omega,$ (8sa) $\displaystyle\frac{\partial
u^{*}}{\partial x}\Big{|}_{x=a}$ $\displaystyle=0$ $\displaystyle\quad t$
$\displaystyle\in(0,T],$ (8sb) $\displaystyle\frac{\partial u^{*}}{\partial
x}\Big{|}_{x=b}$ $\displaystyle=u(\phi)|_{b}-\bar{u}_{b}$ $\displaystyle\quad
t$ $\displaystyle\in(0,T],$ (8sc) $\displaystyle u^{*}(T=0)$ $\displaystyle=0$
$\displaystyle\quad x$ $\displaystyle\in\Omega.$ (8sd)

Collecting (8p), (8r) and (8s), we obtain an expression for the Gâteaux
differential consistent with the Riesz form (8q), namely,
${{\mathcal{J}}}^{\prime}({{\phi}};{\phi}^{\prime})=-\int_{0}^{T}u^{*}\big{|}_{x=a}\,\phi^{\prime}(t)\,dt$.
The gradient defined with respect to the $L^{2}$ topology is then deduced
using the first equality in (8q)

$\nabla_{\phi}^{L^{2}}{{\mathcal{J}}}(t)=-u^{*}(t,x)\big{|}_{x=a}.$ (8t)

The required Sobolev gradient ${\nabla_{\phi}{\mathcal{J}}}$ is then obtained
using the second equality in (8q) and the definition (4) which gives

${\mathcal{J}}^{\prime}(\phi;\phi^{\prime})=\int_{0}^{T}{\nabla_{\phi}^{L^{2}}{\mathcal{J}}}\,\phi^{\prime}\,dt=\int_{0}^{T}{\nabla_{\phi}{\mathcal{J}}}\,\phi^{\prime}\,dt+\ell^{2}\,\int_{0}^{T}\frac{d({\nabla_{\phi}{\mathcal{J}}})}{dt}\,\frac{d\phi^{\prime}}{dt}\,dt.$
(8u)

Performing integration by parts with respect to time $t$ in the second term
and assuming that both the perturbation and the gradient vanish at the
endpoints of the time interval, $\phi^{\prime}(0)=\phi^{\prime}(T)=0$ and
${\nabla_{\phi}{\mathcal{J}}}\big{|}_{t=0}={\nabla_{\phi}{\mathcal{J}}}\big{|}_{{t=T}}=0$,
we obtain the Sobolev gradient as the solution of the elliptic boundary-value
problem [Protas et al., 2004]

$\displaystyle\left[\operatorname{Id}-\ell^{2}\,\frac{d^{2}}{dt^{2}}\right]{\nabla_{\phi}{\mathcal{J}}}(t)$
$\displaystyle={\nabla_{\phi}^{L^{2}}{\mathcal{J}}}(t)\qquad t\in(0,T),$ (8va)
$\displaystyle{\nabla_{\phi}{\mathcal{J}}}\big{|}_{t=0}={\nabla_{\phi}{\mathcal{J}}}\big{|}_{t={T}}$
$\displaystyle=0.$ (8vb)

Given the form of the discrete gradient flow (1), this ensures that the
behavior of the control variable at the endpoints $t=0$ and $t=T$ can be
prescribed in the initial guess $\phi_{0}$ (in other words, with the Sobolev
gradients defined as in (8q) and (8v), the discrete gradient flow (1) does not
affect the boundary values of the control variable, which is desirable in some
applications). It can be shown [Protas et al., 2004] that extraction of
Sobolev gradients via the boundary-value problem (8v) can be interpreted as
application of a low-pass filter to the $L^{2}$ gradient with $\ell$ acting as
the cut-off parameter (it is thus a smoothing operation where Fourier
components of the gradient with wavelengths shorter than $\ell$ are damped).

###  3.2 Projection onto the Subspace Tangent to the Constraint Manifold

Projection onto the subspace ${\mathcal{TM}_{\phi}}$ tangent to the manifold
${\mathcal{M}}$ at a given element $\phi\in{\mathcal{M}}$ is a key step in
enforcing the associated constraint, cf. Figure 1. Since in the present
problem the codimension of the manifold is one, cf. (5), the tangent subspace
can be characterized in terms of an element ${\mathcal{N}_{\phi}}\in
H^{1}(0,T)$ normal to it in a suitable sense (if the codimension of the
manifold is higher, $\operatorname*{codim}({\mathcal{M}})>1$, then there would
be $\operatorname*{codim}({\mathcal{M}})$ such elements and the considerations
below generalize in a natural way). We note that the projection operator
$P_{\mathcal{TM}_{\phi}}$ in the discrete gradient flow (1) needs to be
defined with respect to the $H^{1}$ inner product, cf. (4), as this is the
topology of the ambient space. However, as we did above, we will first obtain
the normal element $\mathcal{N}_{\phi}^{L^{2}}$ defined with respect to the
$L^{2}$ inner product and then deduce its Sobolev counterpart, both as
functions of time.

We begin by considering the Gâteaux differential of the constraint
$\left[E(\cdot;\phi)\right]_{T}=E_{0}$, cf. (5), which yields

$\left[E^{\prime}(\cdot;\phi,\phi^{\prime})\right]_{T}:=\frac{d}{d\epsilon}\
\left[E(\cdot;\phi+\epsilon{\phi}^{\prime})\right]_{T}\bigg{|}_{\epsilon=0}=\frac{1}{T}\int_{0}^{T}\int_{\Omega}u(t,x;\phi)\,u^{\prime}(t,x;\phi,\phi^{\prime})\,dx\,dt=0,$
(8w)

where $u^{\prime}(t,x;\phi,\phi^{\prime})$ is the solution of the perturbation
system (8p). We then seek to express the directional differential (8w), which
is a bounded linear functional of $\phi^{\prime}$ on $L^{2}(0,T)$ and
$H^{1}(0,T)$, in terms of Riesz identities in these spaces, i.e.,

$\left[E^{\prime}(\cdot;\phi,\phi^{\prime})\right]_{T}=\left\langle{\mathcal{N}_{\phi}^{L^{2}}},{\phi}^{\prime}\right\rangle_{L^{2}{(0,T)}}=\Big{\langle}{\mathcal{N}_{\phi}},{\phi}^{\prime}\Big{\rangle}_{H^{1}{(0,T)}}=0.$
(8x)

Noting that the expression on the RHS of (8w) does not explicitly depend on
the perturbation $\phi^{\prime}$ which instead appears in the boundary
condition (8pb), we introduce a new adjoint state
$v^{*}\;:\;[0,T]\times\Omega\rightarrow\mathbb{R}$ that will satisfy the
following adjoint system

$\displaystyle{\mathcal{K}}^{*}v^{*}:=-$ $\displaystyle\frac{\partial
v^{*}}{\partial t}-\Delta v^{*}$ $\displaystyle=u(t,x;\phi)$
$\displaystyle\qquad(t,x)$ $\displaystyle\in(0,T]\times\Omega,$ (8ya)
$\displaystyle\frac{\partial v^{*}}{\partial x}\Big{|}_{x=a}$
$\displaystyle=0$ $\displaystyle\quad t$ $\displaystyle\in(0,T],$ (8yb)
$\displaystyle\frac{\partial v^{*}}{\partial x}\Big{|}_{x=b}$
$\displaystyle=0$ $\displaystyle\quad t$ $\displaystyle\in(0,T],$ (8yc)
$\displaystyle v^{*}(T=0)$ $\displaystyle=0$ $\displaystyle\quad x$
$\displaystyle\in\Omega.$ (8yd)

Then, introducing the duality relation

${\left({\mathcal{K}}u^{\prime},v^{*}\right):=\int_{0}^{T}\int_{\Omega}({\mathcal{K}}u^{\prime})v^{*}\,dxdt=\overbrace{\int_{0}^{T}\int_{\Omega}u^{\prime}({\mathcal{K}^{*}v^{*}})\,dx\,dt}^{\left[E^{\prime}(\cdot;\phi,\phi^{\prime})\right]_{T}}+\int_{0}^{T}v^{*}\big{|}_{x=a}\,\phi^{\prime}(t)\,dt=0,}$
(8z)

we conclude that
$\left[E^{\prime}(\cdot;\phi,\phi^{\prime})\right]_{T}=-\int_{0}^{T}v^{*}\big{|}_{x=a}\,\phi^{\prime}(t)\,dt$.
Using the first equality in (8x) provides an expression for the element normal
to the tangent subspace $\mathcal{TM}_{{\phi}}$ with respect to the $L^{2}$
inner product, namely,

${\mathcal{N}_{\phi}^{L^{2}}}(t)=-v^{*}(t,a),\qquad t\in[0,T].$ (8aa)

Using the second equality in (8x) and following the same steps as above, we
obtain ${\mathcal{N}_{\phi}}\in H^{1}(0,T)$, the element normal to
$\mathcal{TM}_{{\phi}}$ with respect to the $H^{1}$ topology, as the solution
of the elliptic boundary-value problem (8v) with the source term on the RHS in
(8va) replaced with expression (8aa). For a given $\phi\in\mathcal{M}$, the
projection operator then takes the form

${\forall\ {\Phi}\in
H^{1}(0,T)\qquad{P_{\mathcal{TM}_{\phi}}\Phi}:={\Phi-\zeta\,\mathcal{N}_{\phi}},\quad\text{where}\
\
{\zeta:=\frac{\left\langle\Phi,\mathcal{N}_{\phi}\right\rangle_{H^{1}(0,T)}}{\left\langle\mathcal{N}_{\phi},\mathcal{N}_{\phi}\right\rangle_{H^{1}(0,T)}}}}$
(8ab)

can be interpreted as the Lagrange multiplier associated with the condition
$\left\langle
P_{\mathcal{TM}_{\phi}}\Phi,\mathcal{N}_{\phi}\right\rangle_{H^{1}(0,T)}=0$.
We emphasize that each iteration of the gradient descent (1) requires the
solution of two adjoint systems: the solution of system (8s) allows us to
determine the gradient of the objective functional, whereas the solution of
(8y) is needed in order to construct the projection operator (8kb). The two
adjoint systems are defined in terms of the same operator $\mathcal{K}^{*}$,
but are subject to different boundary conditions and source terms.

As already indicated in Section 1, for general nonlinear manifolds
${\mathcal{M}}$ the projected discrete gradient flow (1) produces local
minimizers which do not satisfy the constraint exactly, but only with an error
of order $\mathcal{O}(\tau_{n}^{2})$ at each iteration. If the constraint is
homogeneous as is the case in Problem 1, then one can define the retraction
operator
$\mathcal{R}_{\mathcal{M}}\;:\>{\mathcal{TM}}\rightarrow{\mathcal{M}}$ in
terms of simple normalization

$\mathcal{R}_{{\mathcal{M}}}(\phi):=\sqrt{\frac{E_{0}}{\left[E(\cdot;\phi)\right]_{T}}}\
\phi.$ (8ac)

The resulting Riemannian gradient flow then satisfies the constraint exactly
(i.e., up to round-off errors only) [Absil et al., 2008].

###  3.3 Solution of 3

Here we provide some details describing how the framework introduced in
Sections 3.1 and 3.2 can be used to solve 3. Since the system considered in
this problem is more complicated, to balance completeness and brevity, we only
state the most important steps and refer the reader to [Matharu and Protas,
2022] for further details.

First, we compute the gradient of the objective functional (8l) and following
the same procedure as stated in Section 3.1 begin by determining the Gâteaux
differential of 8l

${\mathcal{J}^{\prime}({\varphi};{\varphi}^{\prime}):=\frac{d}{d\epsilon}\
\mathcal{J}\Big{(}{{\varphi}}+\epsilon{\varphi}^{\prime}\Big{)}\Big{|}_{\epsilon=0}=\frac{1}{D}\,\int_{0}^{T}\int_{\Omega}\left(\mathcal{P}(t)-\widetilde{\mathcal{P}}(t;{\varphi})\right)\boldsymbol{\Delta}\widetilde{\omega}(t,\mathbf{x};{\varphi})\,\widetilde{\omega}^{\prime}(t,\mathbf{x};{\varphi},{\varphi}^{\prime})\,d\mathbf{x}\,dt,}$
(8ad)

where ${\varphi}^{\prime}$ is an arbitrary perturbation of
${\varphi}\in{\mathcal{S}}$, and
$\widetilde{\omega}^{\prime}(t,\mathbf{x};{\varphi},{\varphi}^{\prime})$
satisfies the corresponding perturbation system obtained by linearizing the
governing system (8), namely,

$\displaystyle\mathcal{K}\begin{bmatrix}\widetilde{\omega}^{\prime}\\\ \\\
\widetilde{\psi}^{\prime}\end{bmatrix}:=\begin{bmatrix}\partial_{t}\widetilde{\omega}^{\prime}+\boldsymbol{\nabla}^{\perp}\widetilde{\psi}^{\prime}\cdot\boldsymbol{\nabla}\widetilde{\omega}+\boldsymbol{\nabla}^{\perp}\widetilde{\psi}\cdot\boldsymbol{\nabla}\widetilde{\omega}^{\prime}+\alpha\widetilde{\omega}^{\prime}\qquad\qquad\qquad\\\
-\boldsymbol{\nabla}\cdot\left(2(\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}^{\prime})(\frac{d\nu}{d{s}}\,\varphi\boldsymbol{\nabla}\widetilde{\omega}+\frac{\eta^{3}\,\sqrt{s}+\nu_{0}}{s_{\text{max}}}\,\frac{d\varphi}{d\sigma}\,{\boldsymbol{\nabla}\widetilde{\omega}})+(\nu_{N}+\nu)\boldsymbol{\nabla}\widetilde{\omega}^{\prime}\right)\\\
\boldsymbol{\Delta}\widetilde{\psi}^{\prime}+\widetilde{\omega}^{\prime}\end{bmatrix}$
(8aea)
$\displaystyle\phantom{\mathcal{K}\begin{bmatrix}\widetilde{\omega}^{\prime}\\\
\\\
\widetilde{\psi}^{\prime}\end{bmatrix}:}=\begin{bmatrix}\boldsymbol{\nabla}\cdot\left(\,(\eta^{3}\,\sqrt{s}+\nu_{0}){\varphi}^{\prime}\boldsymbol{\nabla}\widetilde{\omega}\,\right)\\\
\\\ 0\end{bmatrix},$ (8aeb)
$\displaystyle\widetilde{\omega}^{\prime}({t=0,\mathbf{x}})=0,$ (8aec)

where for ease of notation, we have denoted $\sigma:=s/s_{\text{max}}$.
Defining the adjoint system as

$\displaystyle\mathcal{K}^{*}\begin{bmatrix}\widetilde{\omega}^{*}\\\ \\\
\widetilde{\psi}^{*}\end{bmatrix}:=\begin{bmatrix}-\partial_{t}\widetilde{\omega}^{*}-\boldsymbol{\nabla}^{\perp}\widetilde{\psi}\cdot\boldsymbol{\nabla}\widetilde{\omega}^{*}+\alpha\widetilde{\omega}^{*}+\widetilde{\psi}^{*}\qquad\qquad\qquad\qquad\\\
-\boldsymbol{\nabla}\cdot\left(2\,(\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}^{*})\,(\frac{d\nu}{d{s}}\,\varphi\boldsymbol{\nabla}\widetilde{\omega}+\frac{\eta^{3}\,\sqrt{s}+\nu_{0}}{s_{\text{max}}}\,\frac{d\varphi}{d\sigma}\,{\boldsymbol{\nabla}\widetilde{\omega}})+(\nu_{N}+\nu)\boldsymbol{\nabla}\widetilde{\omega}^{*}\right)\\\
\boldsymbol{\Delta}\widetilde{\psi}^{*}-\boldsymbol{\nabla}^{\perp}\cdot(\widetilde{\omega}^{*}\,\boldsymbol{\nabla}\widetilde{\omega})\end{bmatrix}=\begin{bmatrix}W\\\
\\\ 0\end{bmatrix},$ (8afa)
$\displaystyle\widetilde{\omega}^{*}({t=T,\mathbf{x}})=0,$ (8afb)

with source term
$W(t,\mathbf{x}):=\frac{1}{D}\,\left(\mathcal{P}(t)-\widetilde{\mathcal{P}}(t;{\varphi})\right)\boldsymbol{\Delta}\widetilde{\omega}(t,\mathbf{x};{\varphi})$,
we obtain the duality-pairing relation

$\displaystyle\left(\mathcal{K}\begin{bmatrix}\widetilde{\omega}^{\prime}\\\
\widetilde{\psi}^{\prime}\end{bmatrix},\begin{bmatrix}\widetilde{\omega}^{*}\\\
\psi^{*}\end{bmatrix}\right):=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\mathcal{K}\begin{bmatrix}\widetilde{\omega}^{\prime}\\\
\widetilde{\psi}^{\prime}\end{bmatrix}\cdot\begin{bmatrix}\widetilde{\omega}^{*}\\\
\widetilde{\psi}^{*}\end{bmatrix}\,d\mathbf{x}\,dt$ (8ag) $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\begin{bmatrix}\widetilde{\omega}^{\prime}\\\
\widetilde{\psi}^{\prime}\end{bmatrix}\cdot\mathcal{K}^{*}\begin{bmatrix}\widetilde{\omega}^{*}\\\
\widetilde{\psi}^{*}\end{bmatrix}\,d\mathbf{x}\,dt=\left(\begin{bmatrix}\widetilde{\omega}^{\prime}\\\
\widetilde{\psi}^{\prime}\end{bmatrix},\mathcal{K}^{*}\begin{bmatrix}\widetilde{\omega}^{*}\\\
\psi^{*}\end{bmatrix}\right),$

where integration by parts was performed with respect to both space and time
and all boundary terms vanish due to periodic boundary conditions. This
relation together with the adjoint system (8af) allow us to re-express the
Gâteaux differential 8ad as

$\displaystyle\mathcal{J}^{\prime}({\varphi};{\varphi}^{\prime})$
$\displaystyle=\phantom{-}\int_{0}^{T}\int_{\Omega}W(t,\mathbf{x})\,\widetilde{\omega}^{\prime}(t,\mathbf{x};{\varphi},{\varphi}^{\prime})\,d\mathbf{x}\,dt,$
$\displaystyle=-\int_{0}^{T}\int_{\Omega}\,\left(\eta^{3}\,\sqrt{s}+\nu_{0}\right)\,\left(\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}^{*}\right)\,{\varphi}^{\prime}\,d\mathbf{x}\,dt,$
$\displaystyle=\int_{0}^{1}\left[-\int_{0}^{T}\int_{\Omega}\,\delta\left(\frac{\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}}{s_{\text{max}}}-\sigma\right)\,\left(\eta^{3}\,\sqrt{s}+\nu_{0}\right)\,\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}^{*}\,d\mathbf{x}\,dt\right]\,\varphi^{\prime}(\sigma)\,d\sigma,$
(8ah)

where the substitution
$\varphi^{\prime}(\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega})=\int_{0}^{1}\,\delta\left(\frac{\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}}{s_{\text{max}}}-\sigma\right)\,\varphi^{\prime}(\sigma)\,d\sigma$
with $\delta(\cdot)$ denoting the Dirac delta distribution was made and
Fubini’s theorem was used to swap the order of integration. These last two
steps are needed in order to change the integrations variables in (8ah) from
$d\mathbf{x}\,dt$ to $ds$, such that the Gâteaux differential 8ad can be
expressed in the required Riesz form

$\displaystyle\mathcal{J}^{\prime}(\varphi;\varphi^{\prime})=\Big{\langle}\nabla_{\varphi}\mathcal{J},\varphi^{\prime}\Big{\rangle}_{H^{2}({[0,1]})}$
$\displaystyle=\left\langle\nabla_{\varphi}^{L^{2}}\mathcal{J},{\varphi}^{\prime}\right\rangle_{L^{2}({[0,1]})}.$
(8ai)

This then allows us to extract the gradient with respect to the $L^{2}$
topology as

$\nabla_{\varphi}^{L^{2}}{\mathcal{J}}(\sigma)=-\int_{0}^{T}\int_{\Omega}\,\delta\left(\frac{\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}}{s_{\text{max}}}-\sigma\right)\,\left(\eta^{3}\,\sqrt{s}+\nu_{0}\right)\,\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}^{*}\,d\mathbf{x}\,dt,\qquad{\sigma\in[0,1].}$
(8aj)

As was done in Section 3.1, we use the Riesz representation in the Sobolev
space $H^{2}(0,T)$, cf. (8x), and the definition of the subspace $\mathcal{S}$
in (8j) which after integration by parts with respect to $\sigma$, cf. (8u),
give the Sobolev gradient $\nabla_{\varphi}\mathcal{J}$ as the solution of the
following elliptic boundary-value problem

$\displaystyle\left[\operatorname{Id}-\ell_{1}^{2}\,\frac{d^{2}}{d\sigma^{2}}+\ell_{2}^{4}\frac{d^{4}}{d\sigma^{4}}\right]\nabla_{\varphi}\mathcal{J}(\sigma)$
$\displaystyle=\nabla_{\varphi}^{L^{2}}\mathcal{J}(\sigma),\qquad\sigma\in[0,1],$
(8aka)
$\displaystyle\frac{d^{(1)}\,(\nabla_{\varphi}\mathcal{J})}{d\sigma^{(1)}}\Big{|}_{\sigma=0,1}$
$\displaystyle=\frac{d^{(3)}\,(\nabla_{\varphi}\mathcal{J})}{d\sigma^{(3)}}\Big{|}_{\sigma=0,1}=0.$
(8akb)

We add that the boundary conditions in (8akb), see also (8j), were selected in
order to ensure the vanishing of the terms at $\sigma=0,1$ which result from
integration by parts.

Similarly, we determine the normal element $\mathcal{N}_{\varphi}$ by
considering the Gâteaux differential of the constraint 8m with respect to
$\varphi$ and invoking the Riesz representation theorem, i.e.,

$\displaystyle\left[\widetilde{\mathcal{E}}^{\prime}(\cdot;{\varphi},{\varphi}^{\prime})\right]_{T}:=\frac{d}{d\epsilon}\
\left[\widetilde{\mathcal{E}}(\cdot;{\varphi}+\epsilon{\varphi}^{\prime})\right]_{T}\Bigg{|}_{\epsilon=0}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\int_{\Omega}\widetilde{\omega}(t,\mathbf{x};{\varphi})\,\widetilde{\omega}^{\prime}(t,\mathbf{x};{\varphi},{\varphi}^{\prime})\,d\mathbf{x}\,dt,$
$\displaystyle=\Big{\langle}\mathcal{N}_{\varphi},\varphi^{\prime}\Big{\rangle}_{H^{2}({[0,1]})}=\left\langle\mathcal{N}_{\varphi}^{L^{2}}\mathcal{J},{\varphi}^{\prime}\right\rangle_{L^{2}({[0,1]})}=0.$
(8al)

Introducing the new adjoint fields $\breve{{\omega}}^{*}$ and
$\breve{\psi}^{*}$, which satisfy the adjoint system 8af, with the source term
$W(t,\mathbf{x}):=\frac{1}{T}\,\widetilde{\omega}(t,\mathbf{x};{\varphi})$, we
can use the duality pairing 8ag to conclude that

$\displaystyle\left[\widetilde{\mathcal{E}}^{\prime}(\cdot;{\varphi},{\varphi}^{\prime})\right]_{T}$
$\displaystyle=-\int_{0}^{T}\int_{\Omega}\,\left(\eta^{3}\,\sqrt{s}+\nu_{0}\right)\,\left(\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\breve{\omega}^{*}\right)\,{\varphi}^{\prime}\,d\mathbf{x}\,dt,$
$\displaystyle=\int_{0}^{1}\left[-\int_{0}^{T}\int_{\Omega}\,\delta\left(\frac{\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}}{s_{\text{max}}}-\sigma\right)\,\left(\eta^{3}\,\sqrt{s}+\nu_{0}\right)\,\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\breve{\omega}^{*}\,d\mathbf{x}\,dt\right]\,\varphi^{\prime}(\sigma)\,d\sigma.$

Then, using the second equality in 8al, we obtain an expression for the
element normal to the subspace $\mathcal{TM}_{\varphi}$ with respect to the
$L^{2}$ topology as

$\mathcal{N}_{\varphi}^{L^{2}}(\sigma)=-\int_{0}^{T}\int_{\Omega}\,\delta\left(\frac{\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\widetilde{\omega}}{s_{\text{max}}}-\sigma\right)\,\left(\eta^{3}\,\sqrt{s}+\nu_{0}\right)\,\boldsymbol{\nabla}\widetilde{\omega}\cdot\boldsymbol{\nabla}\breve{\omega}^{*}\,d\mathbf{x}\,dt,\qquad{\sigma\in[0,1].}$
(8am)

Finally, employing the last equality in (8al) and performing analogous steps
as described above, we obtain the normal element $\mathcal{N}_{\varphi}\in
H^{2}(0,1)$ as a solution of the boundary-value problem 8ak, but with the
source term in (8aka) replaced with the RHS of 8am. This then allows us to
construct the projection operator $P_{\mathcal{TM}_{\varphi}}$ in the form
(8ab).

##  4 Numerical Discretization and Validation

In this section we briefly describe the numerical discretizations employed to
solve Problems 1, 2 and 3 using the discrete gradient flow (1). Then we
validate key elements of these computations requiring solution of adjoint
systems, namely, evaluation of the cost functional gradients (8t) and (8aj)
and the normal elements (8aa) and (8am). The proposed approach has been
implemented in MATLAB and the code is available in [Matharu, ] (it can be used
to generate the figures shown in Sections 4.1 and 5.1).

###  4.1 Problems 1 and 2

The governing system (2) and the adjoint systems (8s) and (8y) are discretized
in space with a standard second-order finite-difference scheme and integrated
in time using a second-order Crank-Nicolson scheme.

In order to validate the computation of the gradient (8t) and of the normal
element (8aa) we define the following quantities

$\displaystyle\kappa_{1}(\epsilon)$
$\displaystyle:=\frac{\epsilon^{-1}\left[{\mathcal{J}}({\phi}+\epsilon\phi^{\prime})-{\mathcal{J}}({\phi})\right]}{\Big{\langle}{\nabla_{\phi}^{L^{2}}}{\mathcal{J}}(\phi),\phi^{\prime}\Big{\rangle}_{L^{2}(0,T)}},\qquad\epsilon>0$
(8ana) $\displaystyle\kappa_{2}(\epsilon)$
$\displaystyle:=\frac{\epsilon^{-1}\left[E(\cdot;{\phi}+\epsilon\phi^{\prime})-E(\cdot;{\phi})\right]_{T}}{\Big{\langle}{\mathcal{N}_{\phi}^{L^{2}}},\phi^{\prime}\Big{\rangle}_{L^{2}(0,T)}}{,}$
(8anb)

in which the numerators represent first-order finite-difference approximations
of the Gâteaux differentials (8o) and (8w), whereas the denominators are the
corresponding Riesz forms (8q) and (8x). Clearly, we expect that
$\kappa_{i}(\epsilon)\approx 1$, $i=1,2$, and any deviations of these
quantities from unity are a result of errors in the numerical solution of the
PDE systems (2), (8s) and (8y), and in the evaluation of the integrals in (5)
and (6), truncation errors in the finite-difference formula in (8ana)–(8anb)
as well as round-off errors. It is expected that $\kappa_{1}(\epsilon)$ and
$\kappa_{2}(\epsilon)$ should approach unity for a certain range of values of
$\epsilon$ as the numerical discretization parameters are refined.

To fix attention, we focus here on the more general Problem 2 and add that
very similar results were also obtained for Problem 1. In our tests we set
$a=0$, $b=1$, $T=1$ and

$u_{0}(x)=10\,\cos(\pi\,x),\quad
x\in[0,1],\qquad\qquad{\bar{\phi}(t)=e^{-4}\,t\,\left[18-1000\,\cos\left(\frac{15\pi}{2}\,t\right)\right]},\quad
t\in[0,T],$ (8ao)

where $\bar{\phi}(t)$ is the “true” flux defining, via the solution of system
(2), the target profile $\bar{u}_{b}(t)$ appearing in the error functional
(6). In the diagnostic quantities (8ana)–(8anb) we use

${\phi_{0}:=\phi(t)=18\,\sin\left(\frac{\pi}{2}\,t\right)e^{-4t}},\qquad\phi^{\prime}(t)=4t\,e^{(-4+\pi)t},\qquad
t\in[0,T],$ (8ap)

representing, respectively, the “point” (in the function space $H^{1}(0,T)$)
where the Gâteaux differentials are computed and the direction.

Figure 2: [Problem 2] Dependence of 2 $|1-{\kappa_{1}(\epsilon)}|$ and 2
$|1-{\kappa_{2}(\epsilon)}|$ on $\epsilon$ for different spatial and temporal
discretizations $\Delta x$ and $\Delta t$.

In Figures 2a and 2b we show the dependence of the quantities
$|1-{\kappa_{i}(\epsilon)}|$, $i=1,2$, on $\epsilon$ for different indicated
values of the space and time discretization parameters $\Delta x$ and $\Delta
t$. We observe that $\kappa_{1}(\epsilon)$ and $\kappa_{2}(\epsilon)$ deviate
from the unity for very small and very large values of $\epsilon$ which is due
to, respectively, the subtractive cancellation (round-off) errors and the
truncation errors in the finite-difference formula in (8ana)–(8anb), both of
which are well-known effects. However, we also observe that for intermediate
values of $\epsilon$ spanning several orders of magnitude both
$\kappa_{1}(\epsilon)$ and $\kappa_{2}(\epsilon)$ exhibit plateaus with values
decreasing towards unity as the numerical parameters $\Delta x$ and $\Delta t$
are refined, which demonstrates the convergence of the approximations of the
cost functional gradient and the normal element. In the computations presented
in Section 5.1 the Sobolev gradients and normal directions are computed using
the inner product (4) with $\ell_{1}=0.01$.

###  4.2 3

The Navier-Stokes system (7), the corresponding LES system (8) and the adjoint
systems (8af) are approximated using a standard Fourier pseudo-spectral method
in space and using a third-order IMEX scheme introduced in [Alimo et al.,
2021] in time. Evaluation of the closure terms involving the state-dependent
eddy viscosity (8i) requires interpolation from the physical space $\Omega$ to
the state space $\mathcal{I}$ and differentiation with respect to the state
variable $s$, steps which are performed using methods based on Chebyshev
polynomials [Trefethen, 2013]. We refer the reader to [Matharu and Protas,
2022] for further numerical details.

The “target” Navier-Stokes flow $w$ is defined as the solution of system (7)
obtained with the parameters $\nu_{N}=4\times 10^{-4}$, $\alpha=5\times
10^{-3}$, $F=2$, $k_{a}=k_{b}=4$ over the time window $T=50\approx 20t_{e}$,
where $t_{e}:=\left[\int_{0}^{T}\mathcal{E}(t)\,dt/(8\pi^{2}T)\right]^{-1/2}$
is the eddy turnover time [Bracco and McWilliams, 2010]. The direct numerical
simulation (DNS) of 7 and the LES computation solving 8 use, respectively,
$N_{x}=256$ and $N_{x}=64$ equispaced grid points in each spatial direction.
They also use $N_{s}=128$ Chebyshev points to discretize the state space
$\mathcal{I}$. The LES system (8) is obtained using a “box” filter with the
cutoff $k_{c}=8$ [Matharu and Protas, 2022]. The use of this rather
“aggressive” filter is dictated by the desire to define a problem where the
presence of the constraint (8m) has a significant effect on the solutions of
Problem 3, which will in turn allow us to more easily elucidate the properties
of the proposed approach.

In order to validate the computation of the cost functional gradient and the
normal element, we define the following diagnostic quantities analogous to
(8ana)–(8anb)

$\displaystyle\kappa_{3}(\epsilon)$
$\displaystyle:=\frac{\epsilon^{-1}\left[\mathcal{J}({\varphi}+\epsilon\varphi^{\prime})-\mathcal{J}({\varphi})\right]}{\Big{\langle}{\nabla_{\varphi}^{L^{2}}\mathcal{J}(\varphi)},\varphi^{\prime}\Big{\rangle}_{L^{2}(0,1)}},\qquad\epsilon>0$
(8aqa) $\displaystyle\kappa_{4}(\epsilon)$
$\displaystyle:=\frac{\epsilon^{-1}\left[\widetilde{\mathcal{E}}(\cdot;{\varphi}+\epsilon{\varphi}^{\prime})-\widetilde{\mathcal{E}}(\cdot;{\varphi})\right]_{T}}{\Big{\langle}{\mathcal{N}^{L^{2}}_{{\varphi}}},{\varphi}^{\prime}\Big{\rangle}_{L^{2}(0,1)}},$
(8aqb)

where $\varphi=\varphi(s)$ and $\varphi^{\prime}=\varphi^{\prime}(s)$ are the
nondimensional functions corresponding, via ansatz (8i), to the eddy viscosity
in the Leith model [Leith, 1968, 1971, 1996]

$\nu_{L}(s)=(C_{L}\,k_{c})^{3}\,\sqrt{s}$ (8ar)

with the constant $C_{L}=4.702\times 10^{-3}$ chosen to ensure that
$\varphi\in{\mathcal{M}}$, and to the perturbations

$\displaystyle\nu^{\prime}(s)$
$\displaystyle=\nu_{L}(s),\qquad\qquad\qquad\text{with}\quad C_{L}=4.2\times
10^{-3},$ (8asa) $\displaystyle\nu^{\prime}(s)$
$\displaystyle=16\,\cos^{2}\left(\frac{s}{7.5}\right),$ (8asb)
$\displaystyle\nu^{\prime}(s)$
$\displaystyle=\exp\left(\frac{-2s}{30}\right).$ (8asc)

The quantities $|1-\kappa_{3}(\epsilon)|$ and $|1-\kappa_{4}(\epsilon)|$ are
shown as functions of $\epsilon$ for different perturbations
$\varphi^{\prime}$ and two different time steps $\Delta t$ in Figures 3a and
3b, respectively. Since the spatial discretization we use is spectrally
accurate (i.e., it has exponentially small errors), we focus here on assessing
the effect of the time discretization, which has an algebraic only rate of
convergence, on the accuracy of the gradient and the normal element. In
Figures 3a and 3b we see the expected behavior of the quantities
$\kappa_{3}(\epsilon)$ and $\kappa_{4}(\epsilon)$. In particular, they
approach unity over as broad range of intermediate values of $\epsilon$ as the
time step $\Delta t$ is refined, which confirms that the approximations of the
cost functional gradient and of the normal element are convergent. In the
computations presented in Section 5.2 the Sobolev gradients and normal
directions are computed using the inner product (4) with $\ell_{1}=1000$ and
$\ell_{2}=100$.

Figure 3: [Problem 3:] Dependence of (a) $|1-{\kappa_{3}(\epsilon)}|$ and (b)
$|1-{\kappa_{4}(\epsilon)}|$ on $\epsilon$ for different temporal
discretizations with $\Delta t=1\times 10^{-2}$ (empty symbols) and $\Delta
t=5\times 10^{-3}$ (filled symbols). The results are shown for the
perturbations $\varphi^{\prime}$ corresponding to $\nu^{\prime}$ given in 8asa
(blue circles), 8asb (red triangles), and 8asc (yellow squares).

##  5 Results

In this section we study solutions to Problems 1, 2 and 3 focusing our
attention on how well and how efficiently the constraints involved are
satisfied. In iterations (1) the step size $\tau_{n}$ is determined using
Brent’s algorithm to solve the line-minimization problem [Press et al., 2007].
In the unconstrained versions of the problems we use the Polak-Ribière version
of the conjugate gradient method [Nocedal and Wright, 2006] to accelerate
convergence of iterations. In the constrained versions of the problems when no
retraction operator is available, we limit the magnitude of the step size
$\tau_{n}$ to reduce drift away from the constraint manifold $\mathcal{M}$.

###  5.1 Problems 1 and 2

In Problems 1 and 2 we attempt to reconstruct the flux $\phi(t)$ at the left
boundary of the domain by minimizing the error functional (6) such that the
temperature at the right boundary matches the prescribed target temperature
$\bar{u}_{b}$ over the time window $[0,T]$; the latter temperature is obtained
by solving problem (2) with a certain “true” flux $\bar{\phi}(t)$, given in
(8ao), applied at the left boundary. The true flux $\bar{\phi}(t)$ thus
represents the “exact” solution of Problems 1.A and 2.A in the sense that
$\mathcal{J}(\bar{\phi})=0$ in these cases. However, since it does not satisfy
the the nonlinear constraint, $\bar{\phi}\notin\mathcal{M}$, the true flux
cannot be a solution of Problems 1.B and 2.B. In the solution of Problem 1.B
the nonlinear constraint, cf. (5), is enforced exactly using the retraction
operator $\mathcal{R}_{\mathcal{M}}$ whereas in the solution of Problem 2.B
this has to be done approximately using the orthogonal projection on the
tangent subspace (8ab).

The dependence of the normalized error functional
${\mathcal{J}}(\phi^{(n)})/{\mathcal{J}}(\phi_{0})$ and of the normalized
constraint $\left[E(\cdot;\phi^{(n)})\right]_{T}/E_{0}$ on the iteration count
$n$ in the solution of Problems 1.A and 1.B is shown, respectively, in Figures
4 and 4. The corresponding optimal fluxes $\widecheck{\phi}(t)$ are shown as
functions of time $t\in[0,T]$ in Figure 5. Figure 4 shows that, as expected,
in the solution of Problem 1.A the error functional (6) is reduced to a very
low level and this error is significantly higher in the solution of Problem
1.B where an exact solution is not expected to exist. On the other hand, the
nonlinear constraint in (5) is satisfied exactly in Problem 1.B, cf. Figure 4.
In contrast, in the unconstrained problem we see that the difference
$\left[E(\cdot;\phi)\right]_{T}-E_{0}$ becomes large, on the order of $E_{0}$,
after a few iterations, which is not be surprising as this constraint is not
imposed in Problem 1.A. As is evident from Figure 5, the true flux is
reconstructed rather well in Problem 1.A, except for the final times where a
large oscillation is present. The oscillatory character of the optimal flux
$\widecheck{\phi}(t)$ in this case indicates that in its original form Problem
1 is ill-posed. Methods for dealing with this issue are well developed and are
usually based on Tikhonov regularization [Engl et al., 1996] or formulation of
the problem on spaces of sufficiently regular functions [Protas et al., 2004].
Since this issue is only tangential to the main topic of the present study, we
do not consider it further here.

Moving on to discuss Problem 2, the dependence of the normalized error
functional ${\mathcal{J}}(\phi^{(n)})/{\mathcal{J}}(\phi_{0})$ and of the
normalized constraint $\left[E(\cdot;\phi^{(n)})\right]_{T}/E_{0}$ on the
iteration count $n$ in shown in Figures 6 and 6. We emphasize that now the
retraction operator $\mathcal{R}_{\mathcal{M}}$ is no longer available and the
nonlinear constraint in Problem 2.B can only be enforced approximately via
projection on the tangent subspace $\mathcal{TM}_{\phi}$, cf. Figure 1. We can
see that, as a result, the constraint is no longer satisfied exactly, but the
drift away from the constraint manifold $\mathcal{M}$ is slow, much less rapid
than when the constraint is not enforced at all as in Problem 2.A. In this
problem as well we observe a significant reduction of the error functional
(6), cf. Figure 6, leading to a good reconstruction of the true flux
$\bar{\phi}$ in Problem 2.A, cf. Figure 7.

Figure 4: [Problem 1:] Dependence of (a) the normalized error functional
${{\mathcal{J}}}({{\phi}}^{(n)})/{{\mathcal{J}}}({{\phi}}_{0})$ and (b)
normalized constraint
$\left[{E}\left(\cdot;{{\phi}}^{(n)}\right)\right]_{T}/\left[{E}\left(\cdot;{{\phi}}_{0}\right)\right]_{T}$
on the iteration $n$ in the solutions of the unconstrained and constrained
optimization problem. For clarity, in 4, the horizontal axis is split into two
regions with different linear scaling separated by the vertical dotted line.
Figure 5: [Problem 1:] (black, solid line) the true flux $\bar{\phi}(t)$, cf.
(8ao), (green, dotted line) the initial guess $\phi_{0}(t)$ in (1) and 8ap,
(blue, dot-dashed line) the optimal flux $\widecheck{\phi}(t)$ in the
unconstrained problem, and (red, dot-solid line) the optimal flux
$\widecheck{\phi}(t)$ in the constrained problem, all as functions of time
$t\in[0,T]$.

Figure 6: [Problem 2:] Dependence of (a) the normalized error functional
${{\mathcal{J}}}({{\phi}}^{(n)})/{{\mathcal{J}}}({{\phi}}_{0})$ and (b)
normalized constraint
$\left[{E}\left(\cdot;{{\phi}}^{(n)}\right)\right]_{T}/\left[{E}\left(\cdot;{{\phi}}_{0}\right)\right]_{T}$
on the iteration $n$ in the solutions of the unconstrained and constrained
optimization problem. For clarity, in 6, the horizontal axis is split into two
regions with different linear scaling separated by the vertical dotted line.
Figure 7: [Problem 2:] (black, solid line) the true flux $\bar{\phi}(t)$, cf.
(8ao), (green, dotted line) the initial guess $\phi_{0}(t)$ in (1) and 8ap,
(blue, dot-dashed line) the optimal flux $\widecheck{\phi}(t)$ in the
unconstrained problem, and (red, dot-solid line) the optimal flux
$\widecheck{\phi}(t)$ in the constrained problem, all as functions of time
$t\in[0,T]$.

###  5.2 3

Problem 3 is similar to Problem 2 in that the retraction operator is not
available. Solutions to Problems 3.A and 3.B are obtained with iterations (1)
where the eddy viscosity representing the Leith model (8ar) is used as the
initial guess. The dependence of the normalized error functional
$\mathcal{J}(\varphi^{(n)})/\mathcal{J}(\varphi_{0})$ and of the normalized
constraint
$\left[\widetilde{\mathcal{E}}(\cdot;\varphi^{(n)})\right]_{T}/\mathcal{E}_{0}$
on the iteration count $n$ is shown, respectively, in Figures 8a and 8b. In
contrast to Problems 1 and 2, Problem 3 does not admit an “exact” solution and
as a result the values of the objective functional (8l) attained at the minima
are not very small, especially in the constrained Problem 3.B. On the other
hand, in analogy with the results for Problem 2, the minimizers found by
solving the constrained Problem 3.B reveal a much slower drift away from the
constraint manifold $\mathcal{M}$ than is evident in the solution of the
unconstrained Problem 3.A.

The optimal eddy viscosities corresponding, via ansatz (8i), to the solutions
of Problems 3.A and 3.B are shown as functions of $s\in\mathcal{I}$ in Figure
9. Although the eddy viscosities obtained in the two cases appear similar,
they reveal some fundamentally different physical properties. Most notably,
the optimal eddy viscosity obtained in Problem 3.B is strictly dissipative. On
the other hand, this quantity obtained in the problem in which the constraint
(8m) is not imposed is negative for small values of the state variable $s$,
meaning that the closure model in fact injects energy in some regions of the
flow. As discussed in detail in [Matharu and Protas, 2022], such behavior was
already observed in optimal closure models obtained with other formulations.

To close this section, we show the enstrophy $\widetilde{\mathcal{E}}(t)$ and
the palinstrophy $\widetilde{\mathcal{P}}(t)$ in the filtered DNS flow and in
the LES flows with the optimal eddy viscosities obtained by solving Problems
3.A and 3.B, respectively, in Figures 10a and 10b. These figures illustrate
how the matching errors in these two quantities are distributed in time and
are presented on a time window twice longer than the “training” window $T$
considered in the error functional (8l) and the constraint (8m). We see that
while for $t\in[0,T]$ the behavior of the two quantities is consistent with
what was observed in Figures 8a and 8b, for $t\in[T,2T]$ both enstrophy and
palinstrophy in the LES flows deviate more significantly from their values in
the filtered DNS flow.

Figure 8: [Problem 3:] Dependence of (a) the normalized error functional
${\mathcal{J}}({{\varphi}}^{(n)})/{\mathcal{J}}({{\varphi}}_{0})$ and (b) the
normalized constraint
$\left[\widetilde{\mathcal{E}}(\cdot;\varphi^{(n)})\right]_{T}/\mathcal{E}_{0}$
on the iteration $n$ for the solutions of the unconstrained and constrained
optimization problem.

Figure 9: [Problem 3:] Dependence of the optimal eddy viscosity
$\widecheck{\nu}$ on $\sqrt{s}$ for (blue, dashed-dot line)
${\varphi}\in\mathcal{S}$ found by solving Problem 3.A and (red, dashed line)
${\varphi}\in\mathcal{M}$ found by solving Problem 3.B. The initial guess in
(1) is given by the Leith model $\nu_{L}(s)$, cf. 8ar, and is shown as green
dotted line

.

Figure 10: [Problem 3:] Dependence of (a) the enstrophy 8ka and (b) the
palinstrophy 8kb on time for (black, solid line) the filtered DNS flow and the
LES flows with the optimal eddy viscosities corresponding to (blue, dashed-dot
line) $\widecheck{{\varphi}}\in\mathcal{S}$ found by solving Problem 3.A and
(red, dashed line) $\widecheck{{\varphi}}\in\mathcal{M}$ found by solving
Problem 3.B. Thick and thin lines represent to, respectively, time in the
“training window” ($t\in[0,T]$) and in the extended window ($t\in(T,2T]$).

##  6 Discussion & Conclusions

Adjoint-based methods have become a workhorse in the solution of unconstrained
PDE optimization problems [Gunzburger, 2003]. They make it possible to
conveniently determine the gradient (sensitivity) of the objective functional
with respect to a distributed control variable, which can then be used in
various gradient descent algorithms [Nocedal and Wright, 2006]. Unlike most
constraints imposed on the control variable, constraints on the state
variables are generally harder to satisfy since they define, via solutions of
the governing system, complicated manifolds in the space of control variables.
In the present study we demonstrate how the adjoint-based framework can be
extended to handle such constraints approximately, which is done by
constructing a projection of the gradient of the objective functional onto a
subspace tangent to the constraint manifold at the given iteration. This
projection is realized by solving another adjoint problem which is defined in
terms of the same adjoint operator as the adjoint system employed to determine
the gradient, but with different forcing. Thus, for a state constraint
defining a codimension-$m$ manifold, a single iteration of the gradient
algorithm (1) requires a total of $(m+1)$ adjoint solves to evaluate the
gradient and the projection on the tangent subspace. We focus on the
“optimize-then-discretize” paradigm [Gunzburger, 2003] in the infinite-
dimensional setting, where it is more difficult to use “black-box”
optimization software that could otherwise facilitate imposition of such
constraints. Regularity of both the gradient and the projection is carefully
considered, which is done by leveraging the Riesz theorem.

Our approach is illustrated with two test problems, both involving a nonlinear
scalar constraint: a simple problem describing heat conduction in 1D which is
used to introduce the method in a clean setting and a more complicated problem
concerning the optimal design of turbulence closures in 2D incompressible
flows governed by the Navier-Stokes system. This latter problem has a
nonstandard structure [Matharu and Protas, 2022]. For each of the optimization
problems we consider the constrained and unconstrained version with and
without the nonlinear constraint. In addition, for comparison, we also
consider a version of the first problem where the constraint is homogeneous
and hence can be enforced exactly with a retraction operator.

Our computation of the gradients and of the normal elements defining the
tangent subspaces is carefully validated in Section 4. The computational
results presented in Section 5 show that even though in the solutions of
Problems 2.B and 3.B the constraints are not satisfied exactly, the obtained
optimizers lie much closer to the constraint manifold $\mathcal{M}$ than in
Problems 2.A and 3.A where the constraint is not imposed. However, as
expected, the optimizers in the constrained problems correspond to larger
values of the objective functionals. Importantly, enforcement of the
constraint via projection onto the subspace tangent to the constraint manifold
is performed in just as straightforward manner as the computation of the
gradient. Therefore, the proposed approach is likely to benefit practical
applications involving PDE optimization problems.

## Acknowledgements

Partial support for this research was provided through an NSERC (Canada)
Discovery Grant. The first author acknowledges support from the Verg
foundation.

## References

  * Luenberger [1969] D. Luenberger.  _Optimization by Vector Space Methods_.  John Wiley and Sons, 1969.
  * Tröltzsch [2010] F. Tröltzsch.  _Optimal Control of Partial Differential Equations: Theory, Methods and Applications_ , volume 112 of _Graduate Studies in Mathematics_.  American Mathematical Society, 2010.
  * Lions [1968] J.L. Lions.  _Contrôle Optimal de Systèmes Gouvernés par des Équations aux Dérivées Partielles_.  Dunod, Paris, 1968.  (English translation, Springer-Verlag, New-York (1971)).
  * Gunzburger [2003] M. D Gunzburger.  _Perspectives in Flow Control and Optimization_.  SIAM, 2003.  doi: 10.1137/1.9780898718720.  URL https://epubs.siam.org/doi/abs/10.1137/1.9780898718720.
  * Bukshtynov and Protas [2013] V. Bukshtynov and B. Protas.  Optimal reconstruction of material properties in complex multiphysics phenomena.  _J. Comput. Phys._ , 242:889–914, 2013.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2013.02.034.  URL https://doi.org/10.1016/j.jcp.2013.02.034.
  * Vergnault and Sagaut [2014] E. Vergnault and P. Sagaut.  An adjoint-based lattice boltzmann method for noise control problems.  _J. Comp. Phys._ , 276:39–61, 2014.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2014.07.027.  URL https://www.sciencedirect.com/science/article/pii/S0021999114005208.
  * Costanzo et al. [2022] S. Costanzo, T. Sayadi, M. Fosas de Pando, P.J. Schmid, and P. Frey.  Parallel-in-time adjoint-based optimization – application to unsteady incompressible flows.  _J. Comp. Phys._ , 471:111664, 2022.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2022.111664.  URL https://www.sciencedirect.com/science/article/pii/S0021999122007276.
  * Sirignano et al. [2023] J. Sirignano, J. MacArt, and K. Spiliopoulos.  Pde-constrained models with neural network terms: Optimization and global convergence.  _J. Comp. Phys._ , 481:112016, 2023.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2023.112016.  URL https://www.sciencedirect.com/science/article/pii/S0021999123001110.
  * Absil et al. [2008] P.-A. Absil, R. Mahony, and R. Sepulchre.  _Optimization Algorithms on Matrix Manifolds_.  Princeton University Press, 2008.
  * Matharu and Protas [2022] P. Matharu and B. Protas.  Optimal eddy viscosity in closure models for two-dimensional turbulent flows.  _Phys. Rev. Fluids_ , 7:044605, 2022.  doi: 10.1103/PhysRevFluids.7.044605.  URL https://link.aps.org/doi/10.1103/PhysRevFluids.7.044605.
  * Adams and Fournier [2003] R. A. Adams and J. F. Fournier.  _Sobolev Spaces_.  Elsevier/Academic Press, Amsterdam, second edition, 2003.  doi: 10.1016/S0079-8169(03)X8001-0.
  * Lesieur [1993] M. Lesieur.  _Turbulence in Fluids_.  Kluwer Academic Publishers, Dordrecht, Boston, London, 2nd edition, 1993\.
  * Pope [2000] S. B. Pope.  _Turbulent flows_.  Cambridge University Press, Cambridge, 2000.  ISBN 0-521-59886-9.  doi: 10.1017/CBO9780511840531.  URL https://doi.org/10.1017/CBO9780511840531.
  * Davidson [2015] P. A. Davidson.  _Turbulence: An introduction for scientists and engineers_.  Oxford University Press, Oxford, second edition, 2015.  ISBN 978-0-19-872259-5; 978-0-19-872258-8.  doi: 10.1093/acprof:oso/9780198722588.001.0001.  URL https://doi.org/10.1093/acprof:oso/9780198722588.001.0001.
  * Kutz [2017] J. N. Kutz.  Deep learning in fluid dynamics.  _J. Fluid Mech._ , 814:1–4, 2017.  doi: 10.1017/jfm.2016.803.
  * Gamahara and Hattori [2017] Masataka Gamahara and Yuji Hattori.  Searching for turbulence models by artificial neural network.  _Phys. Rev. Fluids_ , 2:054604, 2017.  doi: 10.1103/PhysRevFluids.2.054604.  URL https://link.aps.org/doi/10.1103/PhysRevFluids.2.054604.
  * Jimenez [2018] J. Jimenez.  Machine-aided turbulence theory.  _J. Fluid Mech._ , 854:R1, 2018.  doi: 10.1017/jfm.2018.660.
  * Duraisamy et al. [2019] K. Duraisamy, G. Iaccarino, and H. Xiao.  Turbulence Modeling in the Age of Data.  _Annu. Rev. Fluid Mech._ , 51(1):357–377, 2019\.  doi: 10.1146/annurev-fluid-010518-040547.  URL https://doi.org/10.1146/annurev-fluid-010518-040547.
  * Duraisamy [2021] Karthik Duraisamy.  Perspectives on machine learning-augmented reynolds-averaged and large eddy simulation models of turbulence.  _Phys. Rev. Fluids_ , 6:050504, 2021.  doi: 10.1103/PhysRevFluids.6.050504.  URL https://link.aps.org/doi/10.1103/PhysRevFluids.6.050504.
  * Pawar and San [2021] Suraj Pawar and Omer San.  Data assimilation empowered neural network parametrizations for subgrid processes in geophysical flows.  _Phys. Rev. Fluids_ , 6:050501, 2021.  doi: 10.1103/PhysRevFluids.6.050501.  URL https://link.aps.org/doi/10.1103/PhysRevFluids.6.050501.
  * Waschkowski et al. [2022] F. Waschkowski, Y. Zhao, R. Sandberg, and J. Klewicki.  Multi-objective cfd-driven development of coupled turbulence closure models.  _J. Comp. Phys._ , 452:110922, 2022.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2021.110922.  URL https://www.sciencedirect.com/science/article/pii/S0021999121008172.
  * Matharu and Protas [2020] P. Matharu and B. Protas.  Optimal closures in a simple model for turbulent flows.  _SIAM J. Sci. Comput._ , 42(1):B250–B272, 2020\.  doi: 10.1137/19M1251941.  URL https://doi.org/10.1137/19M1251941.
  * Protas et al. [2004] B. Protas, T. R. Bewley, and G. Hagen.  A computational framework for the regularization of adjoint analysis in multiscale PDE systems.  _J. Comp. Phys._ , 195(1):49 – 89, 2004.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2003.08.031.  URL http://www.sciencedirect.com/science/article/pii/S0021999103005114.
  * Leith [1968] C. E. Leith.  Diffusion approximation for two-dimensional turbulence.  _The Physics of Fluids_ , 11(3):671–672, 1968\.  doi: 10.1063/1.1691968.  URL https://aip.scitation.org/doi/abs/10.1063/1.1691968.
  * Leith [1971] C. E. Leith.  Atmospheric predictability and two-dimensional turbulence.  _Journal of the Atmospheric Sciences_ , 28(2):145–161, 1971.  doi: 10.1175/1520-0469(1971)028<0145:APATDT>2.0.CO;2.  URL https://doi.org/10.1175/1520-0469(1971)028<0145:APATDT>2.0.CO;2.
  * Leith [1996] C.E. Leith.  Stochastic models of chaotic systems.  _Physica D: Nonlinear Phenomena_ , 98(2):481 – 491, 1996.  ISSN 0167-2789.  doi: https://doi.org/10.1016/0167-2789(96)00107-8.  URL http://www.sciencedirect.com/science/article/pii/0167278996001078.  Nonlinear Phenomena in Ocean Dynamics.
  * Berger [1977] M. S. Berger.  _Nonlinearity and Functional Analysis_.  Academic Press, 1977.
  * [28] P. Matharu.  Constraint Optimization - Heat Equation (HeatOptConstr) (2023).  URL https://github.com/pipmath/HeatOptConstr.git.
  * Alimo et al. [2021] Ryan Alimo, Daniele Cavaglieri, Pooriya Beyhaghi, and Thomas R. Bewley.  Design of IMEXRK time integration schemes via Delaunay-based derivative-free optimization with nonconvex constraints and grid-based acceleration.  _J. Glob. Optim._ , 79(3):567–591, 2021.  ISSN 1573-2916.  doi: 10.1007/s10898-019-00855-1.  URL https://doi.org/10.1007/s10898-019-00855-1.
  * Trefethen [2013] L. N. Trefethen.  _Approximation Theory and Approximation Practice_.  SIAM, Philadelphia, 2013.  ISBN 978-1-611972-39-9.  URL https://people.maths.ox.ac.uk/trefethen/ATAP/.
  * Bracco and McWilliams [2010] A. Bracco and J. C. McWilliams.  Reynolds-number dependency in homogeneous, stationary two-dimensional turbulence.  _Journal of Fluid Mechanics_ , 646:517–526, 2010.  doi: 10.1017/S0022112009993661.
  * Press et al. [2007] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery.  _Numerical Recipes 3rd Edition: The Art of Scientific Computations_.  Cambridge University Press, 2007.  ISBN 978-0-521-88068-8.
  * Nocedal and Wright [2006] Jorge Nocedal and Stephen Wright.  _Numerical Optimization_.  Springer Series in Operations Research and Financial Engineering. Springer, second edition, 2006.  ISBN 978-0387-30303-1.  doi: 10.1007/978-0-387-40065-5.
  * Engl et al. [1996] H. Engl, M. Hanke, and A. Neubauer.  _Regularization of Inverse Problems_.  Kluver, Dordrecht, 1996.