Adjoint fluid

In neko, and by extension neko-top, we solve the Navier–Stokes equations

\[ {\frac {\partial \mathbf {u} }{\partial t}} + (\mathbf {u} \cdot \nabla)\mathbf {u} = -\nabla p +{\frac {1}{Re}}\nabla ^{2}\mathbf {u} + \mathbf{f}, \text{ in } \Omega,\\ \nabla \cdot \mathbf {u} = 0, \text{ in } \Omega, \]

where \(\mathbf {u}(\mathbf{x},t)\) denotes the velocity field, \(p(\mathbf{x},t)\) the pressure field, \(\mathbf{f}\) a forcing term and where \(Re\) denotes the Reynolds number. The governing equations are subjected boundary conditions on

\[ \frac{1}{Re} \nabla \mathbf{u} \cdot \mathbf{n} -p \mathbf{n} = 0, \text{ on } \Gamma_O, \\ \mathbf{u} = \mathbf{u}_\text{in}, \text{ on } \Gamma_D, \\ \]

where \(\Gamma_N\) and \(\Gamma_D\) denote outflow and dirichlet boundaries respectively.

A common formulation of the adjoint Navier-Stokes equations reads

\[ {-\frac {\partial \mathbf {u}^\dagger }{\partial t}} + (\nabla \mathbf {u})^T \mathbf {u}^\dagger - (\mathbf {u} \cdot \nabla) \mathbf {u}^\dagger = -\nabla p^\dagger +{\frac {1}{Re}}\nabla ^{2}\mathbf {u}^\dagger + \mathbf{f}^\dagger, \\ \nabla \cdot \mathbf {u} ^\dagger= 0, \]

where where \(\mathbf {u}^\dagger(\mathbf{x},t)\) denotes the adjoint velocity field, \(p^\dagger(\mathbf{x},t)\) the adjoint pressure field and \(\mathbf{f}^\dagger\) denoting a forcing term applied to the adjoint system which generally arises as a consequence of objective functions being evaluated.

The corresponding boundary conditions generally read

\[ \frac{1}{Re} \nabla \mathbf{u}^\dagger \cdot \mathbf{n} -p^\dagger \mathbf{n} + (\mathbf{u} \cdot \mathbf{n}) \mathbf{u}^\dagger = \mathbf{0} \text{ on } \Gamma_N, \\ \mathbf{u} = \mathbf{0}, \text{ on } \Gamma_D. \\ \]

The spectral element method which underpins neko solves the above system of equations using the weak formulation, which has important implications when solving the adjoint system. Primarily, when deriving the adjoint system in strong form, one introduces additional boundary terms on Neumann boundaries, and more importantly, one applies the divergence free condition strongly.

An alternative to the above adjoint system is to remain in weak form and not integrate the convective term by parts, resulting in no additional boundary terms, and more importantly, no pointwise application of the divergence free condition.

This alternative system is what is implemented in neko-top, which reads (in weak form),

\[ -\int_\Omega \mathbf{v}\cdot {\frac {\partial \mathbf {u}^\dagger } {\partial t}} + \int_\Omega \mathbf{v}\cdot (\nabla \mathbf {u})^T \mathbf {u}^\dagger + \int_\Omega \nabla \mathbf{v}\cdot (\mathbf {u} \otimes \mathbf {u}^\dagger ) = -\int_\Omega \mathbf{v}\cdot \nabla p^\dagger +{\frac {1}{Re}}\int_\Omega \nabla \mathbf{v}\cdot \nabla \mathbf {u}^\dagger + \int_\Omega \mathbf{v}\cdot \mathbf{f}^\dagger, \\ \int_\Omega q \nabla \cdot \mathbf {u} ^\dagger= 0, \]

where \(\mathbf {v}\) is a test function.

Adjoint scalar

The adjoint scalar equation takes the following form

\[ {\frac {\partial \phi }{\partial t}} + (\mathbf {u} \cdot \nabla)\phi = {\frac {1}{Pe}}\nabla ^{2}\phi , \text{ in } \Omega,\\ \]

where \(\phi(\mathbf{x},t)\) denotes the scalar field, and where \(Pe\) denotes the Peclet number. In the context of conjugate heat transfer, the velocity equation is often coupled to the scalar equation through the Boussinesq approximation for instance, however, in lieu of this coupling the scalar is often referred to as a "passive scalar" to imply the one way coupling.

