Adjoint of Physics Solution

The fundamental task in solving the adjoint of the modeled physical phenomena is to compute the cost function sensitivity with respect to the mesh, $\frac{d L^{T}}{d X}$ .

The cost function sensitivity with respect to the mesh is the product of the two rightmost terms in Eqn. (5067):

Figure 1. EQUATION_DISPLAY

\frac{d L^{T}}{d X} = \frac{\partial L^{T}}{\partial X} + \frac{\partial Q^{T}}{\partial X} \frac{\partial L^{T}}{\partial Q}

(5075)

where the calculation of this sensitivity takes place in two steps:

Iteratively solving the adjoint solution, providing the sensitivity of the report objective with respect to the residuals of the governing equations, $\frac{d L^{T}}{d R}$ .
Computing the product of the residuals with respect to the mesh, $\frac{\partial R^{T}}{\partial X}$ , with the adjoint solution and adding it to the explicit mesh dependency of the objective, $\frac{\partial L^{T}}{\partial X}$ . This second step is performed as a one-off operation after the adjoint solution is converged.

These steps require the differentiation of the primal solution process, including the solution of the flow equations, energy equations and the turbulence equations when applicable. The derivation of the coupled flow model is described below. The derivation of the coupled solid energy model follows a similar pattern.

Primal Solution

The aim of the iterative solvers within Simcenter STAR-CCM+ is to advance a set of variables so that they may satisfy some set of discretized governing equations. For example, the governing equations for the Coupled Flow Model are the Navier -Stokes equations Eqn. (5076).

The preconditioned Navier-Stokes equations are described as:

Figure 2. EQUATION_DISPLAY

Γ \frac{\partial}{\partial t} \int_{V} Q d V + \oint [F - G] \cdot d A = \int_{V} H d V

(5076)

where $Γ$ is the preconditioning matrix, $F$ and $G$ are the inviscid and viscous fluxes, respectively, and $H$ are the body forces.

For generality, suppose there are $N$ different iterative solvers, each concerned with advancing a sub set of variables such that they satisfy some governing equations. In the case of a turbulent flow problem, these solvers would be the Coupled Implicit Solver, a Turbulence Solver and a Turbulent Viscosity Solver for example. For the $i^{t h}$ solver, the set of variables is referred to as $Q_{i}$ and the discretized governing equations can be represented in residual form as $R_{i}$ . Because the solvers are coupled, this residual can be a function of all variables in the simulation. Using this notation, the primal solution boils down to solving the following set of coupled non-linear equations:

Figure 3. EQUATION_DISPLAY

\begin{matrix} R_{1} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \\ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} R_{2} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \\ ⋮ \end{matrix} \\ R_{i} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \end{matrix} \\ ⋮ \end{matrix} \\ R_{N} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \end{matrix} \end{matrix}

(5077)

For a steady problem, this process takes on the form of a non-linear block Gauss-Seidel method with each solver building its residual based on the latest variable values and advancing its set of variables based on a local preconditioner $[P_{i}]$ . For solver $i$ , the update for non-linear iteration $k$ is given as:

Figure 4. EQUATION_DISPLAY

\begin{matrix} R_{i}^{k} = R_{i} (Q_{1}^{k + 1}, Q_{2}^{k + 1} \dots, Q_{i}^{k}, \dots, Q_{N}^{k}) \\ Q_{i}^{k + 1} = Q_{i}^{k} - {[P_{i}]}^{- 1} R_{i}^{k} \end{matrix}

(5078)

In the case of the Coupled Flow Model, the preconditioner is an approximate linearization of the Navier-Stokes equations w.r.t the flow variables (P, U, V, W, T), defined below:

Figure 5. EQUATION_DISPLAY

[P] = \frac{V}{δ t} Γ + \sum_{j}^{N_{f}} (\frac{\partial F_{j}}{\partial Q} - \frac{\partial G_{j}}{\partial Q}) A_{j}

(5079)

Adjoint Solution

To derive the steady adjoint solution, the coupled non-linear equations defined in the primal solution must be differentiated and transposed. To determine the derivative w.r.t $X$ , the chain rule is applied to each primal residual equation:

Figure 6. EQUATION_DISPLAY

