ENH: refactoring of the sensitivity classes

Before the commit, the sensitivity classes were receiving references of
the (incompressible) primal and adjoint variables.  However, if
additional physics was added (energy equation, multiphase, etc), the
infrastructure wasn't convenient for accommodating (new terms in the FI
and E-SI formulations, new terms in the sensitivity map, etc).

Now, the sensitivity classes receive a reference to an
incompressibleAdjointSolver and receive the terms for the FI and
sensitivity maps through there. The latter is still WIP.

Modified adjointSimple to incorporate these changes as well.

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C

@@ -79,19 +79,10 @@ FIBase::FIBase
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    shapeSensitivities
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    shapeSensitivities(mesh, dict, adjointSolver),
     gradDxDbMult_
     (
         IOobject

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.H

@@ -109,9 +109,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.C

@@ -79,19 +79,10 @@ SIBase::SIBase
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    shapeSensitivities
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    shapeSensitivities(mesh, dict, adjointSolver),
     surfaceSensitivity_
     (
         mesh,
@@ -100,9 +91,7 @@ SIBase::SIBase
         // and the dict returned by subOrEmptyDict (if found)
         // does not know its parent, optionalSubDict is used
         dict.optionalSubDict("surfaceSensitivities"),
-        primalVars,
-        adjointVars,
-        objectiveManager
+        adjointSolver
     ),
     includeObjective_(true),
     writeSensitivityMap_(true)

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.H

@@ -106,9 +106,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointSensitivity/adjointSensitivityIncompressible.C

@@ -53,16 +53,15 @@ adjointSensitivity::adjointSensitivity
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
     sensitivity(mesh, dict),
     derivatives_(0),
-    primalVars_(primalVars),
-    adjointVars_(adjointVars),
-    objectiveManager_(objectiveManager)
+    adjointSolver_(adjointSolver),
+    primalVars_(adjointSolver.getPrimalVars()),
+    adjointVars_(adjointSolver.getAdjointVars()),
+    objectiveManager_(adjointSolver.getObjectiveManager())
 {}
 
 
@@ -72,9 +71,7 @@ autoPtr<adjointSensitivity> adjointSensitivity::New
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 {
     const word modelType(dict.get<word>("type"));
@@ -96,14 +93,7 @@ autoPtr<adjointSensitivity> adjointSensitivity::New
 
     return autoPtr<adjointSensitivity>
     (
-        ctorPtr
-        (
-            mesh,
-            dict,
-            primalVars,
-            adjointVars,
-            objectiveManager
-        )
+        ctorPtr(mesh, dict, adjointSolver)
     );
 }
 
@@ -142,165 +132,7 @@ void adjointSensitivity::write(const word& baseName)
 
 tmp<volTensorField> adjointSensitivity::computeGradDxDbMultiplier()
 {
-    // Term depending on the adjoint turbulence model
-    autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS
-    (
-        adjointVars_.adjointTurbulence()
-    );
-    tmp<volTensorField> tturbulenceTerm(adjointRAS->FISensitivityTerm());
-    volTensorField& turbulenceTerm = tturbulenceTerm.ref();
-
-    // nu effective
-    tmp<volScalarField> tnuEff(adjointRAS->nuEff());
-    const volScalarField& nuEff = tnuEff();
-
-    tmp<volTensorField> tflowTerm
-    (
-        new volTensorField
-        (
-            IOobject
-            (
-               "flowTerm",
-               mesh_.time().timeName(),
-               mesh_,
-               IOobject::NO_READ,
-               IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedTensor(sqr(dimLength)/pow3(dimTime), Zero)
-        )
-    );
-    volTensorField& flowTerm = tflowTerm.ref();
-
-    const volScalarField& p = primalVars_.p();
-    const volVectorField& U = primalVars_.U();
-    const volScalarField& pa = adjointVars_.pa();
-    const volVectorField& Ua = adjointVars_.Ua();
-    volTensorField gradU(fvc::grad(U));
-    volTensorField gradUa(fvc::grad(Ua));
-
-    // Explicitly correct the boundary gradient to get rid of
-    // the tangential component
-    forAll(mesh_.boundary(), patchI)
-    {
-        const fvPatch& patch = mesh_.boundary()[patchI];
-        if (isA<wallFvPatch>(patch))
-        {
-            tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
-            const vectorField& nf = tnf();
-            gradU.boundaryFieldRef()[patchI] =
-                nf*U.boundaryField()[patchI].snGrad();
-            //gradUa.boundaryField()[patchI] =
-            //    nf*Ua.boundaryField()[patchI].snGrad();
-        }
-    }
-
-    volTensorField stress(nuEff*(gradU + T(gradU)));
-    autoPtr<volVectorField> stressXPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressX", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressYPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressY", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressZPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressZ", stress.dimensions())
-    );
-
-    stressXPtr().replace(0, stress.component(0));
-    stressXPtr().replace(1, stress.component(1));
-    stressXPtr().replace(2, stress.component(2));
-
-    stressYPtr().replace(0, stress.component(3));
-    stressYPtr().replace(1, stress.component(4));
-    stressYPtr().replace(2, stress.component(5));
-
-    stressZPtr().replace(0, stress.component(6));
-    stressZPtr().replace(1, stress.component(7));
-    stressZPtr().replace(2, stress.component(8));
-
-    volTensorField gradStressX(fvc::grad(stressXPtr()));
-    volTensorField gradStressY(fvc::grad(stressYPtr()));
-    volTensorField gradStressZ(fvc::grad(stressZPtr()));
-
-    // Contribution from objective functions and constraints
-    volTensorField objectiveContributions
-    (
-        IOobject
-        (
-            "objectiveContributions",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedTensor(sqr(dimLength)/pow3(dimTime), Zero)
-    );
-    PtrList<objective>& functions(objectiveManager_.getObjectiveFunctions());
-    forAll(functions, funcI)
-    {
-        objectiveContributions +=
-            functions[funcI].weight()
-           *functions[funcI].gradDxDbMultiplier();
-    }
-
-    // Note:
-    // term4 (Ua & grad(stress)) is numerically tricky.  Its div leads to third
-    // order spatial derivs in E-SI based computations Applying the product
-    // derivative rule (putting Ua inside the grad) gives better results in
-    // NACA0012, SA, WF.  However, the original formulation should be kept at
-    // the boundary in order to respect the Ua boundary conditions (necessary
-    // for E-SI to give the same sens as FI).  A mixed approach is hence
-    // followed
-    volTensorField term4
-    (
-      - nuEff*(gradUa & (gradU + T(gradU)))
-      + fvc::grad(nuEff * Ua & (gradU + T(gradU)))
-    );
-
-    forAll(mesh_.boundary(), pI)
-    {
-        if (!isA<coupledFvPatch>(mesh_.boundary()[pI]))
-        {
-            term4.boundaryFieldRef()[pI] =
-                Ua.component(0)().boundaryField()[pI]
-               *gradStressX.boundaryField()[pI]
-              + Ua.component(1)().boundaryField()[pI]
-               *gradStressY.boundaryField()[pI]
-              + Ua.component(2)().boundaryField()[pI]
-               *gradStressZ.boundaryField()[pI];
-        }
-    }
-
-    const autoPtr<ATCModel>& ATCModel =
-        mesh_.lookupObject<incompressibleAdjointSolver>
-        (
-            objectiveManager_.adjointSolverName()
-        ).getATCModel();
-
-    // Compute dxdb multiplier
-    flowTerm =
-        // Term 1, ATC
-        ATCModel->getFISensitivityTerm()
-        // Term 2
-      - fvc::grad(p) * Ua
-        // Term 3
-      - nuEff*(gradU & (gradUa + T(gradUa)))
-        // Term 4
-      + term4
-        // Term 5
-      + (pa * gradU)
-        // Term 6, from the adjoint turbulence model
-      + turbulenceTerm.T()
-        // Term 7, term from objective functions
-      + objectiveContributions;
-
-    flowTerm.correctBoundaryConditions();
-
-    return (tflowTerm);
+    return adjointSolver_.computeGradDxDbMultiplier();
 }
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointSensitivity/adjointSensitivityIncompressible.H

@@ -66,6 +66,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration
+class incompressibleAdjointSolver;
+
 namespace incompressible
 {
 
@@ -82,7 +85,8 @@ protected:
     // Protected data
 
         scalarField derivatives_;
-        incompressibleVars& primalVars_;
+        incompressibleAdjointSolver& adjointSolver_;
+        const incompressibleVars& primalVars_;
         incompressibleAdjointVars& adjointVars_;
         objectiveManager& objectiveManager_;
 
@@ -114,16 +118,12 @@ public:
             (
                 const fvMesh& mesh,
                 const dictionary& dict,
-                incompressibleVars& primalVars,
-                incompressibleAdjointVars& adjointVars,
-                objectiveManager& objectiveManager
+                incompressibleAdjointSolver& adjointSolver
             ),
             (
                 mesh,
                 dict,
-                primalVars,
-                adjointVars,
-                objectiveManager
+                adjointSolver
             )
         );
 
@@ -135,9 +135,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
     // Selectors
@@ -147,9 +145,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.C

@@ -55,19 +55,10 @@ sensitivityBezier::sensitivityBezier
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    SIBase
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    SIBase(mesh, dict, adjointSolver),
     //Bezier_(mesh, dict), // AJH Read locally?
     Bezier_(mesh, mesh.lookupObject<IOdictionary>("optimisationDict")),
     sens_(Bezier_.nBezier(), Zero),

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.H

@@ -107,9 +107,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.C

@@ -122,19 +122,10 @@ sensitivityBezierFI::sensitivityBezierFI
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    FIBase
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    FIBase(mesh, dict, adjointSolver),
     //Bezier_(mesh, dict), // AJH Read locally?
     Bezier_(mesh, mesh.lookupObject<IOdictionary>("optimisationDict")),
     flowSens_(3*Bezier_.nBezier(), Zero),

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.H

@@ -135,9 +135,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityMultiple/sensitivityMultipleIncompressible.C

@@ -54,19 +54,10 @@ sensitivityMultiple::sensitivityMultiple
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    adjointSensitivity
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    adjointSensitivity(mesh, dict, adjointSolver),
     sensTypes_(dict.subDict("sensTypes").toc()),
     sens_(sensTypes_.size())
 {
@@ -79,9 +70,7 @@ sensitivityMultiple::sensitivityMultiple
             (
                 mesh,
                 dict.subDict("sensTypes").subDict(sensTypes_[sI]),
-                primalVars,
-                adjointVars,
-                objectiveManager
+                adjointSolver
             )
         );
     }

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityMultiple/sensitivityMultipleIncompressible.H

@@ -90,9 +90,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C

@@ -28,6 +28,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sensitivitySurfaceIncompressible.H"
+#include "incompressibleAdjointSolver.H"
 #include "PrimitivePatchInterpolation.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
@@ -405,19 +406,10 @@ sensitivitySurface::sensitivitySurface
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    adjointSensitivity
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    adjointSensitivity(mesh, dict, adjointSolver),
     shapeSensitivitiesBase(mesh, dict),
     includeSurfaceArea_(false),
     includePressureTerm_(false),
@@ -798,6 +790,10 @@ void sensitivitySurface::accumulateIntegrand(const scalar dt)
         )*dt;
     }
 
