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 1/3] 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 44dec84580..f3041856b6 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 90bcce43aa..e17c5365a8 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 ec070cb0da..baff1c1ecf 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;

From 90d62c33e432fb1b662756a039de3cfa6ccf9662 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Mon, 29 Nov 2021 15:10:19 +0000
Subject: [PATCH 2/3] INT: Minor integration updates

---
 .../sensitivitySurfaceIncompressible.C        | 48 +++++++++----------
 1 file changed, 23 insertions(+), 25 deletions(-)

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 e17c5365a8..7c5cafa6fd 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
@@ -31,10 +31,11 @@ License
 #include "PrimitivePatchInterpolation.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
-#include "faCFD.H"
-#include "fixedValueFaPatchFieldsFwd.H"
+#include "faMatrices.H"
+#include "famSup.H"
+#include "famLaplacian.H"
+#include "volSurfaceMapping.H"
 #include "fixedValueFaPatchFields.H"
-#include "zeroGradientFaPatchFieldsFwd.H"
 #include "zeroGradientFaPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -79,7 +80,7 @@ void sensitivitySurface::addGeometricSens()
         for (const label patchI : sensitivityPatchIDs_)
         {
             const fvPatch& patch = mesh_.boundary()[patchI];
-            vectorField nf(patch.nf());
+            const vectorField nf(patch.nf());
 
             // point sens result for patch
             vectorField& patchdSdb = pointSensdSdb()[patchI];
@@ -87,14 +88,10 @@ void sensitivitySurface::addGeometricSens()
 
             vectorField dSdbMultiplierTot(patch.size(), Zero);
             vectorField dndbMultiplierTot(patch.size(), Zero);
-            forAll(functions, funcI)
+            for (auto& fun : functions)
             {
-                dSdbMultiplierTot +=
-                    functions[funcI].weight()
-                   *functions[funcI].dSdbMultiplier(patchI);
-                dndbMultiplierTot +=
-                    functions[funcI].weight()
-                   *functions[funcI].dndbMultiplier(patchI);
+                dSdbMultiplierTot += fun.weight()*fun.dSdbMultiplier(patchI);
+                dndbMultiplierTot += fun.weight()*fun.dndbMultiplier(patchI);
             }
             // Correspondence of local point addressing to global point
             // addressing
@@ -114,13 +111,12 @@ void sensitivitySurface::addGeometricSens()
             forAll(meshPoints, ppI)
             {
                 const labelList& pointFaces = patchPointFaces[ppI];
-                forAll(pointFaces, pfI)
+                for (label localFaceIndex : pointFaces)
                 {
-                    label localFaceIndex = pointFaces[pfI];
                     label globalFaceIndex = patchStartIndex + localFaceIndex;
                     const face& faceI = faces[globalFaceIndex];
                     // Point coordinates. All indices in global numbering
-                    pointField p(faceI.points(mesh_.points()));
+                    const pointField p(faceI.points(mesh_.points()));
                     tensorField p_d(faceI.size(), Zero);
                     forAll(faceI, facePointI)
                     {
@@ -129,8 +125,10 @@ void sensitivitySurface::addGeometricSens()
                             p_d[facePointI] = tensor::I;
                         }
                     }
-                    tensorField deltaNormals =
-                        dBoundary.makeFaceCentresAndAreas_d(p, p_d);
+                    const tensorField deltaNormals
+                    (
+                        dBoundary.makeFaceCentresAndAreas_d(p, p_d)
+                    );
 
                     // Element [1] is the variation in the (dimensional) normal
                     const tensor& deltaSf = deltaNormals[1];
@@ -139,7 +137,7 @@ void sensitivitySurface::addGeometricSens()
 
                     // Element [2] is the variation in the unit normal
                     const tensor& deltaNf = deltaNormals[2];
-                    patchdndb[ppI]   +=
+                    patchdndb[ppI] +=
                         dndbMultiplierTot[localFaceIndex] & deltaNf;
                 }
             }
@@ -231,8 +229,8 @@ void sensitivitySurface::smoothSensitivities()
         "faceLabels",
         mesh_.time().findInstance
         (
-            mesh_.dbDir()/faMesh::meshSubDir, 
-            "faceLabels", 
+            mesh_.dbDir()/faMesh::meshSubDir,
+            "faceLabels",
             IOobject::READ_IF_PRESENT
         ),
         faMesh::meshSubDir,