\begin{matrix} \frac{d R_{1}}{d X} + \frac{\partial R_{1}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{1}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{1}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{1}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \\ \frac{d R_{2}}{d X} + \frac{\partial R_{2}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{2}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{2}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{2}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \\ \begin{matrix} \begin{matrix} \begin{matrix} ⋮ \\ \frac{d R_{i}}{d X} + \frac{\partial R_{i}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{i}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{i}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{i}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \end{matrix} \\ ⋮ \end{matrix} \\ \frac{d R_{N}}{d X} + \frac{\partial R_{N}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{N}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{N}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{N}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \end{matrix} \end{matrix}

(5080)

The result of the above differentiation is a large, fully coupled linear system between all the solvers:

Figure 7. EQUATION_DISPLAY

[\frac{\partial R}{\partial Q}] \frac{\partial Q}{\partial X} = - \frac{\partial R}{\partial X}

(5081)

Unlike the preconditioner matrix used to advance the primal, the linearizations in the above process are exact and account for 2nd order terms.

Figure 8. EQUATION_DISPLAY

\frac{\partial Q^{T}}{\partial X} = - \frac{\partial R^{T}}{\partial X} {[\frac{\partial R}{\partial Q}]}^{- T}

(5082)

As the inverse calculation for ${[\frac{\partial R}{\partial Q}]}^{- T}$ is computationally expensive, Eqn. (5081) is multiplied with the previously calculated result of ${\frac{\partial L}{\partial Q}}^{T}$ (see Solution Analysis):

Figure 9. EQUATION_DISPLAY

{\frac{d L}{d X}}^{T} = {\frac{\partial L}{\partial X}}^{T} - {\frac{\partial R}{\partial X}}^{T} {[\frac{\partial R}{\partial Q}]}^{- T} {\frac{\partial L}{\partial Q}}^{T}

(5083)

where ${\frac{\partial R}{\partial X}}^{T}$ is known because $R$ is known. The last two terms define the adjoint solution variable, $Λ_{Q}$ :

Figure 10. EQUATION_DISPLAY

{Λ_{Q} = [\frac{\partial R}{\partial Q}]}^{- T} {\frac{\partial L}{\partial Q}}^{T}

(5084)

which can be interpreted as ${\frac{d L}{d R}}^{T}$ .

To avoid calculating and assembling a large, complex matrix, the adjoint solution process utilizes algorithms which only require matrix vector products, namely Restarted GMRES and Defect correction. The defection-correction can be written as:

Figure 11. EQUATION_DISPLAY

{[P]}^{T} {δ Λ}_{Q} = \frac{\partial L}{\partial Q} - [\frac{\partial R}{\partial Q}] Λ_{Q} = - R_{Λ_{Q}} (Λ_{b f Q})

(5085)

For GMRES, the system is solved directly and Eqn. (5085) is used as a preconditioner within the GMRES algorithm. To implement a sweep of defect correction across multiple solvers, the primal solution process is applied in reverse order and each adjoint solver operates on a column of the overall matrix. At each iteration, the solver computes an adjoint residual and advances its adjoint variables using the transpose of the local preconditioner. The only departure from the primal structure is that the adjoint is advanced in a block Jacobi-like manner, where all solvers use the adjoint solution from the previous iteration when computing their adjoint residuals. Advancing the solution in this way saves memory when solving for multiple adjoint solutions.

Figure 12. EQUATION_DISPLAY

\begin{matrix} R_{Λ_{Q, i}}^{k} = \frac{\partial R_{1}}{\partial Q_{i}} Λ_{Q, 1}^{k} + \frac{\partial R_{2}}{\partial Q_{i}} Λ_{Q, 2}^{k} + \dots + \frac{\partial R_{i}}{\partial Q_{i}} Λ_{Q, i}^{k} + \dots + \frac{\partial R_{N}}{\partial Q_{i}} Λ_{Q, N}^{k} \\ Λ_{Q, i}^{k + 1} = Λ_{Q, i}^{k} - {[P_{i}]}^{- T} R_{Λ_{Q, i}}^{k} \end{matrix}

(5086)

The above process represents the fully coupled adjoint solution for the primal process. For example, if the primal consists of a flow solution using the SA turbulence model, the above represents the implementation when the Adjoint SA model is selected. It is also possible to introduce a frozen treatment for a particular solver. This is the treatment applied to turbulence solvers when the Adjoint Frozen Turbulence Model is selected. When a solver is frozen, the row representing that solvers residual and the column representing that solvers variables are removed from the linear system of Eqn. (5081).

The practical consequence of this removal is that the corresponding adjoint solvers can be left out of the process, enhancing the robustness and performance of the adjoint solution process at the expense of accuracy.