From 3ef0d19ef42bd210ea17114eb9e0c782ba3ce9e5 Mon Sep 17 00:00:00 2001
From: Vaggelis Papoutsis <vaggelisp@gmail.com>
Date: Thu, 22 Jul 2021 19:04:03 +0300
Subject: [PATCH] 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.
---
.../adjointOptimisation/adjoint/Make/options | 2 +
.../sensitivitySurfaceIncompressible.C | 202 +++++++++++++++++-
.../sensitivitySurfaceIncompressible.H | 40 +++-
3 files changed, 238 insertions(+), 6 deletions(-)
diff --git a/src/optimisation/adjointOptimisation/adjoint/Make/options b/src/optimisation/adjointOptimisation/adjoint/Make/options
index 44dec845809..f3041856b6b 100644
--- a/src/optimisation/adjointOptimisation/adjoint/Make/options
+++ b/src/optimisation/adjointOptimisation/adjoint/Make/options
@@ -1,5 +1,6 @@
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 \
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
index 90bcce43aa6..e17c5365a87 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
@@ -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();
+ }
}
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
index ec070cb0da4..baff1c1ecf7 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
@@ -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;
--
GitLab