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.
- enhance POSIX compliance
- apply distinct colours and dash type for each line
- standardize the frame size to 1200x627
- dynamically replace the title with <function-object-name>/<file-name>
- address underscore character issues
- introduce legend components for tensors
- resolve a bug caused by parentheses in tensor files
BUG: particleTrackProperties: correct the typo (fixes#3050)
- on large memory systems (eg, 6TB) the process information
exceeds an 'int' range, so adjust parsing of the /proc/..
to use int64
ENH: update/modernize OSspecific system information
ENH: minor update of profiling code
- std::string, noexcept, lazier evaluations
STYLE: use direct call of memInfo
- static version of polyMesh::meshDir(), which takes a region name
polyMesh::meshDir(regionName)
vs
polyMesh::regionName(regionName)/polyMesh::meshSubDir
STYLE: use polyMesh::regionName(..) instead of comparing to defaultRegion
STYLE: use getOrDefault when retrieving various -region options
FIX: polyMesh::dbDir() now checks registry name, not full path (#3033)
- in most cases a parallel-consistent order is required.
Even when the order is not important, it will generally require
fewer allocations to create a UPtrList of entries instead of a
HashTable or even a wordList.
- use typeHeaderOk<regIOobject>(false) for some generic file existence
checks. Often had something like labelIOField as a placeholder, but
that may be construed to have a particular something.
- this complements the whichPatch(meshFacei) method [binary search]
and the list of patchID() by adding internal range checks.
eg,
Before
~~~~~~
if (facei >= mesh.nInternalFaces() && facei < mesh.nFaces())
{
patchi = pbm.patchID()[facei - mesh.nInternalFaces()];
...
}
After
~~~~~
patchi = pbm.patchID(facei);
if (patchi >= 0)
{
...
}
- skip loading of fields with -no-internal, -no-boundary
- suppress reporting fields with -no-internal, -no-boundary
- cache loaded volume field for reuse with point interpolation.
Trade off some memory overhead against reading twice.
NOTE: this issue will not be evident with foamToEnsight since there
it only handles cell data *or* point data (not both), so a field is
only ever loaded/processed once.
- since ensight format is always float and also always written
component-wise, perform the double -> float narrowing when
extracting the components. This reduces the amount of data
transferred between processors.
ENH: avoid vtk/ensight parallel communication of empty messages
- since ensight writes by element type (eg, tet, hex, polyhedral) the
individual written field sections will tend to be relatively sparse.
Skip zero-size messages, which should help reduce some of the
synchronization bottlenecks.
ENH: use 'data chunking' when writing ensight files in parallel
- since ensight fields are written on a per-element basis, the
corresponding segment can become rather sparsely distributed. With
'data chunking', we attempt to get as many send/recv messages in
before flushing the buffer for writing. This should make the
sequential send/recv less affected by the IO time.
ENH: allow use of an external buffer when writing ensight components
STYLE: remove last vestiges of autoPtr<ensightFile> for output routines
- consistent with sumOp
ENH: globalIndex with gatherNonLocal tag, and use leading dispatch tags
- useful for gather/write where the master data can be written
separately. Leading vs trailing dispatch tags for more similarity to
other C++ conventions.
- replaced ad hoc handling of formatOptions with coordSetWriter and
surfaceWriter helpers.
Accompanying this change, it is now possible to specify "default"
settings to be inherited, format-specific settings and have a
similar layering with surface-specific overrides.
- snappyHexMesh now conforms to setFormats
Eg,
formatOptions
{
default
{
verbose true;
format binary;
}
vtk
{
precision 10;
}
}
surfaces
{
surf1
{
...
formatOptions
{
ensight
{
scale 1000;
}
}
}
}
Header information now includes, e.g.
f [Hz] vs P(f) [Pa]
Lower frequency: 2.500000e+01
Upper frequency: 5.000000e+03
Window model: Hanning
Window number: 2
Window samples: 512
Window overlap %: 5.000000e+01
dBRef : 2.000000e-05
Area average: false
Area sum : 6.475194e-04
Number of faces: 473
Note: output files now have .dat extension
- include -no-libs option by default, similar to '-lib',
which makes it available to all solvers/utilities.
Add argList allowLibs() method to query it.
- relocate with/no functionObjects logic from Time to argList
itself as argList allowFunctionObjects()
- add libs/functionObjects override handling to decomposePar etc
ENH: report the stream relativeName for IOerrors (see c9333a5ac8)
- consistent with defining IO of int32_t/int64_t and with recent
changes to ensightFile. Using the primitives directly instead of
typedefs to them makes the code somewhat less opaque.
STYLE: qualify format/version/compression with IOstreamOption not IOstream
STYLE: reduce number of lookups when scanning {fa,fv}Solution
STYLE: call IOobject::writeEndDivider as static
thermoTools is a relocation of various existing tools:
- src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/
- src/semiPermeableBaffle/derivedFvPatchFields/
- src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/
ENH: Allwmake: reordering various compilation steps
Co-authored-by: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
- in various situations with mesh regions it is also useful to
filter out or remove the defaultRegion name (ie, "region0").
Can now do that conveniently from the polyMesh itself or as a static
function. Simply use this
const word& regionDir = polyMesh::regionName(regionName);
OR mesh.regionName()
instead of
const word& regionDir =
(
regionName != polyMesh::defaultRegion
? regionName
: word::null
);
Additionally, since the string '/' join operator filters out empty
strings, the following will work correctly:
(polyMesh::regionName(regionName)/polyMesh::meshSubDir)
(mesh.regionName()/polyMesh::meshSubDir)
- previously filtered on the existence of area fields, but with
faMesh::TryNew this is not required anymore.
STYLE: enable -verbose for various parallel utilities (consistency)
- specifies the number of consecutive cells to assign to the same
randomly chosen processor. Can be used to have a less extremely
random distribution for testing possible breaking points.
Eg,
method random;
coeffs
{
agglom 4;
}
- Add finiteArea cellID (actually face ids) / faceLabel and procID
for foamToVTK with -write-ids. Useful when this type of information
is needed.
- 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
- also disables PointData if manifold cells are detected.
This is a partial workaround for volPointInterpolation problems
with handling manifold cells.
- additional verbosity option for conversions
- ignore old `-finite-area` option and always convert available
finiteArea mesh/fields unless `-no-finite-area` is specified (#2374)
ENH: simplify point offset handling for ensight output
- extend writing to include compact face/cell lists
- a try/catch approach is not really robust enough (or even possible)
since read failures likely do not occur on all ranks simultaneously.
This leads to situations where the master has thrown an exception
(and thus exiting the current routine) while other ranks are still
waiting to receive data and the program blocks completely.
Since this primarily affects data conversion routines such as
foamToEnsight etc, treat similarly to lagrangian: check for the
existence of essential files before proceeding or not. This is
wrapped into a TryNew factory method:
autoPtr<faMesh> faMeshPtr(faMesh::TryNew(mesh));
if (faMeshPtr) ...
- the very old 'writer' class was fully stateless and always templated
on an particular output type.
This is now replaced with a 'coordSetWriter' with similar concepts
as previously introduced for surface writers (#1206).
- writers change from being a generic state-less set of routines to
more properly conforming to the normal notion of a writer.
- Parallel data is done *outside* of the writers, since they are used
in a wide variety of contexts and the caller is currently still in
a better position for deciding how to combine parallel data.
ENH: update sampleSets to sample on per-field basis (#2347)
- sample/write a field in a single step.
- support for 'sampleOnExecute' to obtain values at execution
intervals without writing.
- support 'sets' input as a dictionary entry (as well as a list),
which is similar to the changes for sampled-surface and permits use
of changeDictionary to modify content.
- globalIndex for gather to reduce parallel communication, less code
- qualify the sampleSet results (properties) with the name of the set.
The sample results were previously without a qualifier, which meant
that only the last property value was actually saved (previous ones
overwritten).
For example,
```
sample1
{
scalar
{
average(line,T) 349.96521;
min(line,T) 349.9544281;
max(line,T) 350;
average(cells,T) 349.9854619;
min(cells,T) 349.6589286;
max(cells,T) 350.4967271;
average(line,epsilon) 0.04947733869;
min(line,epsilon) 0.04449639927;
max(line,epsilon) 0.06452856475;
}
label
{
size(line,T) 79;
size(cells,T) 1720;
size(line,epsilon) 79;
}
}
```
ENH: update particleTracks application
- use globalIndex to manage original parcel addressing and
for gathering. Simplify code by introducing a helper class,
storing intermediate fields in hash tables instead of
separate lists.
ADDITIONAL NOTES:
- the regionSizeDistribution largely retains separate writers since
the utility of placing sum/dev/count for all fields into a single file
is questionable.
- the streamline writing remains a "soft" upgrade, which means that
scalar and vector fields are still collected a priori and not
on-the-fly. This is due to how the streamline infrastructure is
currently handled (should be upgraded in the future).
GIT: relocate globalIndex (is independent of mesh)
STYLE: include label/scalar Fwd in contiguous.H
STYLE: unneed commSchedule include in GeometricField
