ENH: add the infrastructure for computing and utilising

the Jacobian of an objective function, defined at the boundary, wrt nut
and gradU. Also modified the current objectives that include such
contributions

src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.C

@@ -103,6 +103,18 @@ tmp<scalarField> boundaryAdjointContribution::adjointTMVariable2Source()
 }
 
 
+tmp<scalarField> boundaryAdjointContribution::dJdnut()
+{
+    return tmp<scalarField>::New(patch_.size(), Zero);
+}
+
+
+tmp<tensorField> boundaryAdjointContribution::dJdGradU()
+{
+    return tmp<tensorField>::New(patch_.size(), Zero);
+}
+
+
 tmp<scalarField> boundaryAdjointContribution::TMVariable1Diffusion()
 {
     return tmp<scalarField>::New(patch_.size(), Zero);

src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.H

@@ -138,6 +138,8 @@ public:
         virtual tmp<vectorField> normalVelocitySource() = 0;
         virtual tmp<scalarField> adjointTMVariable1Source();
         virtual tmp<scalarField> adjointTMVariable2Source();
+        virtual tmp<scalarField> dJdnut();
+        virtual tmp<tensorField> dJdGradU();
         virtual tmp<scalarField> energySource() = 0;
 
         virtual tmp<scalarField> momentumDiffusion() = 0;

src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C

@@ -223,6 +223,30 @@ boundaryAdjointContributionIncompressible::adjointTMVariable2Source()
 }
 
 
+tmp<scalarField>
+boundaryAdjointContributionIncompressible::dJdnut()
+{
+    return
+        sumContributions
+        (
+            objectiveManager_.getObjectiveFunctions(),
+            &objectiveIncompressible::boundarydJdnut
+        );
+}
+
+
+tmp<tensorField>
+boundaryAdjointContributionIncompressible::dJdGradU()
+{
+    return
+        sumContributions
+        (
+            objectiveManager_.getObjectiveFunctions(),
+            &objectiveIncompressible::boundarydJdGradU
+        );
+}
+
+
 tmp<scalarField> boundaryAdjointContributionIncompressible::momentumDiffusion()
 {
     tmp<scalarField> tnuEff(new scalarField(patch_.size(), Zero));

src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.H

@@ -139,6 +139,8 @@ public:
         tmp<scalarField> energySource();
         tmp<scalarField> adjointTMVariable1Source();
         tmp<scalarField> adjointTMVariable2Source();
+        tmp<scalarField> dJdnut();
+        tmp<tensorField> dJdGradU();
 
         tmp<scalarField> momentumDiffusion();
         tmp<scalarField> laminarDiffusivity();

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

@@ -113,7 +113,8 @@ objectiveForce::objectiveForce
     bdJdpPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
     bdSdbMultPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
     bdxdbMultPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
-    bdJdStressPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
+    bdJdnutPtr_.reset(createZeroBoundaryPtr<scalar>(mesh_));
+    bdJdGradUPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
 }
 
 
@@ -268,14 +269,35 @@ void objectiveForce::update_dxdbMultiplier()
 }
 
 
-void objectiveForce::update_dJdStressMultiplier()
+void objectiveForce::update_boundarydJdnut()
 {
+    const volVectorField& U = vars_.U();
+    volSymmTensorField devGradU(dev(twoSymm(fvc::grad(U))));
+
     for (const label patchI : forcePatches_)
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
         tmp<vectorField> tnf = patch.nf();
         const vectorField& nf = tnf();
-        bdJdStressPtr_()[patchI] = (forceDirection_ * nf)/denom();
+        bdJdnutPtr_()[patchI] =
+           -((devGradU.boundaryField()[patchI] & forceDirection_) & nf)/denom();
+    }
+}
+
+
+void objectiveForce::update_boundarydJdGradU()
+{
+    const autoPtr<incompressible::RASModelVariables>&
+        turbVars = vars_.RASModelVariables();
+    const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
+    volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
+    for (const label patchI : forcePatches_)
+    {
+        const fvPatch& patch = mesh_.boundary()[patchI];
+        const vectorField& Sf = patch.Sf();
+        bdJdGradUPtr_()[patchI] =
+          - nuEff.boundaryField()[patchI]
+           *dev(forceDirection_*Sf + Sf*forceDirection_);
     }
 }
 

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

