The solution of the QP subproblem can become quite expensive, especially
for cases with many design variables (e.g. topology optimisation).
A (potentially dense) matrix with the size of the design variables is
solved using a matrix-free CG solver. The convergence speed greatly
depends on the used preconditioner. This commit adds
preconditioner-vector products based on the L-BFGS inverse Hessian and,
more importantly, a preconditioner computed using the Sherman-Morrison
formula. The latter is applicable here since the LHS of the QP problem
is computed as the sum of rank-2 L-BFGS updates, a sum of rank-1 updates
(as many as the flow-related constraints) and a diagonal matrix
depending on the bound constraints.
Additionally, the QP subproblem could have no feasible points. To relax
this, constraints can be applied gradually through the
targetConstraintReduction enty (typical value of 0.1 for topology
optimisation).
Parts of the adjoint optimisation library were re-designed to generalise
the way sensitivity derivatives (SDs) are computed and to allow easier
extension to primal problems other than the ones governed by
incompressible flows. In specific:
- the adjoint solver now holds virtual functions returning the part of
SDs that depends only on the primal and the adjoint fields.
- a new class named designVariables was introduced which, apart from
defining the design variables of the optimisation problem and
providing hooks for updating them in an optimisation loop, provides
the part of the SDs that affects directly the flow residuals (e.g.
geometric variations in shape optimisation, derivatives of source
terms in topology optimisation, etc). The final assembly of the SDs
happens here, with the updated sensitivity class acting as an
intermediate.
With the new structure, when the primal problem changes (for instance,
passive scalars are included), the same design variables and sensitivity
classes can be re-used for all physics, with additional contributions to
the SDs being limited (and contained) to the new adjoint solver to be
implemented. The old code structure would require new SD classes for
each additional primal problem.
As a side-effect, setting up a case has arguably become a bit easier and
more intuitive.
Additional changes include:
---------------------------
- Changes in the formulation and computation of shape sensitivity derivatives
using the E-SI approach. The latter is now derived directly from the
FI approach, with proper discretization for the terms and boundary
conditions that emerge from applying the Gauss divergence theorem used
to transition from FI to E-SI. When E-SI and FI are based on the same
Laplace grid displacement model, they are now numerically equivalent
(the previous formulation proved the theoretical equivalence of the
two approaches but numerical results could differ, depending on the
case).
- Sensitivity maps at faces are now computed based (and are deriving
from) sensitivity maps at points, with a constistent point-to-face
interpolation (requires the differentiation of volPointInterpolation).
- The objective class now allocates only the member pointers that
correspond to the non-zero derivatives of the objective w.r.t. the
flow and geometric quantities, leading to a reduced memory footprint.
Additionally, contributions from volume-based objectives to the
adjoint equations have been re-worked, removing the need for
objectiveManager to be virtual.
- In constrained optimisation, an adjoint solver needs to be present for
each constraint function. For geometric constraints though, no adjoint
equations need to solved. This is now accounted for through the null
adjoint solver and the geometric objectives which do not allocate
adjoint fields for this kind of constraints, reducing memory
requirements and file clutter.
- Refactoring of the updateMethod to collaborate with the new
designVariables. Additionally, all updateMethods can now read and
write restart data in binary, facilitating exact continuation.
Furthermore, code shared by various quasi-Newton methods (BFGS, DBFGS,
LBFGS, SR1) has been organised in the namesake class. Over and above,
an SQP variant capable of tackling inequality constraints has been
added (ISQP, with I indicating that the QP problem in the presence of
inequality constraints is solved through an interior point method).
Inequality constraints can be one-sided (constraint < upper-value)
or double-sided (lower-value < constraint < upper-value).
- Bounds can now be defined for the design variables.
For volumetricBSplines in specific, these can be computed as the
mid-points of the control points and their neighbouring ones. This
usually leads to better-defined optimisation problems and reduces the
chances of an invalid mesh during optimisation.
- Convergence criteria can now be defined for the optimisation loop
which will stop if the relative objective function reduction over
the last objective value is lower than a given threshold and
constraints are satisfied within a give tolerance. If no criteria are
defined, the optimisation will run for the max. given number of cycles
provided in controlDict.
- Added a new grid displacement method based on the p-Laplacian
equation, which seems to outperform other PDE-based approaches.
TUT: updated the shape optimisation tutorials and added a new one
showcasing the use of double-sided constraints, ISQP, applying
no-overlapping constraints to volumetric B-Splines control points
and defining convergence criteria for the optimisation loop.
- avoids redundant dictionary searching
STYLE: remove dictionary lookupOrDefaultCompat wrapper
- deprecated and replaced by getOrDefaultCompat (2019-05).
The function is usually specific to internal keyword upgrading
(version compatibility) and unlikely to exist in any user code.
- 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.
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.
- bundles frequently used 'gather/scatter' patterns more consistently.
- combineAllGather -> combineGather + broadcast
- listCombineAllGather -> listCombineGather + broadcast
- mapCombineAllGather -> mapCombineGather + broadcast
- allGatherList -> gatherList + scatterList
- reduce -> gather + broadcast (ie, allreduce)
- The allGatherList currently wraps gatherList/scatterList, but may be
replaced with a different algorithm in the future.
STYLE: PstreamCombineReduceOps.H is mostly unneeded now
- controlPointsDefinition is now controled by a class with
runTimeSelection.
- Added a new controlPointsDefinition option that translates, rotates
and scales a given box. The required entries have the same meaning as
in the Paraview 'Transform' filter, facilitating the transition between the
visual placement of control boxes (e.g. in Paraview) and their setup
in the code.
- Improved performance during the parameterization, sensitivity
computation and grid displacement phases by re-using already computed
basis functions.
The if(Pstream::master()) clause in NURBS3DVolume::writeCpsInDict() was
causing the fileName of the regIOobject not to be allocated in all
processors, giving problems when masterUncollatedFileOperation::masterOp
was called by collatedFileOperation::writeObject for the mkDirOp.
- with '&&' conditions, often better to check for non-null autoPtr
first (it is cheap)
- check as bool instead of valid() method for cleaner code, especially
when the wrapped item itself has a valid/empty or good.
Also when handling multiple checks.
Now
if (ptr && ptr->valid())
if (ptr1 || ptr2)
instead
if (ptr.valid() && ptr->valid())
if (ptr1.valid() || ptr2.valid())
- previously introduced `getOrDefault` as a dictionary _get_ method,
now complete the transition and use it everywhere instead of
`lookupOrDefault`. This avoids mixed usage of the two methods that
are identical in behaviour, makes for shorter names, and promotes
the distinction between "lookup" access (ie, return a token stream,
locate and return an entry) and "get" access (ie, the above with
conversion to concrete types such as scalar, label etc).
When more than one volumetric B-Splines control boxes are present, the
sensitivity constituents corresponding to the non-active design
variables were not bounded(zeroed) correctly. The resultant
sensitivities, used in the optimization, were bounded correctly, so this
was more a bug pertaining to the output file of the sensitivities rather
than a functional one.
- Failed due to double*Matrix<float> multiplication.
Style changes
- use SquareMatrix with Identity on construction
- use Zero in constructors
- remove trailing space and semi-colons
The adjoint library is enhanced with new functionality enabling
automated shape optimisation loops. A parameterisation scheme based on
volumetric B-Splines is introduced, the control points of which act as
the design variables in the optimisation loop [1, 2]. The control
points of the volumetric B-Splines boxes can be defined in either
Cartesian or cylindrical coordinates.
The entire loop (solution of the flow and adjoint equations, computation
of sensitivity derivatives, update of the design variables and mesh) is
run within adjointOptimisationFoam. A number of methods to update the
design variables are implemented, including popular Quasi-Newton methods
like BFGS and methods capable of handling constraints like loop using
the SQP or constraint projection.
The software was developed by PCOpt/NTUA and FOSS GP, with contributions from
Dr. Evangelos Papoutsis-Kiachagias,
Konstantinos Gkaragounis,
Professor Kyriakos Giannakoglou,
Andy Heather
[1] E.M. Papoutsis-Kiachagias, N. Magoulas, J. Mueller, C. Othmer,
K.C. Giannakoglou: 'Noise Reduction in Car Aerodynamics using a
Surrogate Objective Function and the Continuous Adjoint Method with
Wall Functions', Computers & Fluids, 122:223-232, 2015
[2] E. M. Papoutsis-Kiachagias, V. G. Asouti, K. C. Giannakoglou,
K. Gkagkas, S. Shimokawa, E. Itakura: ‘Multi-point aerodynamic shape
optimization of cars based on continuous adjoint’, Structural and
Multidisciplinary Optimization, 59(2):675–694, 2019
