Newton-Krylov primal iterations & Krylov discrete adjoint

[pcarruscag]

Proposed Changes

In a nutshell, we still solve a linear system to update the solution, but the approximate Jacobian is replaced by matrix-free products with the "real" Jacobian which are obtained by finite differences (each product requires the computation of the residual for a perturbed solution).
This linear system is much more ill-conditioned, and so in addition to the classical linear preconditioners, it is possible to use the traditional linear system as the preconditioner (you gotta love the kind of stuff SU2 lets you get away with sometimes).
This nesting makes each iteration potentially very expensive, it is crucial not to "over solve" the linear systems, and to use adaptive CFL to run always at the highest CFL possible.

Notwithstanding the preconditioner madness, this also works for explicit solvers, but that may only be good enough for dual time methods with relatively small time steps.

This seems to be good at converging problems that would otherwise not get past a limit cycle.
I do not expect it to be faster when that is not a problem.

Some references:
https://arc.aiaa.org/doi/abs/10.2514/6.2021-0857
https://arc.aiaa.org/doi/10.2514/1.J052255

Example config for the RANS ONERA M6 in the repo: turb_ONERAM6_nk.txt
It contains the explanation and recommendations for the different parameters.

The initial idea was to couple all solvers (flow, turbulence, etc.) but that makes the problem even more ill-conditioned and so I took a step back (as that also makes the implementation simpler...).
The skeleton of the coupled iteration exists in commit c97221b and it was functional... if anyone wants to take it further, it's there.

PR Checklist

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags, or simply --warnlevel=2 when using meson).
• My contribution is commented and consistent with SU2 style.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp) , if necessary.

[oleburghardt]

Hi Pedro, seems like nice stuff is coming into being here.. 👍

just had a look at both papers, especially the second. From reading through the new integration class (and assuming you are working towards Newton-like iterations for FSI), could you point me to the coupled part? I'm expecting one has to assemble a "bigger" residual across the disciplines (or a routine that has the same functionality as a multiplication with it, like Algorithm 4 on page 943 in the paper from Kenway, Kennedy and Martins or what you are doing in CNewtonIntegration::MatrixFreeProduct, I guess) :-)

[pcarruscag]

Hi Ole, the first paper is actually closer to what I am aiming for with this implementation.
Mostly I am looking for an alternative to multigrid to get good convergence of the fluid side, that remains effective when many MPI ranks are used.
To apply something like this to multiphysics would require implementing it at the driver level, to include the transfer routines. However, I think for FSI this would not be competitive due to elasticity-based mesh deformation being so expensive compared to one fluid iteration. Martins et. al. use an explicit deformation method that I presume to be very quick as it operates on block structured grids.
The "coupled" part in this title refers to flow, turbulence, and other scalars, and even that I am not sure if it is a good idea :)

[bigfooted]

This is a great new development! We would definitely need this for stiffer equations. At the moment most of our combustion models have 'OK' convergence, but this could really improve things! Is it functional for the incompressible flow testcases?

[pcarruscag]

In principle yes, but there are a few complications for coupling everything which made me take a step back:

• The way the different solvers treat the linear system solution is not consistent, e.g. SST does not use it directly, which means the Jacobian-vector product would not be consistent.
• Scaling, this uses finite differences to compute a directed derivative, if the variables have very different magnitudes the error will be large.
• Conditioning, the residual Jacobian is very ill-conditioned it is difficult to solve systems with it at high CFL, and if CFL is low there is not much point in having the true Jacobian...
• Speed, each matrix-free product means one evaluation of the convective and viscous residuals, for schemes that are not vectorized yet this costs a lot which means the quasi-Newton approach will be faster even if it takes more iterations.

The current implementation seems capable of pushing through some limit cycles that would otherwise just last forever.
But coupling everything is a bigger bite than I can take right now.

[pcarruscag]

I also moved the MPI comms of CSysVectors out of the matrix, in time we can make this even more generic to communicate anything that "looks like" a CSysVector which would allow us to communicate solution variables without having to define enum types for them.

[pcarruscag]

If there are no major comments I would add a testcase and merge this for the next release (and address any post mortem comments as they appear).

[pcarruscag]

The "Krylov discrete adjoint" part is for the inner iterations of multizone problems.
It is based on manipulating the fixed point lambda = Ext + lambda (dx G) into lambda (dx G - I) = -Ext, and calling the LHS a "product" which is passed to GMRES.
This is not a strictly valid thing to do since there is a kind of flexible preconditioner in G, but it works better than the quasi-Newton thing (which is still available).
This manipulation is not possible for single zone because the RHS of the adjoint problem (i.e. the OF gradient) is not separate from the adjoints of the iteration (G).

[oleburghardt]

Please keep the quasi-Newton thing for a while, I'm still building the "real"-Newton thing (through a class inheriting from CQuasiNewtonInvLeastSquares; one would have to figure out a common base class). Interestingly the real one requires the computation of Krylov subspaces, too, so there might be new or other overlaps..

[pcarruscag]

Of course, I'm a hoarder of features, I keep everything I can.

[oleburghardt]

Alright. From what I see, most of the implementation happens within the new integration class with little interference with the solver code. I'd have to read carefully through the first paper but I can also do it post mortem ;-)
I'll create a WIP for the Newton corrector in the next few days so that you can check what I'm up to there. I think it's quite different from what is done here though :-)