From 8127ec6a05d994a3da232aab45faa9c246feab27 Mon Sep 17 00:00:00 2001
From: Vaggelis Papoutsis <vaggelisp@gmail.com>
Date: Wed, 23 Nov 2022 15:51:15 +0200
Subject: [PATCH] ENH: added some auxiliary functions in objective

to help allocate pointers related to contributions to the adjoint
turbulence model PDEs, populate them and check the validity of the
cellZones provided for cellZone-based objectives
---
 .../objectiveIncompressible.C                 | 111 ++++++++++++++++++
 .../objectiveIncompressible.H                 |  22 ++++
 2 files changed, 133 insertions(+)

diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C
index 6c71f1e66a..bf59c4d3b4 100644
--- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C
@@ -29,6 +29,7 @@ License
 
 #include "objectiveIncompressible.H"
 #include "incompressiblePrimalSolver.H"
+#include "incompressibleAdjointSolver.H"
 #include "createZeroField.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -611,6 +612,116 @@ void objectiveIncompressible::nullify()
 }
 
 
+void objectiveIncompressible::allocatedJdTurbulence()
+{
+    // Figure out the availability of the adjoint turbulence model variables
+    // through their primal counterparts, since the contructor of the adjoint
+    // solver has not been completed yet
+    const incompressiblePrimalSolver& primSolver =
+        mesh_.lookupObject<incompressiblePrimalSolver>(primalSolverName_);
+    const autoPtr<incompressible::RASModelVariables>& rasVars =
+        primSolver.getIncoVars().RASModelVariables();
+
+    if (rasVars().hasTMVar1())
+    {
+        const dimensionSet primalVarDims = rasVars->TMVar1Inst().dimensions();
+        const dimensionSet sourceDims = dimArea/pow3(dimTime)/primalVarDims;
+        dJdTMvar1Ptr_.reset
+        (
+            createZeroFieldPtr<scalar>
+            (
+                mesh_,
+                "ATMSource1",
+                sourceDims
+            )
+        );
+    }
+    if (rasVars().hasTMVar2())
+    {
+        const dimensionSet primalVarDims = rasVars->TMVar2Inst().dimensions();
+        const dimensionSet sourceDims = dimArea/pow3(dimTime)/primalVarDims;
+        dJdTMvar2Ptr_.reset
+        (
+            createZeroFieldPtr<scalar>
+            (
+                mesh_,
+                "ATMSource2",
+                sourceDims
+            )
+        );
+    }
+}
+
+
+void objectiveIncompressible::checkCellZonesSize
+(
+    const labelList& zoneIDs
+) const
+{
+    label nCells(0);
+    for (const label zI : zoneIDs)
+    {
+        nCells += mesh_.cellZones()[zI].size();
+    }
+    reduce(nCells, sumOp<label>());
+    if (!nCells)
+    {
+        FatalErrorInFunction
+            << "Provided cellZones include no cells"
+            << exit(FatalError);
+    }
+}
+
+
+void objectiveIncompressible::update_dJdTMvar
+(
+    autoPtr<volScalarField>& dJdTMvarPtr,
+    tmp<volScalarField>
+        (incompressibleAdjoint::adjointRASModel::*JacobianFunc)() const,
+    const volScalarField& JacobianMultiplier,
+    const labelList& zones
+)
+{
+    if (dJdTMvarPtr)
+    {
+        // nut Jacobians are currently computed in the adjoint turbulence
+        // models, though they would be better placed within the primal
+        // turbulence model.
+        // Safeguard the computation until the adjoint solver is complete
+        if (mesh_.foundObject<incompressibleAdjointSolver>(adjointSolverName_))
+        {
+            const incompressibleAdjointSolver& adjSolver =
+                mesh_.lookupObject<incompressibleAdjointSolver>
+                    (adjointSolverName_);
+            const autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS =
+                adjSolver.getAdjointVars().adjointTurbulence();
+
+            tmp<volScalarField> tnutJacobian = (adjointRAS->*JacobianFunc)();
+            const volScalarField& nutJacobian = tnutJacobian();
+
+            volScalarField& dJdTMvar = dJdTMvarPtr();
+
+            for (const label zI : zones)
+            {
+                const cellZone& zoneI = mesh_.cellZones()[zI];
+                for (const label cellI : zoneI)
+                {
+                    dJdTMvar[cellI] =
+                        JacobianMultiplier[cellI]*nutJacobian[cellI];
+                }
+            }
+        }
+        else
+        {
+            WarningInFunction
+                << "Skipping the computation of nutJacobian until "
+                << "the adjoint solver is complete"
+                << endl;
+        }
+    }
+}
+
+
 bool objectiveIncompressible::write(const bool valid) const
 {
     return objective::write(valid);
diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H
index b7a01154c1..eb2b89fae9 100644
--- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H
@@ -42,6 +42,7 @@ SourceFiles
 
 #include "objective.H"
 #include "incompressibleVars.H"
+#include "adjointRASModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -252,6 +253,27 @@ public:
         //- Update objective function derivatives
         virtual void nullify();
 
+        // Funtions used in objectives including turbulent quantities
+
+            //- Allocate fields related to the differentiation of turbulence
+            //- models, if necessary
+            void allocatedJdTurbulence();
+
+            //- Check if cellZones provided include at least one cell
+            void checkCellZonesSize(const labelList& zoneIDs) const;
+
+            //- Compute dJdTMVar{1,2}
+            void update_dJdTMvar
+            (
+                autoPtr<volScalarField>& dJdTMvarPtr,
+                tmp<volScalarField>
+                    (incompressibleAdjoint::adjointRASModel::*JacobianFunc)()
+                    const,
+                const volScalarField& JacobianMultiplier,
+                const labelList& zones
+            );
+
+
         //- Update vol and boundary fields and derivatives
         //  Do nothing in the base. The relevant ones should be overwritten
         //  in the child objective functions