@@ -111,8 +111,11 @@ public:
         //- Update delta(x)/delta b multiplier
         void update_dxdbMultiplier();
 
-        //- Update dJ/dStress multiplier
+        //- Update dJ/dnut multiplier
-        void update_dJdStressMultiplier();
+        void update_boundarydJdnut();
+
+        //- Update dJ/dGradU multiplier
+        void update_boundarydJdGradU();
 
         //- Return denominator, without density
         virtual scalar denom() const;

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

@@ -84,7 +84,9 @@ objectiveIncompressible::objectiveIncompressible
     bdJdpPtr_(nullptr),
     bdJdTPtr_(nullptr),
     bdJdTMvar1Ptr_(nullptr),
-    bdJdTMvar2Ptr_(nullptr)
+    bdJdTMvar2Ptr_(nullptr),
+    bdJdnutPtr_(nullptr),
+    bdJdGradUPtr_(nullptr)
 {
     weight_ = dict.get<scalar>("weight");
     computeMeanFields_ = vars_.computeMeanFields();
@@ -182,6 +184,14 @@ void objectiveIncompressible::doNormalization()
         {
             bdJdTMvar2Ptr_() *= oneOverNorm;
         }
+        if (hasBoundarydJdnut())
+        {
+            bdJdnutPtr_() *= oneOverNorm;
+        }
+        if (hasBoundarydJdGradU())
+        {
+            bdJdGradUPtr_() *= oneOverNorm;
+        }
 
         // Normalize geometric fields
         objective::doNormalization();
@@ -375,6 +385,32 @@ const fvPatchScalarField& objectiveIncompressible::boundarydJdTMvar2
 }
 
 
+const fvPatchScalarField& objectiveIncompressible::boundarydJdnut
+(
+    const label patchI
+)
+{
+    if (!bdJdnutPtr_)
+    {
+        bdJdnutPtr_.reset(createZeroBoundaryPtr<scalar>(mesh_));
+    }
+    return bdJdnutPtr_()[patchI];
+}
+
+
+const fvPatchTensorField& objectiveIncompressible::boundarydJdGradU
+(
+    const label patchI
+)
+{
+    if (!bdJdGradUPtr_)
+    {
+        bdJdGradUPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
+    }
+    return bdJdGradUPtr_()[patchI];
+}
+
+
 const boundaryVectorField& objectiveIncompressible::boundarydJdv()
 {
     if (!bdJdvPtr_)
@@ -445,6 +481,26 @@ const boundaryScalarField& objectiveIncompressible::boundarydJdTMvar2()
 }
 
 