In neko-top the adjoint scalar reads (in weak form),

\[ -\int_\Omega \psi\cdot {\frac {\partial \phi^\dagger } {\partial t}} + \int_\Omega (\nabla \psi \cdot \mathbf {u}) \phi^\dagger = {\frac {1}{Pe}}\int_\Omega \nabla \psi \cdot \nabla \phi^\dagger,\\ \]

where \(\phi(\mathbf{x},t)^\dagger\) denotes the adjoint scalar field and where \(\psi\) is a test function. In addition, the perturbation of the term \((\mathbf {u} \cdot \nabla)\phi\) results in an additional term in adjoint velocity momentum equation, which now reads

\[ -\int_\Omega \mathbf{v}\cdot {\frac {\partial \mathbf {u}^\dagger } {\partial t}} + \int_\Omega \mathbf{v}\cdot (\nabla \mathbf {u})^T \mathbf {u}^\dagger + \int_\Omega \nabla \mathbf{v}\cdot (\mathbf {u} \otimes \mathbf {u}^\dagger ) \underline{+ \int_\Omega (\mathbf{v}\cdot \nabla \phi ) \phi^\dagger} = -\int_\Omega \mathbf{v}\cdot \nabla p^\dagger +{\frac {1}{Re}}\int_\Omega \nabla \mathbf{v}\cdot \nabla \mathbf {u}^\dagger + \int_\Omega \mathbf{v}\cdot \mathbf{f}^\dagger. \]

It can be seen from the above equations that due to the one way coupling between \(\mathbf{u}\) and \(\phi\), when solving the forward problem one must first solve for \(\mathbf{u}\) and then solve for \(\phi\). However, in the adjoint the opposite occurs as the equation for \(\mathbf{u}^\dagger\) depends on \(\phi^\dagger\), hence in neko-top we first solve for the adjoint scalar and then solve for the adjoint velocity.

Multiple scalars

In the case of using multiple scalar transport equations involving multiple scalar fields \(\phi_i(\mathbf{x},t)\), the assumption is made that all scalars are passive, that is, coupled in one way to the velocity, and, in no way coupled to one another. If this is not the case, additional care must be taken by the user when deriving their adjoint equations, and will likely introduce additional coupling terms in the adjoint equations. In lieu of this, each adjoint equation will take an identical form to that stated above, and the adjoint velocity momentum equation will contain the summation over all the additional coupling terms

\[ - \int_\Omega \mathbf{v} \cdot \frac{\partial \mathbf{u}^\dagger}{\partial t} + \int_\Omega \mathbf{v} \cdot (\nabla \mathbf{u})^T \mathbf{u}^\dagger + \int_\Omega \nabla \mathbf{v} \cdot (\mathbf{u} \otimes \mathbf{u}^\dagger) + \underline{\sum_i \int_\Omega (\mathbf{v} \cdot \nabla \phi_i) \phi_i^\dagger} = - \int_\Omega \mathbf{v} \cdot \nabla p^\dagger + \frac{1}{Re} \int_\Omega \nabla \mathbf{v} \cdot \nabla \mathbf{u}^\dagger + \int_\Omega \mathbf{v} \cdot \mathbf{f}^\dagger. \]

Note: If a given \(\phi_i(\mathbf{x},t)\) is not contained in the objective, there will be no resulting forcing to drive the adjoint scalar equation, effectively yielding \(\phi_i(\mathbf{x},t)^\dagger = \mathbf{0}\), and hence, will in theory have no contribution to the adjoint momentum equation.

Immersed Boundary Methods

Following a Brinkman style immersed boundary method, the presence of an immersed object is imposed by the Brinkman forcing term \(\mathbf{f} = - \chi \mathbf{u}\), where \(\chi\) is the spatially dependent Brinkman coefficient, satisfying

\[ \chi = \begin{cases} 0 & \text{in the fluid region,} \\ \overline{\chi} & \text{in the solid region.} \end{cases} \]

Considering \(\overline{\chi}\) to be a large value, this discontinuous forcing term models a momentum loss in solid region, and thereby simulating porous media with very low permeability. It is worth noting that while the Brinkman penalization method is rooted in the idea of modelling solid regions as porous media with vanishing permeability, the approach proposed by Goldstein frames the interaction of the solid on the fluid as a control problem, where the feedback force is tuned to drive the velocity to zero in the solid region. Regardless of their different motivation, both methods result in a similar mathematical structure. More information regarding the mapping of \(\chi\) can be found in Mapping cascade.

