Commit 230d7681 authored by Vaggelis Papoutsis's avatar Vaggelis Papoutsis Committed by Andrew Heather
Browse files

ENH: added functionality for smoothing the sensitivity derivatives

A Helmholtz-like filter is applied to the original field of sensitivity
derivatives. The corresponding PDE is solved on the sensitivity patches,
using the finite area infrastructure. A smoothing radius is needed,
which is computed based on the average 'length' of the boundary faces,
if not provided by the user explicitly.

If an faMesh is provided, it will be used; otherwise it will be created
on the fly based on either an faMeshDefinition dictionary in system or
one constructed internally based on the sensitivity patches.
parent 2880a543
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/finiteArea/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
......@@ -14,6 +15,7 @@ EXE_INC = \
LIB_LIBS = \
-lfiniteVolume \
-lfiniteArea \
-lmeshTools \
-lsurfMesh \
-lsampling \
......
......@@ -5,8 +5,8 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2020 PCOpt/NTUA
Copyright (C) 2013-2020 FOSS GP
Copyright (C) 2007-2021 PCOpt/NTUA
Copyright (C) 2013-2021 FOSS GP
Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
......@@ -31,6 +31,11 @@ License
#include "PrimitivePatchInterpolation.H"
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
#include "faCFD.H"
#include "fixedValueFaPatchFieldsFwd.H"
#include "fixedValueFaPatchFields.H"
#include "zeroGradientFaPatchFieldsFwd.H"
#include "zeroGradientFaPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -42,7 +47,7 @@ namespace incompressible
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defineTypeNameAndDebug(sensitivitySurface, 0);
defineTypeNameAndDebug(sensitivitySurface, 1);
addToRunTimeSelectionTable
(
adjointSensitivity,
......@@ -214,6 +219,188 @@ void sensitivitySurface::setSuffixName()
}
void sensitivitySurface::smoothSensitivities()
{
// Read in parameters
const label iters(dict().getOrDefault<label>("iters", 500));
const scalar tolerance(dict().getOrDefault<scalar>("tolerance", 1.e-06));
autoPtr<faMesh> aMeshPtr(nullptr);
IOobject faceLabels
(
"faceLabels",
mesh_.time().findInstance
(
mesh_.dbDir()/faMesh::meshSubDir,
"faceLabels",
IOobject::READ_IF_PRESENT
),
faMesh::meshSubDir,
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
);
// If the faMesh already exists, read it
if (faceLabels.typeHeaderOk<labelIOList>(false))
{
Info<< "Reading the already constructed faMesh" << endl;
aMeshPtr.reset(new faMesh(mesh_));
}
else
{
// Dictionary used to construct the faMesh
dictionary faMeshDefinition;
IOobject faMeshDefinitionDict
(
"faMeshDefinition",
mesh_.time().caseSystem(),
mesh_,
IOobject::MUST_READ,
IOobject::NO_WRITE
);
// If the faMeshDefinitionDict exists, use it to construct the mesh
if (faMeshDefinitionDict.typeHeaderOk<IOdictionary>(false))
{
Info<< "Reading faMeshDefinition from system " << endl;
faMeshDefinition = IOdictionary(faMeshDefinitionDict);
}
// Otherwise, faMesh is generated from all patches on which we compute
// sensitivities
else
{
Info<< "Constructing faMeshDefinition from sensitivity patches"
<< endl;
wordList polyMeshPatches(sensitivityPatchIDs_.size());
label i(0);
for (const label patchID : sensitivityPatchIDs_)
{
polyMeshPatches[i++] = mesh_.boundary()[patchID].name();
}
faMeshDefinition.add<wordList>("polyMeshPatches", polyMeshPatches);
dictionary& boundary = faMeshDefinition.subDictOrAdd("boundary");
}
// Construct faMesh
aMeshPtr.reset(new faMesh(mesh_, faMeshDefinition));
}
faMesh& aMesh = aMeshPtr.ref();
// Physical radius of the smoothing, provided either directly or computed
// based on the average 'length' of boundary faces
const scalar Rphysical
(dict().getOrDefault<scalar>("radius", computeRadius(aMesh)));
DebugInfo
<< "Physical radius of the sensitivity smoothing "
<< Rphysical << nl << endl;
// Radius used as the diffusivity in the Helmholtz filter, computed as a
// function of the physical radius
const dimensionedScalar RpdeSqr
(
"RpdeSqr", dimArea, sqr(Rphysical/(2.*::sqrt(3.)))
);
dimensionedScalar one("1", dimless, 1.);
// Mapping engine
volSurfaceMapping vsm(aMesh);
// Source term in faMatrix needs to be an areaField
areaScalarField sens
(
IOobject
(
"sens",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
aMesh,
dimensionedScalar(dimless, Zero),
zeroGradientFaPatchField<scalar>::typeName
);
// Copy sensitivities to area field
sens.primitiveFieldRef() =
vsm.mapToSurface<scalar>(wallFaceSensNormalPtr_());
// Initialisation of the smoothed sensitivities field based on the original
// sensitivities
areaScalarField smoothedSens("smoothedSens", sens);
for (label iter = 0; iter < iters; ++iter)
{
Info<< "Sensitivity smoothing iteration " << iter << endl;
faScalarMatrix smoothEqn
(
fam::Sp(one, smoothedSens)
- fam::laplacian(RpdeSqr, smoothedSens)
==
sens
);
smoothEqn.relax();
const scalar residual(mag(smoothEqn.solve().initialResidual()));
DebugInfo
<< "Max smoothSens " << gMax(mag(smoothedSens)()) << endl;
// Print execution time
mesh_.time().printExecutionTime(Info);
// Check convergence
if (residual < tolerance)
{
Info<< "\n***Reached smoothing equation convergence limit, "
"iteration " << iter << "***\n\n";
break;
}
}
// Field used to write the smoothed sensitivity field to file
volScalarField volSmoothedSens
(
IOobject
(
"smoothedSurfaceSens" + surfaceFieldSuffix_,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar(dimless, Zero)
);
// Transfer result back to volField and write
vsm.mapToVolume(smoothedSens, volSmoothedSens.boundaryFieldRef());
volSmoothedSens.write();
}
scalar sensitivitySurface::computeRadius(const faMesh& aMesh)
{
scalar averageArea(gAverage(aMesh.S().field()));
const Vector<label>& geometricD = mesh_.geometricD();
const boundBox& bounds = mesh_.bounds();
forAll(geometricD, iDir)
{
if (geometricD[iDir] == -1)
{
averageArea /= bounds.span()[iDir];
}
}
scalar mult = dict().getOrDefault<scalar>("meanRadiusMultiplier", 10);
return mult*pow(averageArea, scalar(1)/scalar(mesh_.nGeometricD() - 1));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
sensitivitySurface::sensitivitySurface
......@@ -244,6 +431,7 @@ sensitivitySurface::sensitivitySurface
includeMeshMovement_(false),
includeObjective_(false),
writeGeometricInfo_(false),
smoothSensitivities_(false),
eikonalSolver_(nullptr),
meshMovementSolver_(nullptr),
......@@ -346,6 +534,8 @@ void sensitivitySurface::read()
dict().getOrDefault<bool>("includeObjectiveContribution", true);
writeGeometricInfo_ =
dict().getOrDefault<bool>("writeGeometricInfo", false);
smoothSensitivities_ =
dict().getOrDefault<bool>("smoothSensitivities", false);
// Allocate new solvers if necessary
if (includeDistance_ && !eikonalSolver_)
......@@ -700,6 +890,12 @@ void sensitivitySurface::assembleSensitivities()
}
nPassedFaces += patch.size();
}
// Smooth sensitivities if needed
if (smoothSensitivities_)
{
smoothSensitivities();
}
}
......
......@@ -5,8 +5,8 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2007-2020 PCOpt/NTUA
Copyright (C) 2013-2020 FOSS GP
Copyright (C) 2007-2021 PCOpt/NTUA
Copyright (C) 2013-2021 FOSS GP
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
......@@ -31,6 +31,27 @@ Class
Description
Calculation of adjoint based sensitivities at wall faces
Reference:
\verbatim
The computation of the sensitivity derivatives follows the (E-)SI
formulation of
Kavvadias, I. S., Papoutsis-Kiachagias, E. M.,
Giannakoglou, K. C. (2015).
On the proper treatment of grid sensitivities in continuous adjoint
methods for shape optimization.
Journal of computational physics, 301, 1-18.
https://doi.org/10.1016/j.jcp.2015.08.012
whereas their smoothing based on the computation of the 'Sobolev
gradient' is derived from
Vassberg J. C., Jameson A. (2006).
Aerodynamic Shape Optimization Part I: Theoretical Background.
VKI Lecture Series, Introduction to Optimization and
Multidisciplinary Design, Brussels, Belgium, 8 March, 2006.
\endverbatim
SourceFiles
sensitivitySurface.C
......@@ -44,6 +65,7 @@ SourceFiles
#include "adjointEikonalSolverIncompressible.H"
#include "adjointMeshMovementSolverIncompressible.H"
#include "deltaBoundary.H"
#include "faMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -97,6 +119,10 @@ protected:
//- Write geometric info for use by external programs
bool writeGeometricInfo_;
//- Smooth sensitivity derivatives based on the computation of the
//- 'Sobolev gradient'
bool smoothSensitivities_;
autoPtr<adjointEikonalSolver> eikonalSolver_;
autoPtr<adjointMeshMovementSolver> meshMovementSolver_;
......@@ -116,6 +142,14 @@ protected:
//- Set suffix name for sensitivity fields
void setSuffixName();
//- Smooth sensitivity derivatives based on the computation of the
//- 'Sobolev gradient'
void smoothSensitivities();
//- Compute the physical smoothing radius based on the average boundary
//- face 'length'
scalar computeRadius(const faMesh& aMesh);
private:
......@@ -177,7 +211,7 @@ public:
//- Write sensitivity maps
virtual void write(const word& baseName = word::null);
// Inline geters and setters
// Inline getters and setters
//- Get access to the includeObjective bool
inline bool getIncludeObjective() const;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment