ENH: adjoint code review

Summary

An enhancement and review of the adjoint optimisation library, featuring

• the adjoint to the k-ω SST turbulence model

• easier and more accurate restarted runs

• reduced turnaround times

• reduced peak memory consumption

  Resolved bugs

Partially resolves #2502, allowing the update of the nearWallDist field throughout an optimisation loop.

Details of the code enhancement

Added the adjoint to the k-ω SST turbulence model for incompressible flows. This allows for avoiding the "frozen turbulence" assumption when using this turbulence model within an optimisation loop or when computing sensitivity maps. Making the "frozen turbulence" assumption, i.e. assuming that the turbulent viscosity field does not change when the shape changes throughout the optimisation, can lead to erroneously computed sensitivity derivatives. As an example, the drag sensitivity maps computed on the surface of the Ahmed body using the "frozen turbulence" assumption and the fully differentiated k-ω SST model are depicted below

featuring areas where the sign of the sensitivity map changes when turbulence is not differentiated.

Taking this a step further and executing optimisations using "frozen turbulence" (FT) and "differentiated turbulence" (DT) highlights the need for differentiating the k-ω SST model in this case

The DT approach reduces drag by more the 6% within 20 optimisation cycles whereas the optimisation based on the FT assumption diverges after the fourth cycle.

The work is based on [1] with changes in the discretisation of a number of differential operators and the formulation of the adjoint to the wall functions employed by the primal model.

Source code

$FOAM_SRC/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST

Examples

$FOAM_TUTORIALS/incompressible/adjointOptimisationFoam/shapeOptimisation/naca0012/kOmegaSST/lift $FOAM_TUTORIALS/incompressible/adjointOptimisationFoam/shapeOptimisation/sbend/turbulent/kOmegaSST/opt

References

[1] Kavvadias, I., Papoutsis-Kiachagias, E., Dimitrakopoulos, G., & Giannakoglou, K. (2014). The continuous adjoint approach to the k–

$omega$

SST turbulence model with applications in shape optimization. Engineering Optimization, 47(11), 1523-1542. https://doi.org/10.1080/0305215X.2014.979816

Attribution

The initial implementation, as described in [1], was performed by Dr. Ioannis Kavvadias.

Details of the code review

Three main aspects are addressed

• Continuation

  Restarting a (partially) already ran optimisation is now easier and more accurate.

  1. Adjoint solvers write/read their sensitivity derivatives under the 'uniform' folder, to avoid potential loss of accuracy due to I/O.
  2. Volumetric B-Splines control points of each optimisation cycle are written under 'uniform' and support also binary I/O. As a consequence, the controlPointsDefinition in constant/dynamicMeshDict does not need to be changed to 'fromFile' anymore in order to perform the continuation, removing a potential source of miss-setups.
  3. The adjoint grid displacement field (ma) is now appended by the name of the adjoint solver, if more than one exist. This facilitates continuation since, before the change, only the ma field of the last adjoint solver was written to file. No changes to fvSchemes/fvSolution are necessary though.

• Reduced turnaround times

  A number of changes to reduce the solution turnaround time of the adjoint equations (see bac1d8ba for details). In brief, these include caching of some expensive but constant quantities and removing some coding shortcomings. All cases should see some benefit (e.g. around 9.5% reduction in the turnaround time of the adjoint solver for the motorbike tutorial) but the latter can be even more pronounced in cases with many outlet boundaries.

• Reduced peak memory consumption

  The peak memory consumption of the adjoint code is observed during the computation of sensitivity derivatives, when using either the FI or the E-SI approach, due to the need of manipulating a number of volTensorFields necessary to compute the multiplier of the spatial gradient of the grid sensitivities. This part of the code has been re-written to reduce this peak memory consumption.

  Risks

  No changes to the user input, apart from the second item in the 'Continuation' listing above, which should make the setup of continuation runs easier/more straightforward.

  Additionally, sensitivities are now computed at the end of each adjoint solver, instead of being computed when all adjoint solvers are finished, but this should not affect the code behaviour and/or setup in any way.

  Finally, even though 5937c37a resolves the main issue in #2502, it practically disables caching of gradients after the first mesh update within an optimisation loop (see #2502 for a discussion).

changed milestone to %v2112

added community enhancement labels

assigned to @andy

added 1 commit

• ef9df63f - TUT: added a tutorial showcasing the use of the nutSqr objective

Compare with previous version

marked this merge request as ready

changed title from Draft: adjoint code review to ENH: adjoint code review

changed milestone to %v2206