Note: In the future we will provide a full adjoint derivation in this section of the theory guide. This will tie together all aspects from the mapping, to how these terms arise etc. For now we are simply documenting the equations being solved.

Time integration

This section summarizes the discrete adjoint for the pressure-correction time integration used in Neko-TOP. We start from a right-handed Riemann sum for the objective,

\[ F \approx \Delta t \sum_{n=1}^{N} f(\mathbf{u}_n), \qquad \mathrm{d}F = \Delta t \sum_{n=1}^{N} \int_\Omega \nabla_{\mathbf{u}_n} f(\mathbf{u}_n)\cdot \delta \mathbf{u}_n \, \mathrm{d}\Omega. \]

The velocity-pressure splitting currently implemented Neko follows that of Karniadakis et al. 1991 and is commonly referred to as the \( \mathbb{P}_N-\mathbb{P}_N\) formulation, whereby we introduce an intermediate velocity \(\mathbf{u}_n^{\star}\) and write the forward time step as three sub-problems:

\[ \int_\Omega \mathbf{v}_n \cdot \frac{\mathbf{u}_n^{\star}}{\Delta t} \, \mathrm{d}\Omega = \int_\Omega \mathbf{v}_n \cdot \frac{\mathbf{u}_{n-1}}{\Delta t} \, \mathrm{d}\Omega - \int_\Omega \mathbf{v}_n \cdot \mathcal{C}_{\mathbf{u}_{n-1}}\,\mathbf{u}_{n-1} \, \mathrm{d}\Omega - \int_\Omega \mathbf{v}_n \cdot \chi\,\mathbf{u}_{n-1} \, \mathrm{d}\Omega, \]

\[ \int_\Omega \nabla q_n \cdot \nabla p_n \, \mathrm{d}\Omega = \int_\Omega \nabla q_n \cdot \frac{\mathbf{u}_n^{\star}}{\Delta t} \, \mathrm{d}\Omega - \int_\Omega \nabla q_n \cdot \nabla \times (\nabla \times \mathbf{u}_{n-1}) \, \mathrm{d}\Omega, \]

\[ \int_\Omega \mathbf{w}_n \cdot \frac{\mathbf{u}_n}{\Delta t} \, \mathrm{d}\Omega + \int_\Omega \mathbf{w}_n \cdot \mathcal{L}\mathbf{u}_n \, \mathrm{d}\Omega = \int_\Omega \mathbf{w}_n \cdot \frac{\mathbf{u}_n^{\star}}{\Delta t} \, \mathrm{d}\Omega - \int_\Omega \mathbf{w}_n \cdot \nabla p_n \, \mathrm{d}\Omega. \]

Here \(\mathcal{C}_{\mathbf{u}_{n-1}}\) denotes the linearized convection operator about \(\mathbf{u}_{n-1}\), \(\chi\) denotes any explicit linear operator (e.g. Brinkman), and \(\mathcal{L}\) is the implicit diffusion/Helmholtz operator.

Perturbing these equations and collecting terms yields a discrete adjoint that is solved in three steps (backward in time). First, solve for the adjoint velocity associated with the implicit step,

\[ -\int_\Omega \frac{\mathbf{u}_n^{*\dagger}}{\Delta t} \, \mathrm{d}\Omega + \int_\Omega \mathcal{L}^\dagger \mathbf{u}_n^{*\dagger} \, \mathrm{d}\Omega = \int_\Omega \nabla_{\mathbf{u}_n} f(\mathbf{u}_n)\,\Delta t \, \mathrm{d}\Omega + \int_\Omega \frac{\mathbf{u}_{n+1}^{\dagger}}{\Delta t} \, \mathrm{d}\Omega - \int_\Omega \mathcal{C}^\dagger_{\mathbf{u}_n}\,\mathbf{u}_{n+1}^{\dagger} \, \mathrm{d}\Omega - \int_\Omega \chi^\dagger\, \mathbf{u}_{n+1}^{\dagger} \, \mathrm{d}\Omega. \]

Second, solve the adjoint pressure equation,

\[ \int_\Omega \nabla q_n \cdot \nabla p_n^{\dagger} \, \mathrm{d}\Omega = \int_\Omega -\mathbf{u}_n^{*\dagger} \cdot \nabla p_n^{\dagger} \, \mathrm{d}\Omega. \]