+    // Add terms from physics other than the typical incompressible flow eqns
+    adjointSolver_.additionalSensitivityMapTerms
+        (wallFaceSensVecPtr_(), sensitivityPatchIDs_, dt);
+
     // Add the sensitivity part corresponding to changes of the normal vector
     // Computed at points and mapped to faces
     addGeometricSens();

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H

@@ -175,9 +175,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C

@@ -28,6 +28,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sensitivitySurfacePointsIncompressible.H"
+#include "incompressibleAdjointSolver.H"
 #include "addToRunTimeSelectionTable.H"
 #include "syncTools.H"
 
@@ -333,19 +334,10 @@ sensitivitySurfacePoints::sensitivitySurfacePoints
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    adjointSensitivity
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    adjointSensitivity(mesh, dict, adjointSolver),
     shapeSensitivitiesBase(mesh, dict),
     includeSurfaceArea_(false),
     includePressureTerm_(false),
@@ -611,6 +603,10 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt)
           + dxdbMultiplierTot
         )*dt;
     }
+
+    // Add terms from physics other than the typical incompressible flow eqns
+    adjointSolver_.additionalSensitivityMapTerms
+        (wallFaceSens_(), sensitivityPatchIDs_, dt);
 }
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.H

@@ -159,9 +159,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.C

