Coupled Aero-Structural Adjoint Analysis with SU2 and External FEM Solvers

[patelha57]

Background

@aa-g and myself have been working to extend the SU2 Python API pysu2 for use in aero-structural MDO with external structural solvers like NASTRAN and externally defined design variables (see issue #1262 and pull request draft #1300).

To set up such a coupled discrete-adjoint sensitivity analysis correctly, we have found that we need a better understanding of the definition of the partial derivatives, adjoint variables, and related terms appearing in the discrete adjoint flow solver of SU2.

Discussion

 A good place to start would be the variables which are written to file by the adjoint solver:

• the adjoint variables (denoted Adjoint_*) in the solution files of SU2_CFD_AD
• the volume and surface sensitivities (denoted Sensitivity_*) in the solution files of SU2_CFD_AD
• the volume and surface sensitivities (denoted Sensitivity_*) in the solution files of SU2_DOT_AD (after SU2_CFD_AD)

Would it be possible to identify or relate these values to the partial derivatives and adjoints variables in the following formulation of the coupled aero-structural adjoint equations, which we have assumed so far in our work?

Theoretical Derivation

Combining the residual equations from the aerodynamic and structural disciplines, the residual of the multidisciplinary system can be expressed as

where q are the flow state variables, u the structural displacements, x the design variables, A the residuals (or pseudo-residual based on fixed-point form) of the governing flow equations, and S the residuals of the governing structural equations. If there are more design variables than functions of interest, it is computationally more efficient to use the adjoint method

where psi is a vector of flow adjoint variables, phi a vector of structural adjoint variables, and I the objective function. Once all adjoint vectors have converged, we can compute the final sensitivities of the objective function as

[TobiKattmann]

Harsh and I now talked a bit about some of these topics. I prob need to structure this summary a bit more but this hopefully contains the relevant points

@patelha57 please attach the overview related to the output of CFD_AD and DOT_AD so we can conclude/document some findings here.

Wrt computing the matrix-vector product dA/dq^T \psi (derivative of residual wrt to flow state times adjoint variables): One of the major points not done was to register the required output (i.e. residuals (instead of the updated flow state which SU2_CFD_AD does usually)) after the DirectRun.

The DirectRun is probably best adapted such the Iterate (Linear solver call) is not done, as this is not required. The call to the seeding (InitializeAdjoints) does not have to be adapted if the correct Output has been registered and if the correct adjoint solution \psi is in the Solution of the adjoint solver.

Another question is how to handle the Extraction of the AD-adjoints of the inputs (i.e. the initial flow state) which holds the mat-vec-product that has to be communicated to OpenMDAO. One could run the IterateDiscAdj function because this would put the correct result into the Solution of the adjoint solver which would prob be easier to extract / expose to OpenMDAO.

Note that OpenMDAO has to provide the current correct \psi value each time calling SU2_CFD_AD stuff because OpenMDAO of course iterates in \psi

Now it is also necessary to compute the derivative dI/dq / dI/dx separately. As the SingleZoneDriver might be easier to digest one needs to adapt the logic there as the SingleZoneDriver always tapes through the fixed-point and then tapes the obj-func computation in one run. We need to be able to record the obj-func evaluation separately. If for the Residual computation the Iterate step is dismissed the correct (restarted) primal values will remain in the Solution. I am not sure one could use the linear combination of mat-vec product and Obj-func derivative (so that this whole thing would not necessary) : therefore I would need to understand how OpenMDAO solves the above system for the adjont equations

The next questions would be the computation of dA/du^Y \psi as it requires the Volume deformation step as well. One can make use of the DEFORM_MESH option in the SU2_CFD_AD` run which already makes the volume deformation step. At least that would be a starting point to check how one could incorporate that.

[patelha57]

In this table, xv are the fluid volume mesh coordinates, xa are the surface mesh coordinates, and x are the design variables defined in SU2 (e.g., FFD box coordinates). The surface sensitivities are output for a subset of those volume nodes, evaluated only at the FSI boundary or the Outer Mold Line (OML).

Additional Notes

SU2 Implementation Outline

I am adding a few more notes about the adjoint implementation in SU2 and specific method calls from our discussions:

A. Taping with CDiscAdjSinglezoneDriver::SetRecording()

1. Register q
2. Compute residual A
3. Register output A (i.e., set in AD output index container)

B. Matrix-vector product with CDiscAdjSinglezoneDriver::Run()

1. Initialize AD adjoint of A with psi as the seed vector with InitializeAdjoint()
2. Evaluate tape with AD::ComputeAdjoint()
3. Extract dA/dq^T psi with IterateDiscAdj()

OpenMDAO Implementation

The system of coupled adjoint equations is iteratively solved via a Linear Block Gauss-Seidel (LBGS) solver. @TobiKattmann, is there anything in particular you want me to expand on? Please let me know and I can provide more details.

[TobiKattmann]

Thanks!
Just one thing: For the surface_sens* output of SU2_DOT_AD I would personally introduce a surface mesh coordinate vector e.g. x_a. I know, that you know what the surface_sens of SU2_DOT_AD actually means, but as it is written down the term is equivalent to the second term of the surface_adjoint from SU2_CFD_AD (except to the transpose) ... and I would avoid this equal terminology.

[patelha57]

That's a good point. I've updated the table to reflect that change.

[aa-g]

@patelha57 and I have been working on implementing a new driver (based on the CDiscAdjSinglezoneDriver) that handles adjoint solutions and calculations of sensitivities based on the residual formulation as above (in contrast to the fixed-point iteration as it is currently done).

The key changes and differences so far can be summarized as follows:

1. The outputs to be registered are now the residuals, instead of the solution (i.e. the conservative state):

   For this, we added a member Residual to CVariable, which is registered using a method RegisterResidual much like it is done for the Solution.

2. The direct run no longer needs to include the iteration on the linear system solution, and now only runs the calculation of the residuals.

   For this, we added a method SetResidual to CSolver, which at the moment is only actually implemented in CFVMFlowSolverBase, where it calls the protected method EdgeFluxResidual to evaluate the residual at each vertex. Once this is done, the value of LinSysRes is used to set the Residual which was registered as an output.

3. The next step would be to split the recording of the residual and objective function computation, both w.r.t. the solution ("main recording") and the geometry ("secondary recording"), resulting in four separate recordings: dA/dq, dA/du_a, dI/dq, dI/du_a.

   In the CSinglezoneDriver, the method SetRecording distinguishes between the main and secondary recordings by the input, specified as kind_recording which is of type RECORDING. One possibility could be to add a similar logic for the output (i.e. distinguish between RECORDING_INPUT and RECORDING_OUTPUT).

Since we are both inexperienced with the AD tool (and in general, contributing to SU2), we're always looking to make sure that we are taking the right approach in our implementation. So, @pcarruscag, @TobiKattmann, please let us know if any of these ideas seem concerning to you :) or if you have suggestions for improvement

[pcarruscag]

Hi Harsh, I think that is going the wrong way, why copy the residuals into CVariable if you have them in LinSysRes...
Have you looked at the Newton-Krylov implementation? In particular you have CNewtonIntegration::ComputeResiduals.
You also need the contribution of the boundaries to the residuals, these happen in SpaceIntegration (which the method above calls).
Forget about the DiscAdjDriver, you will end up having to modify the iteration that way, you have this in the CDiscAdjSolver:

void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) {

  /*--- Register variables as output of the solver iteration. Boolean false indicates that an output is registered ---*/

  direct_solver->GetNodes()->RegisterSolution(false);
}

Just change it to register LinSysRes, it should work to build a conventional adjoint for Euler problems.
Then you can expose some method in the wrapper that seeds the residual and evaluates the tape.
Getting the gradient of the OF is a bit more difficult, but I think it would be enough to hack the implicit Euler iteration into not solving the linear system.
You can see what steps need to be captured in the recording of the OF in the DiscAdjMultizoneDriver.

[pcarruscag]

Sorry Alessandro I misread who the author was.

[patelha57]

Regarding "Forget about the DiscAdjDriver", are you saying creating a new driver based on residuals is not the correct approach? If that is the case, how do we bypass the call to direct_iteration->Iterate in CDiscAdjSinglezoneDriver::DirectRun()?

Also can you please elaborate on the implementation of having 4 separate recordings (for residuals and obj. functions, each w.r.t to flow states and surface displacements) as we outlined above? How do we handle this because each needs to be addressed separately?

[pcarruscag]

Creating a "new" driver will result in duplication, at most you need to override the strategy for recording.
You don't need 4 necessarily, you will record the computation of residuals but not the update of the solution, thus you might be able to compute both the residuals and the OF (the latter based on "un"-updated valued) and then depending on which outputs you seed you get either the gradient of the OF or the jacobian-vector product.
You only need 2 recordings (state and coords) for efficiency (if you record more operations the tape gets more expensive to evaluate).
The reason you should look at the multizone rather than singlezone is because of the mesh deformation. You want to be able to evaluate that cross term separately from the block diagonal term, because it will include mesh deformation and that is expensive.
I or Ole or Tobi can give you the walkthrough next Wednesday.