@@ -260,7 +258,7 @@ void sensitivitySurface::smoothSensitivities()
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         );
-        
+
         // If the faMeshDefinitionDict exists, use it to construct the mesh
         if (faMeshDefinitionDict.typeHeaderOk<IOdictionary>(false))
         {
@@ -271,16 +269,16 @@ void sensitivitySurface::smoothSensitivities()
         // sensitivities
         else
         {
-            Info<< "Constructing faMeshDefinition from sensitivity patches" 
+            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();    
+                polyMeshPatches[i++] = mesh_.boundary()[patchID].name();
             }
             faMeshDefinition.add<wordList>("polyMeshPatches", polyMeshPatches);
-            dictionary& boundary = faMeshDefinition.subDictOrAdd("boundary");
+            (void)faMeshDefinition.subDictOrAdd("boundary");
         }
 
         // Construct faMesh
@@ -293,10 +291,10 @@ void sensitivitySurface::smoothSensitivities()
     const scalar Rphysical
         (dict().getOrDefault<scalar>("radius", computeRadius(aMesh)));
     DebugInfo
-        << "Physical radius of the sensitivity smoothing " 
+        << "Physical radius of the sensitivity smoothing "
         << Rphysical << nl << endl;
 
-    // Radius used as the diffusivity in the Helmholtz filter, computed as a 
+    // Radius used as the diffusivity in the Helmholtz filter, computed as a
     // function of the physical radius
     const dimensionedScalar RpdeSqr
     (

From f6dfa76f6d3737694327209e777288e4f6de1300 Mon Sep 17 00:00:00 2001
From: Vaggelis Papoutsis <vaggelisp@gmail.com>
Date: Wed, 8 Dec 2021 16:05:59 +0200
Subject: [PATCH 3/3] TUT: updated the motorbike tutorial

to also compute smoothed sensitivity maps
---
 .../motorBike/system/faSchemes                | 47 +++++++++++++++
 .../motorBike/system/faSolution               | 33 ++++++++++
 .../motorBike/system/optimisationDict         | 60 ++++++++++++++++---
 3 files changed, 133 insertions(+), 7 deletions(-)
 create mode 100644 tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSchemes
 create mode 100644 tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSolution

diff --git a/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSchemes b/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSchemes
new file mode 100644
index 0000000000..46b178e110
--- /dev/null
+++ b/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSchemes
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2108                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      faSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default            steadyState;
+}
+
+gradSchemes
+{
+    default            Gauss linear;
+}
+
+divSchemes
+{
+    default            Gauss linear;
+}
+
+laplacianSchemes
+{
+    default            Gauss linear limited 0.333;
+}
+
+interpolationSchemes
+{
+    default            linear;
+}
+
+snGradSchemes
+{
+    default            limited 0.333;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSolution b/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSolution
new file mode 100644
index 0000000000..a939d23203
--- /dev/null
+++ b/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/faSolution
@@ -0,0 +1,33 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2108                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      faSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    smoothedSens
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-08;
+        relTol          0.01;
+    }
+}
+
+relaxationFactors
+{
+    smoothedSens    0.4;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/optimisationDict b/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/optimisationDict
index e48a9498d3..94effd916b 100644
--- a/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/optimisationDict
+++ b/tutorials/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike/system/optimisationDict
@@ -108,14 +108,60 @@ optimisation
 {
     sensitivities
     {
-        type               surfacePoints;
-        patches            (motorBikeGroup);
-        includeSurfaceArea false;
-        adjointEikonalSolver
+        // Used to compute a number of variants of the sensitivity map
+        // at once
+        type            multiple;
+        patches         (motorBikeGroup);
+        sensTypes
         {
-            tolerance 1.e-5;
-            iters     1000;
-            epsilon   0.1;
+            pointBased
+            {
+                type               surfacePoints;
+                patches            (motorBikeGroup);
+                includeSurfaceArea false;
+                adjointEikonalSolver
+                {
+                    tolerance 1.e-5;
+                    iters     1000;
+                    epsilon   0.1;
+                }
+            }
+            faceBased-unsmoothed
+            {
+                type                 surface;
+                patches              (motorBikeGroup);
+                includeSurfaceArea   false;
+            }
+            faceBased-RMult_2
+            {
+                type                 surface;
+                patches              (motorBikeGroup);
+                includeSurfaceArea   false;
+                smoothSensitivities  true;
+                meanRadiusMultiplier 2;
+                suffix               Rmult2; // suffix of the sensitivity map output files
+                iters                2000;
+            }
+            faceBased-RMult_5
+            {
+                type                 surface;
+                patches              (motorBikeGroup);
+                includeSurfaceArea   false;
+                smoothSensitivities  true;
+                meanRadiusMultiplier 5;
+                suffix               Rmult5; // suffix of the sensitivity map output files
+                iters                2000;
+            }
+            faceBased-RMult_10
+            {
+                type                 surface;
+                patches              (motorBikeGroup);
+                includeSurfaceArea   false;
+                smoothSensitivities  true;
+                meanRadiusMultiplier 10;
+                suffix               Rmult10; // suffix of the sensitivity map output files
+                iters                2000;
+            }
         }
     }
 }