marked this merge request as draft

added 399 commits

• ef9df63f...6e644a1e - 377 commits from branch develop
• d7846b90 - ENH: add the infrastructure for computing and utilising
• 58b4df0c - ENH: added the adjoint to the kOmega SST turbulence model
• 6dfda354 - TUT: added a tutorial for the adjointkOmegaSST turbulence model
• 58c047b0 - COMP: Minor changes for clean complilation
• 0b455bbe - ENH: changes towards a machine-accurate continuation
• 29c834c3 - ENH: refactoring of the sensitivity classes
• 72969eee - ENH: enable writing volumetric B-Splines control points in binary
• e7ead700 - ENH: the adjoint grid displacement field (ma)
• 828bd8c8 - ENH: the adjoint eikonal equation grabs the epsilon value
• 04246e83 - ENH: adjustments to the efficiency of the adjoint code
• 91438ac2 - ENH: changes reducing the peak memory consumption of shape sensitivities
• 3fb4c4ab - ENH: when using (E)SI sensitivities and a symmetry(Plane) is included
• 757cdad8 - ENH: made the boundControlPointMovement method of
• eb3e95a0 - TUT: updated the BFGS continuation tutorial
• 2d7daca5 - TUT: added a tutorial showcasing the use of the nutSqr objective
• 9c5dc08e - ENH: useSolverNameForFields is now set to true automatically
• 741bdf6b - ENH: primal and adjoint solvers utilised by adjointOptimisationFoam
• 566562b4 - TUT: added a tutorial checking/showcasing the usage of
• 67bd3b04 - BUG: the functions reporting the existance of turbulence fields
• 8ea77f25 - Merge branch 'feature-adjointkOmegaSST' into adjoint-code-review
• e81b4937 - ENH: exposed wallFunctionCoefficients in nutWallFunction
• 032c8665 - ENH: use of wallFunctionCoefficients in adjointkOmegaSST

Compare with previous version

Hi @andy

I have pushed some more updates to the adjoint framework and the adjoint to k-omega SST to this branch. All tutorials and some additional tests seem to run OK, with the exception of one (see point 2 below). Points that might require your attention:

1. Automated addition of divSchemes When more than one adjoint solvers exist (e.g. when constraints are present) or when we are running a multi-point optimization problem (i.e. more than one primal and adjoint solvers are present), the names of the primal/adjoint fields are appended by the name of the corresponding solver. This creates the need to change the entries in divSchemes,something that might easily be forgotten by the user and lead to poor convergence. One solution would be to always expect the same entry in divSchemes, irrespective of the names of the fields, by using fvm::div(phi,yourField,"yourEntryHere"), but this prevents us from giving different discreritzation schemes for different primal/adjoint solvers, which is a bit restrictive. Instead, I used the non-const reference to divSchemes, added by @mark in 3a1a160c, to add the necessary schemes on the fly, based on a fall-back scheme that corresponds to the entry we would have, if only one primal/adjoint solver was present (see 741bdf6b). This works fine, but there is one possible pitfall: if the user edits and writes fvSchemes manually, then the aforementioned divSchemes added by the code, which reside in memory only, will be lost. One solution to this problem is to write fvSchemes through the code when the divSchemes are added. When attempting to do this by grabbing a reference to fvSchemes, which is an IOdictionary, and calling write through regIOobject, what I get is a call to objectRegistry.writeObject() instead of IOdictionary.writeObject(), which is kind of weird, at least to me (see commented lines 108-110 in variablesSet.C; compiled with gcc 9.3.0). Casting doesn't seem to help either. I have used a work-around of creating a local fvSchemes IOdictionary and writing it, but you might want a more elegant solution. This could be based on schemesLookup::writeDicts(), but there seems to be no hook to it (for instance, through writeData); even when I added that (locally, not pushed), this still wasn't effective since objectRegistry seems to bypass IOdictionary here. This isn't major by any means, just wanted you opinion.

2. The tutorial under $FOAM_TUTORIALS/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike fails when attempting to solve the smoothing equation of the sensitivity map, based on an faMesh. The faMesh seems to be generated without a problem but the code fails inFoam::faMesh::calcEdgeTransformTensors() (see attached log for a trace log.adjointOptimisationFoam).The problem is not related to what is added in this merge request and should be reproducible in develop. It didn't exist in v2112, so it must have been introduced some time between then and now in the develop branch. I cannot pin-point it further for the moment since we are basing our in-house developments on master and I haven't merged with develop since v2112. @mark, any thoughts? If you prefer, I can open a new issue for this one.