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.

src/optimisation/adjointOptimisation/adjoint/Make/options

 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 \

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();
+    }
 }
 
 

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;