ENH: adjoint code review

See merge request !518

Vaggelis Papoutsis authored 2 years ago and Andrew Heather committed 2 years ago

Necessary when grad(Ua) is cached

Vaggelis Papoutsis authored 2 years ago and Andrew Heather committed 2 years ago

Partially fixes #2502 (see discussion there)

Vaggelis Papoutsis authored 2 years ago and Andrew Heather committed 2 years ago

in RASModelVariables were doing this by checking whether the
corresponding pointer was allocated. In some cases, however, even if the
field does not exist, the pointer is not null, leading to the wrong
output. Made the correspding functions virtual and overwritten their
return values in the derived classes. Kept the initial implementation in
base to facilitate the clone function.

Vaggelis Papoutsis authored 2 years ago and Andrew Heather committed 2 years ago

in cases with more than one primal or adjoint solvers

TUT: removed all occurances of useSolverNameForFields
from the optimisation tutorials since it is now set
automatically.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

to be in line with the changes in 00f46a2c52 and added the same tutorial
which performs the optimisation in one-go, for comparison.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

volBSplinesBase const, since the no private/protected
variable is changed there.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

in the sensitivity patches, symmetry::evaluate() needs access to the
internalField which does exist, leading to wrong memory access.

Fixed by specifying a calculated type fvPatchField for all patches when
creating a boundaryField<Type>

Using a symmetry(Plane) as a sensitivity patch is quite rare and
borderline wrong, but this provides a fix nonetheless.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

The multiplier of grad(dxdb) is a volTensorField which, by itself, is
memory consuming. The function computing it though was sloppy in terms
of memory management, constituting the peak memory consumption during an
adjoint optimisation. Initial changes to remedy the problem include the
deallocation of some of the volTensorFields included in the computation
of grad(dxdb) once unneeded, the utilisation of volSymmTensorFields
instead of volTensorFields where possible and avoiding allocating some
unnecessary intermediate fields.

Actions to further reduce memory consumption:

- For historical reasons, the code computes/stores the transpose of
  grad(dxdb), which is then transposed when used in the computation of
  the FI or the ESI sensitivity derivatives.  This redundant
  transposition can be avoid, saving the allocation of an additional
  volTensorField, but the changes need to permeate a number of places in
  the code that contribute to grad(dxdb) (e.g.  ATC, adjoint turbulence
  models, adjoint MRF, etc).
- Allocation of unnecessary pointers in the objective class should be
  avoided.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

- ATCstandard, ATCUaGradU:
  the ATC is now added as a dimensioned field and not as an fvMatrix
  to UaEqn. This get rid of many unnecessary allocations.

- ATCstandard:
  gradU is cached within the class to avoid its re-computation in
  every adjoint iteration of the steady state solver.

- Inlined a number of functions within the primal and adjoint solvers.
  This probably has a negligible effect since they likely were inlined
  by the compiler either way.

- The momentum diffusivity at the boundary, used by the adjoint boundary
  conditions, was computed for the entire field and, then, only the
  boundary field of each adjoint boundary condition was used. If many
  outlet boundaries exist, the entire nuEff field would be computed as
  many times as the number of boundaries, leading to an unnecessary
  computational overhead.

- Outlet boundary conditions (both pressure and velocity) use the local
  patch gradient to compute their fluxes. This patch gradient requires
  the computation of the adjacent cell gradient, which is done on the
  fly, on a per patch basis. To compute this patch adjacent gradient
  however, the field under the grad sign is interpolated on the entire
  mesh. If many outlets exist, this leads to a huge computational
  overhead. Solved by caching the interpolated field to the database and
  re-using it, in a way similar to the caching of gradient fields (see
  fvc::grad).

WIP: functions returning references to primal and adjoint boundary
fields within boundaryAdjointContributions seem to have a non-negligible
overhead for cases with many patches. No easy work-around here since
these are virtual and cannot be inlined.

WIP: introduced the code structure for caching the contributions to
the adjoint boundary conditions that depend only on the primal fields
and reusing. The process needs to be completed and evaluated, to make
sure that the extra code complexity is justified by gains in
performance.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

from the equivalent entry of the primal eikonal equation solver, unless
specified in the corresponding adjoint dictionary.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

is now appended by the name of the adjoint solver, if more than one
exist. This was necessary for an accurate continuation since, before
these changes, only the ma field of the last solver was written. As a
result, when restarting the first adjoint solver was reading the ma
field of the last one. No changes are needed in fvSolution and fvSchemes
w.r.t. the previous code version.

Vaggelis Papoutsis authored 4 years ago and Andrew Heather committed 2 years ago

as a step towards machine-accuracy continuation of the optimisation
loop.