+const boundaryScalarField& objectiveIncompressible::boundarydJdnut()
+{
+    if (!bdJdnutPtr_)
+    {
+        bdJdnutPtr_.reset(createZeroBoundaryPtr<scalar>(mesh_));
+    }
+    return bdJdnutPtr_();
+}
+
+
+const boundaryTensorField& objectiveIncompressible::boundarydJdGradU()
+{
+    if (!bdJdGradUPtr_)
+    {
+        bdJdGradUPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
+    }
+    return *bdJdGradUPtr_;
+}
+
+
 void objectiveIncompressible::update()
 {
     // Objective function value
@@ -472,13 +528,14 @@ void objectiveIncompressible::update()
     update_boundarydJdT();
     update_boundarydJdTMvar1();
     update_boundarydJdTMvar2();
+    update_boundarydJdnut();
+    update_boundarydJdGradU();
     update_boundarydJdb();
     update_dSdbMultiplier();
     update_dndbMultiplier();
     update_dxdbMultiplier();
     update_dxdbDirectMultiplier();
     update_boundaryEdgeContribution();
-    update_dJdStressMultiplier();
 
     // Divide everything with normalization factor
     doNormalization();
@@ -539,6 +596,14 @@ void objectiveIncompressible::nullify()
         {
             bdJdTMvar2Ptr_() == scalar(0);
         }
+        if (hasBoundarydJdnut())
+        {
+            bdJdnutPtr_() == scalar(0);
+        }
+        if (hasBoundarydJdGradU())
+        {
+            bdJdGradUPtr_() == tensor::zero;
+        }
 
         // Nullify geometric fields and sets nullified_ to true
         objective::nullify();

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

@@ -97,6 +97,12 @@ protected:
         //- Adjoint outlet turbulence model var 2
         autoPtr<boundaryScalarField> bdJdTMvar2Ptr_;
 
+        //- Jacobian wrt to nut
+        autoPtr<boundaryScalarField> bdJdnutPtr_;
+
+        //- Term multiplying gradU variations
+        autoPtr<boundaryTensorField> bdJdGradUPtr_;
+
 
 private:
 
@@ -207,6 +213,12 @@ public:
         //- patch
         const fvPatchScalarField& boundarydJdTMvar2(const label);
 
+        //- Objective partial deriv wrt nut for a specific patch
+        const fvPatchScalarField& boundarydJdnut(const label);
+
+        //- Objective partial deriv wrt stress tensor
+        const fvPatchTensorField& boundarydJdGradU(const label);
+
         //- Objective partial deriv wrt velocity for all patches
         const boundaryVectorField& boundarydJdv();
 
@@ -228,6 +240,12 @@ public:
         //- Objective partial deriv wrt turbulence model var 2 for all patches
         const boundaryScalarField& boundarydJdTMvar2();
 
+        //- Objective partial deriv wrt nut for all patches
+        const boundaryScalarField& boundarydJdnut();
+
+        //- Objective partial deriv wrt gradU
+        const boundaryTensorField& boundarydJdGradU();
+
         //- Update objective function derivatives
         virtual void update();
 
@@ -282,6 +300,12 @@ public:
         virtual void update_boundarydJdTMvar2()
         {}
 
+        virtual void update_boundarydJdnut()
+        {}
+
+        virtual void update_boundarydJdGradU()
+        {}
+
         virtual void update_boundarydJdb()
         {}
 
@@ -321,6 +345,8 @@ public:
         inline bool hasBoundarydJdT() const;
         inline bool hasBoundarydJdTMVar1() const;
         inline bool hasBoundarydJdTMVar2() const;
+        inline bool hasBoundarydJdnut() const;
+        inline bool hasBoundarydJdGradU() const;
 };
 
 

src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressibleI.H

@@ -102,4 +102,16 @@ inline bool Foam::objectiveIncompressible::hasBoundarydJdTMVar2() const
 }
 
 
+inline bool Foam::objectiveIncompressible::hasBoundarydJdnut() const
+{
+    return bool(bdJdnutPtr_);
+}
+
+
+inline bool Foam::objectiveIncompressible::hasBoundarydJdGradU() const
+{
+    return bool(bdJdGradUPtr_);
+}
+
+
 // ************************************************************************* //

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

@@ -120,6 +120,8 @@ objectiveMoment::objectiveMoment
     bdSdbMultPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
     bdxdbMultPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
     bdxdbDirectMultPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
+    bdJdnutPtr_.reset(createZeroBoundaryPtr<scalar>(mesh_));
+  //bdJdGradUPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
 }
 
 
@@ -306,6 +308,25 @@ void objectiveMoment::update_dxdbDirectMultiplier()
 }
 
 
+void objectiveMoment::update_boundarydJdnut()
+{
+    const volVectorField& U = vars_.U();
+    volSymmTensorField devGradU(dev(twoSymm(fvc::grad(U))));
+
+    for (const label patchI : momentPatches_)
+    {
+        const fvPatch& patch = mesh_.boundary()[patchI];
+        tmp<vectorField> nf(patch.nf());
+        const vectorField dx(patch.Cf() - rotationCentre_);
+        bdJdnutPtr_()[patchI] =
+           -rhoInf_
+           *(
+                (dx^(devGradU.boundaryField()[patchI] & nf)) & momentDirection_
+            )*invDenom_;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace objectives

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

@@ -118,6 +118,14 @@ public:
         //- Update delta(x)/delta b multiplier coming directly from the
         //- objective
         void update_dxdbDirectMultiplier();
+
+        //- Update dJ/dnut multiplier
+        void update_boundarydJdnut();
+
+        //- Update dJ/dGradU multiplier
+        /* WIP
+        void update_boundarydJdGradU();
+        */
 };
 
 

src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.C

@@ -156,7 +156,6 @@ objective::objective
     bdxdbMultPtr_(nullptr),
     bdxdbDirectMultPtr_(nullptr),
     bEdgeContribution_(nullptr),
-    bdJdStressPtr_(nullptr),
     divDxDbMultPtr_(nullptr),
     gradDxDbMultPtr_(nullptr),
 
@@ -376,10 +375,6 @@ void objective::doNormalization()
         {
             gradDxDbMultPtr_() *= oneOverNorm;
         }
-        if (hasBoundarydJdStress())
-        {
-            bdJdStressPtr_() *= oneOverNorm;
-        }
     }
 }
 
