ENH: adjointRASModel will not crash if set to laminar

and the Jacobian of the objective function wrt the turbulence variables
is called (rare/unorthodox case).

Additionally, objectivePowerDissipation dissipation can now be used in
topology optimisation, adding the necessary blockage dependency to it.

src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C

@@ -373,7 +373,7 @@ void objectiveIncompressible::update_dJdTMvar
     const labelList& zones
 )
 {
-    if (dJdTMvarPtr)
+    if (dJdTMvarPtr.good())
     {
         // nut Jacobians are currently computed in the adjoint turbulence
         // models, though they would be better placed within the primal

src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveNutSqr/objectiveNutSqr.C

@@ -51,6 +51,8 @@ addToRunTimeSelectionTable
 );
 
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
 void objectiveNutSqr::populateFieldNames()
 {
     if (adjointTurbulenceNames_.empty())

src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.C

@@ -29,6 +29,7 @@ License
 #include "objectivePowerDissipation.H"
 #include "incompressibleAdjointSolver.H"
 #include "createZeroField.H"
+#include "topOVariablesBase.H"
 #include "IOmanip.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -51,6 +52,8 @@ addToRunTimeSelectionTable
 );
 
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
 void objectivePowerDissipation::populateFieldNames()
 {
     if (fieldNames_.size() == 1)
@@ -134,6 +137,9 @@ objectivePowerDissipation::objectivePowerDissipation
             sqr(dimLength)/pow3(dimTime)
         )
     );
+
+    // Allocate direct sensitivity contributions for topology optimisation
+    dJdbPtr_.reset(createZeroFieldPtr<scalar>(mesh_, "dJdb", dimless));
 }
 
 
@@ -157,6 +163,19 @@ scalar objectivePowerDissipation::J()
         scalarField integrandZone(integrand.primitiveField(), zoneI);
 
         J_ += 0.5*gSum(integrandZone*VZone);
+        if (mesh_.foundObject<topOVariablesBase>("topoVars"))
+        {
+            const topOVariablesBase& vars =
+                mesh_.lookupObject<topOVariablesBase>("topoVars");
+            const volScalarField& beta = vars.beta();
+            scalar porosityContr = Zero;
+            for (const label cellI : zoneI)
+            {
+                porosityContr += beta[cellI]*magSqr(U[cellI])*V[cellI];
+            }
+            porosityContr *= vars.getBetaMax();
+            J_ += returnReduce(porosityContr, sumOp<scalar>());
+        }
     }
 
     return J_;
@@ -204,6 +223,23 @@ void objectivePowerDissipation::update_dJdv()
             dJdvPtr_()[cellI] += integrand[cellI];
         }
     }
+
+    // Add source from porosity dependencies
+    if (mesh_.foundObject<topOVariablesBase>("topoVars"))
+    {
+        const topOVariablesBase& vars =
+            mesh_.lookupObject<topOVariablesBase>("topoVars");
+        const volScalarField& beta = vars.beta();
+        const scalar betaMax = vars.getBetaMax();
+        for (const label zI : zones_)
+        {
+            const cellZone& zoneI = mesh_.cellZones()[zI];
+            for (const label cellI : zoneI)
+            {
+                dJdvPtr_()[cellI] += 2*betaMax*beta[cellI]*U[cellI];
+            }
+        }
+    }
 }
 
 
@@ -276,6 +312,30 @@ void objectivePowerDissipation::update_gradDxDbMultiplier()
 }
 
 
+void objectivePowerDissipation::update_dJdb()
+{
+    if (mesh_.foundObject<topOVariablesBase>("topoVars"))
+    {
+        scalarField& dJdb = dJdbPtr_().primitiveFieldRef();
+        dJdb = Zero;
+        const volVectorField& U = vars_.UInst();
+        const topOVariablesBase& vars =
+            mesh_.lookupObject<topOVariablesBase>("topoVars");
+        const scalar betaMax = vars.getBetaMax();
+        for (const label zI : zones_)
+        {
+            const cellZone& zoneI = mesh_.cellZones()[zI];
+            for (const label cellI : zoneI)
+            {
+                dJdb[cellI] += betaMax*magSqr(U[cellI]);
+            }
+        }
+    }
+}
+
+
+
+
 void objectivePowerDissipation::addSource(fvScalarMatrix& matrix)
 {
     populateFieldNames();

src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.H

@@ -114,6 +114,9 @@ public:
         //- Update grad(dx/db multiplier). Volume-based sensitivity term
         virtual void update_gradDxDbMultiplier();
 
+        //- Contribution to field sensitivities
+        virtual void update_dJdb();
+
         //- Add source terms to the adjoint turbulence model equations
         virtual void addSource(fvScalarMatrix& matrix);
 };

src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C

@@ -380,6 +380,13 @@ const wordList& adjointRASModel::getAdjointTMVariablesBaseNames() const
 
 tmp<volScalarField> adjointRASModel::nutJacobianTMVar1() const
 {
+    dimensionSet dims
+    (
+        bool(adjointTMVariable1Ptr_)
+      ? dimViscosity/adjointTMVariable1Ptr_().dimensions()
+      : dimless
+    );
+
     return
         tmp<volScalarField>::New
         (
@@ -392,17 +399,20 @@ tmp<volScalarField> adjointRASModel::nutJacobianTMVar1() const
                 IOobject::NO_WRITE
             ),
             mesh_,
-            dimensionedScalar
-            (
-               dimViscosity/adjointTMVariable1Ptr_().dimensions(),
-               Zero
-            )
+            dimensionedScalar(dims, Zero)
         );
 }
 
 
 tmp<volScalarField> adjointRASModel::nutJacobianTMVar2() const
 {
+    dimensionSet dims
+    (
+        bool(adjointTMVariable2Ptr_)
+      ? dimViscosity/adjointTMVariable2Ptr_().dimensions()
+      : dimless
+    );
+
     return
         tmp<volScalarField>::New
         (
@@ -415,11 +425,7 @@ tmp<volScalarField> adjointRASModel::nutJacobianTMVar2() const
                 IOobject::NO_WRITE
             ),
             mesh_,
-            dimensionedScalar
-            (
-                dimViscosity/adjointTMVariable2Ptr_().dimensions(),
-                Zero
-            )
+            dimensionedScalar(dims, Zero)
         );
 }