Additionally, control points are now written under the time/uniform
folder, to be in-line with rest of the code structure for continuation.
As a side-effect, the controlPointsDefinition in
constant/dynamicMeshDict does not need to be changed to 'fromFile'
anymore in order to perform the continuation. The 'fromFile' option is
still valid if the user wants to supply the control points manually but,
as with all other controlPointsDefinitions, it will be disregarded if the
proper file exists under the time/uniform/volumetricBSplines folder.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

Before the commit, the sensitivity classes were receiving references of
the (incompressible) primal and adjoint variables.  However, if
additional physics was added (energy equation, multiphase, etc), the
infrastructure wasn't convenient for accommodating (new terms in the FI
and E-SI formulations, new terms in the sensitivity map, etc).

Now, the sensitivity classes receive a reference to an
incompressibleAdjointSolver and receive the terms for the FI and
sensitivity maps through there. The latter is still WIP.

Modified adjointSimple to incorporate these changes as well.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

Each solver now writes its sensitivity derivatives to its dictionary,
enabling also a binary format. If present, the sensitivities are then
re-read from the dictionary, avoiding thus possible loss of information
due to re-computation.

As a side-effect, sensitivities are computed after the completion of
each adjoint solver, instead of being computed after all adjoint solvers
have been completed.

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

- computation of lift sensitivity derivatives for the NACA0012 airfoil
- optimisation of an s-bend duct for min. total pressure losses

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

for incompressible flows. The typical convention of appending the primal
field name with 'a' to form the adjoint field is followed for the
adjoint turbulent kinetic energy (i.e. 'ka') but since this would produce
an ugly variable name for the adjoint to omega (i.e. omegaa), the latter
is abbreviated to 'wa'.

The work is based on

\verbatim
    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
\endverbatim

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.

Regarding the latter, the code assumes (and differentiates) the default
behaviour of nutkWallFunction (i.e. nutWallFunction::blendingType::STEPWISE)
and omegaWallFunction (i.e. omegaWallFunction::blendingType::BINOMIAL2).

Due to the availability of a number of terms required for the
formulation of the wall function for ka, the latter is implemented
within adjointkOmegaSST itself, with contributions from objective functions
implemented within kaqRWallFunction. Wall functions for wa are
implemented within waWallFunction.

The initial implementation of the above-mentioned reference was
performed by Dr. Ioannis Kavvadias

Vaggelis Papoutsis authored 3 years ago and Andrew Heather committed 2 years ago

the Jacobian of an objective function, defined at the boundary, wrt nut
and gradU. Also modified the current objectives that include such
contributions

- update the area-centres processor/processor information as part of
  faMesh::init() after all of the global data and geometry data is
  setup.

- improve flattenEdgeField helper to properly handle empty patches.
  This change removes the false fails when testing edge-centre
  redistribution (FULLDEBUG mode).

TUT: add filmPanel (rivulet) tutorial

- include constant/faMesh cleanup (cleanFaMesh) as part of standard
  cleanCase

- simplify cleanPolyMesh function to now just warn about old
  constant/polyMesh/blockMeshDict but not try to remove anything

- cleanup cellDist.vtu (decomposePar -dry-run) as well

ENH: foamRunTutorials - fallback to Allrun-parallel, Allrun-serial

TUT: call m4 with file argument instead of redirected stdin

TUT: adjust suffixes on decomposeParDict variants

Reverted edfault behaviour

Explicitly unset the mesh.moving flag whilst doing the topo change.

Changed behaviour to make it the default

ENH: binField: new function object for data binning

See merge request !433

Andrew Heather authored 2 years ago and Andrew Heather committed 2 years ago

- enables runtime selection of operand coefficients by 'coefficients' entry
- removes binning - now handled using the new 'binField' FO

Co-authored-by: Kutalmis Bercin <kutalmis.bercin@esi-group.com>

Andrew Heather authored 2 years ago and Andrew Heather committed 2 years ago

The new 'binField' function object calculates binned data,
where specified patches are divided into segments according to
various input bin characteristics, so that spatially-localised
information can be output for each segment.

Co-authored-by: Kutalmis Bercin <kutalmis.bercin@esi-group.com>

- previously only defined for cell-cutting version, now for
  iso-surface version too

TUT: remove old transform/coordinateSystem syntax

- simpler to write for sampled cutting planes etc.
  For example,

      slice
      {
          type        cuttingPlane;
          point       (0 0 0);
          normal      (0 0 1);
          interpolate true;
      }

  instead of

      slice
      {
          type        cuttingPlane;
          planeType   pointAndNormal;
          pointAndNormalDict
          {
              point   (0 0 0);
              normal  (0 0 1);
          }
          interpolate true;
      }

STYLE: add noexcept to some plane methods