@@ -507,16 +502,6 @@ const vectorField& objective::boundaryEdgeMultiplier
 }
 
 
-const fvPatchTensorField& objective::boundarydJdStress(const label patchI)
-{
-    if (!bdJdStressPtr_)
-    {
-        bdJdStressPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
-    }
-    return bdJdStressPtr_()[patchI];
-}
-
-
 const boundaryVectorField& objective::boundarydJdb()
 {
     if (!bdJdbPtr_)
@@ -580,16 +565,6 @@ const vectorField3& objective::boundaryEdgeMultiplier()
 }
 
 
-const boundaryTensorField& objective::boundarydJdStress()
-{
-    if (!bdJdStressPtr_)
-    {
-        bdJdStressPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
-    }
-    return *bdJdStressPtr_;
-}
-
-
 const volScalarField& objective::divDxDbMultiplier()
 {
     if (!divDxDbMultPtr_)
@@ -666,10 +641,6 @@ void objective::nullify()
                 field = vector::zero;
             }
         }
-        if (hasBoundarydJdStress())
-        {
-            bdJdStressPtr_() == tensor::zero;
-        }
         if (hasDivDxDbMult())
         {
             divDxDbMultPtr_() ==

src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.H

@@ -127,9 +127,6 @@ protected:
         //- NURBS surfaces G1 continuity
         autoPtr<vectorField3> bEdgeContribution_;
 
-        //- For use with discrete-like sensitivities
-        autoPtr<boundaryTensorField> bdJdStressPtr_;
-
         // Contribution to volume-based sensitivities from volume-based
         // objective functions
         //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -298,9 +295,6 @@ public:
             const label edgeI
         );
 
-        //- Objective partial deriv wrt stress tensor
-        const fvPatchTensorField& boundarydJdStress(const label);
-
         //- Contribution to surface sensitivities for all patches
         const boundaryVectorField& boundarydJdb();
 
@@ -319,9 +313,6 @@ public:
         //- Multiplier located at patch boundary edges
         const vectorField3& boundaryEdgeMultiplier();
 
-        //- Objective partial deriv wrt stress tensor
-        const boundaryTensorField& boundarydJdStress();
-
         //- Multiplier of grad( delta(x)/delta b) for volume-based sensitivities
         const volScalarField& divDxDbMultiplier();
 
@@ -363,10 +354,6 @@ public:
         virtual void update_boundaryEdgeContribution()
         {}
 
-        //- Update dJ/dStress field
-        virtual void update_dJdStressMultiplier()
-        {}
-
         //- Update div( dx/db multiplier). Volume-based sensitivity term
         virtual void update_divDxDbMultiplier()
         {}
@@ -415,7 +402,6 @@ public:
         inline bool hasdxdbMult() const;
         inline bool hasdxdbDirectMult() const;
         inline bool hasBoundaryEdgeContribution() const;
-        inline bool hasBoundarydJdStress() const;
         inline bool hasDivDxDbMult() const;
         inline bool hasGradDxDbMult() const;
 

src/optimisation/adjointOptimisation/adjoint/objectives/objective/objectiveI.H

@@ -90,12 +90,6 @@ inline bool Foam::objective::hasGradDxDbMult() const
 }
 
 
-inline bool Foam::objective::hasBoundarydJdStress() const
-{
-    return bool(bdJdStressPtr_);
-}
-
-
 inline bool Foam::objective::hasIntegrationStartTime() const
 {
     return bool(integrationStartTimePtr_);