Finally, recover the adjoint velocity associated with the explicit step,

\[ \mathbf{u}_n^{\dagger} = \mathbf{u}_n^{*\dagger} - \nabla p_n^{\dagger}. \]

With these adjoint variables, the discrete gradient contribution associated with an explicit linear operator \(\chi\) is

\[ \mathrm{d}F(\chi,\delta \chi) = -\sum_{n=1}^{N} \int_\Omega \delta \chi\, \mathbf{u}_{n+1}^{\dagger} \mathbf{u}_n \, \mathrm{d}\Omega. \]

The order above matches the implementation: solve the implicit adjoint step, then the adjoint pressure correction, then form the projection to obtain the adjoint velocity used in the sensitivity accumulation.

The curl–curl term

The forward pressure equation above contains the curl–curl term \(-\int_\Omega \nabla q_n \cdot \nabla \times (\nabla \times \mathbf{u}_{n-1})\,\mathrm{d}\Omega\), which Neko evaluates on the extrapolated velocity \(\mathbf{u}_e\) (the same extrapolation as for the explicit terms). Its adjoint therefore enters the adjoint velocity equation as an explicit load in the adjoint pressure. In the continuous setting, \(\nabla \times \nabla p^\dagger = 0\) reduces this load to a boundary integral. The discrete operators do not satisfy that identity, so a boundary-only form is not the transpose of what the forward solver computes and biases the gradient. The implementation (`adjoint_curl_curl`) instead uses the exact transpose of the discrete term, which acts on the whole volume.

Before gather-scatter, the forward solver adds to the pressure residual

\[ r_{cc}(\mathbf{u}_e) = -\frac{\mu}{\rho}\, G^T M \mathbf{W}, \qquad \mathbf{W} = M D_c M D_c\, \mathbf{u}_e, \]

where \(D_c\) is the element-local, pointwise curl, \(M\mathbf{y} = \bar{B}^{-1}\,\mathrm{gs}(B\mathbf{y})\) is the mass-weighted average applied by Neko's curl (with \(B\) the element-local diagonal mass matrix, \(\mathrm{gs}\) the gather-scatter sum and \(\bar{B} = \mathrm{gs}(B)\) the assembled mass matrix), and \(G\) is the weak gradient, \(G_i p = B D_i p\). Because \(\mathrm{gs}\) replaces every degree of freedom in a shared group by the sum over that group, \(M\mathbf{y}\) already carries the same value on every copy of a shared degree of freedom; a second application of \(M\) multiplies that group-uniform value by the local \(B\), sums it back to \(\bar{B}\) times the value, then divides by \(\bar{B}\) again, so it returns the input unchanged. \(M\) is therefore a projection, \(M^2 = M\), and the \(M^2\) that would otherwise appear where the leading \(M\) of \(r_{cc}\) meets the leading \(M\) of \(\mathbf{W}\) collapses to a single \(M\). Since \(r_{cc}\) is linear in \(\mathbf{u}_e\), the adjoint load is

\[ \mathbf{L} = -\Big(\frac{\partial r_{cc}}{\partial \mathbf{u}_e}\Big)^T p^\dagger = \frac{\mu}{\rho}\, D_c^T M^T D_c^T M^T G\, p^\dagger, \qquad M^T\mathbf{y} = B\,\mathrm{gs}(\bar{B}^{-1}\mathbf{y}), \]

carrying two \(M^T\) factors rather than the three a term-by-term transpose of \(\mathbf{W}\) would suggest. It is added to the adjoint forcing before the explicit extrapolation. That step scales the forcing by \(\rho\) and weights the lagged values, so the adjoint velocity right-hand side receives \(\mu\, D_c^T M^T D_c^T M^T G\, p^\dagger\), weighted over the adjoint pressures of the later forward steps exactly as the forward extrapolation weights \(\mathbf{u}_e\).

Warning: A mesh that is not three-dimensional is rejected before the adjoint fluid is initialised: the curl–curl transpose is exact in three dimensions only. The forward solver also applies the curl–curl term through the symmetry-surface term of the pressure residual and through the rotations at cyclic boundaries; neither transpose is implemented, so a symmetry velocity boundary condition in either case.fluid.boundary_conditions or case.adjoint_fluid.boundary_conditions, and case.fluid.cyclic set to true, are rejected when the adjoint is set up.