@@ -143,19 +143,10 @@ sensitivityVolBSplines::sensitivityVolBSplines
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    SIBase
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    SIBase(mesh, dict, adjointSolver),
     volBSplinesBase_
     (
         const_cast<volBSplinesBase&>(volBSplinesBase::New(mesh))

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.H

@@ -115,9 +115,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C

@@ -56,19 +56,10 @@ sensitivityVolBSplinesFI::sensitivityVolBSplinesFI
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    FIBase
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    FIBase(mesh, dict, adjointSolver),
     volBSplinesBase_
     (
         const_cast<volBSplinesBase&>(volBSplinesBase::New(mesh))

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.H

@@ -118,9 +118,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C

@@ -132,19 +132,10 @@ shapeSensitivities::shapeSensitivities
 (
     const fvMesh& mesh,
     const dictionary& dict,
-    incompressibleVars& primalVars,
-    incompressibleAdjointVars& adjointVars,
-    objectiveManager& objectiveManager
+    incompressibleAdjointSolver& adjointSolver
 )
 :
-    adjointSensitivity
-    (
-        mesh,
-        dict,
-        primalVars,
-        adjointVars,
-        objectiveManager
-    ),
+    adjointSensitivity(mesh, dict, adjointSolver),
     shapeSensitivitiesBase(mesh, dict),
     dSfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),
     dnfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.H

@@ -105,9 +105,7 @@ public:
         (
             const fvMesh& mesh,
             const dictionary& dict,
-            incompressibleVars& primalVars,
-            incompressibleAdjointVars& adjointVars,
-            objectiveManager& objectiveManager
+            incompressibleAdjointSolver& adjointSolver
         );
 
 

src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C

@@ -146,9 +146,7 @@ Foam::adjointSimple::adjointSimple
             (
                 mesh,
                 optDict.subDict("optimisation").subDict("sensitivities"),
-                primalVars_,
-                adjointVars_,
-                objectiveManagerPtr_()
+                *this
             ).ptr()
         );
         // Read stored sensitivities, if they exist

