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
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., PapoutsisKiachagias, 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), 15231542. 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.
 Adjoint solvers write/read their sensitivity derivatives under the 'uniform' folder, to avoid potential loss of accuracy due to I/O.
 Volumetric BSplines 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 misssetups.
 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 ESI 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 rewritten to reduce this peak memory consumption.
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).