src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.C

@@ -29,6 +29,7 @@ License
 
 #include "incompressibleAdjointSolver.H"
 #include "incompressiblePrimalSolver.H"
+#include "wallFvPatch.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -166,4 +167,175 @@ void Foam::incompressibleAdjointSolver::updatePrimalBasedQuantities()
 }
 
 
+Foam::tmp<Foam::volTensorField>
+Foam::incompressibleAdjointSolver::computeGradDxDbMultiplier()
+{
+    // Term depending on the adjoint turbulence model
+    autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS
+    (
+        getAdjointVars().adjointTurbulence()
+    );
+    tmp<volTensorField> tturbulenceTerm(adjointRAS->FISensitivityTerm());
+    volTensorField& turbulenceTerm = tturbulenceTerm.ref();
+
+    // nu effective
+    tmp<volScalarField> tnuEff(adjointRAS->nuEff());
+    const volScalarField& nuEff = tnuEff();
+
+    tmp<volTensorField> tflowTerm
+    (
+        new volTensorField
+        (
+            IOobject
+            (
+               "flowTerm",
+               mesh_.time().timeName(),
+               mesh_,
+               IOobject::NO_READ,
+               IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedTensor(sqr(dimLength)/pow3(dimTime), Zero)
+        )
+    );
+    volTensorField& flowTerm = tflowTerm.ref();
+
+    const volScalarField& p = primalVars_.p();
+    const volVectorField& U = primalVars_.U();
+    const volScalarField& pa = getAdjointVars().pa();
+    const volVectorField& Ua = getAdjointVars().Ua();
+    volTensorField gradU(fvc::grad(U));
+    volTensorField gradUa(fvc::grad(Ua));
+
+    // Explicitly correct the boundary gradient to get rid of
+    // the tangential component
+    forAll(mesh_.boundary(), patchI)
+    {
+        const fvPatch& patch = mesh_.boundary()[patchI];
+        if (isA<wallFvPatch>(patch))
+        {
+            tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
+            const vectorField& nf = tnf();
+            gradU.boundaryFieldRef()[patchI] =
+                nf*U.boundaryField()[patchI].snGrad();
+            //gradUa.boundaryField()[patchI] =
+            //    nf*Ua.boundaryField()[patchI].snGrad();
+        }
+    }
+
+    volTensorField stress(nuEff*(gradU + T(gradU)));
+    autoPtr<volVectorField> stressXPtr
+    (
+        createZeroFieldPtr<vector>(mesh_, "stressX", stress.dimensions())
+    );
+    autoPtr<volVectorField> stressYPtr
+    (
+        createZeroFieldPtr<vector>(mesh_, "stressY", stress.dimensions())
+    );
+    autoPtr<volVectorField> stressZPtr
+    (
+        createZeroFieldPtr<vector>(mesh_, "stressZ", stress.dimensions())
+    );
+
+    stressXPtr().replace(0, stress.component(0));
+    stressXPtr().replace(1, stress.component(1));
+    stressXPtr().replace(2, stress.component(2));
+
+    stressYPtr().replace(0, stress.component(3));
+    stressYPtr().replace(1, stress.component(4));
+    stressYPtr().replace(2, stress.component(5));
+
+    stressZPtr().replace(0, stress.component(6));
+    stressZPtr().replace(1, stress.component(7));
+    stressZPtr().replace(2, stress.component(8));
+
+    volTensorField gradStressX(fvc::grad(stressXPtr()));
+    volTensorField gradStressY(fvc::grad(stressYPtr()));
+    volTensorField gradStressZ(fvc::grad(stressZPtr()));
+
+    // Contribution from objective functions and constraints
+    volTensorField objectiveContributions
+    (
+        IOobject
+        (
+            "objectiveContributions",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedTensor(sqr(dimLength)/pow3(dimTime), Zero)
+    );
+    PtrList<objective>& functions
+        (objectiveManagerPtr_->getObjectiveFunctions());
+    forAll(functions, funcI)
+    {
+        objectiveContributions +=
+            functions[funcI].weight()
+           *functions[funcI].gradDxDbMultiplier();
+    }
+
+    // Note:
+    // term4 (Ua & grad(stress)) is numerically tricky.  Its div leads to third
+    // order spatial derivs in E-SI based computations Applying the product
+    // derivative rule (putting Ua inside the grad) gives better results in
+    // NACA0012, SA, WF.  However, the original formulation should be kept at
+    // the boundary in order to respect the Ua boundary conditions (necessary
+    // for E-SI to give the same sens as FI).  A mixed approach is hence
+    // followed
+    volTensorField term4
+    (
+      - nuEff*(gradUa & (gradU + T(gradU)))
+      + fvc::grad(nuEff * Ua & (gradU + T(gradU)))
+    );
+
+    forAll(mesh_.boundary(), pI)
+    {
+        if (!isA<coupledFvPatch>(mesh_.boundary()[pI]))
+        {
+            term4.boundaryFieldRef()[pI] =
+                Ua.component(0)().boundaryField()[pI]
+               *gradStressX.boundaryField()[pI]
+              + Ua.component(1)().boundaryField()[pI]
+               *gradStressY.boundaryField()[pI]
+              + Ua.component(2)().boundaryField()[pI]
+               *gradStressZ.boundaryField()[pI];
+        }
+    }
+
+    // Compute dxdb multiplier
+    flowTerm =
+        // Term 1, ATC
+        ATCModel_->getFISensitivityTerm()
+        // Term 2
+      - fvc::grad(p) * Ua
+        // Term 3
+      - nuEff*(gradU & (gradUa + T(gradUa)))
+        // Term 4
+      + term4
+        // Term 5
+      + (pa * gradU)
+        // Term 6, from the adjoint turbulence model
+      + turbulenceTerm.T()
+        // Term 7, term from objective functions
+      + objectiveContributions;
+
+    flowTerm.correctBoundaryConditions();
+
+    return (tflowTerm);
+}
+
+
+void Foam::incompressibleAdjointSolver::additionalSensitivityMapTerms
+(
+    boundaryVectorField& sensitivityMap,
+    const labelHashSet& patchIDs,
+    const scalar dt
+)
+{
+    // Does nothing in base
+}
+
+
 // ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.H

@@ -168,6 +168,26 @@ public:
             virtual void updatePrimalBasedQuantities();
 
 
+        // Sensitivity related functions
+
+            //- Compute the multiplier for grad(dxdb)
+            //  Used in shape sensitivity derivatives, computed with the FI
+            //  and E-SI approaches
+            virtual tmp<volTensorField> computeGradDxDbMultiplier();
+
+            //- Terms to be added to the sensitivity map, depending
+            //- on the adjoint solver
+            //  The typical terms associated with the typical incompressible
+            //  flow equations are included in the sensitivity classes.
+            //  This is to be used whenever additional physics is added
+            virtual void additionalSensitivityMapTerms
+            (
+                boundaryVectorField& sensitivityMap,
+                const labelHashSet& patchIDs,
+                const scalar dt
+            );
+
+
         // IO
 
             //- In case of multi-point runs with turbulent flows,