ENH: adjustments to the efficiency of the adjoint code

- ATCstandard, ATCUaGradU:
  the ATC is now added as a dimensioned field and not as an fvMatrix
  to UaEqn. This get rid of many unnecessary allocations.

- ATCstandard:
  gradU is cached within the class to avoid its re-computation in
  every adjoint iteration of the steady state solver.

- Inlined a number of functions within the primal and adjoint solvers.
  This probably has a negligible effect since they likely were inlined
  by the compiler either way.

- The momentum diffusivity at the boundary, used by the adjoint boundary
  conditions, was computed for the entire field and, then, only the
  boundary field of each adjoint boundary condition was used. If many
  outlet boundaries exist, the entire nuEff field would be computed as
  many times as the number of boundaries, leading to an unnecessary
  computational overhead.

- Outlet boundary conditions (both pressure and velocity) use the local
  patch gradient to compute their fluxes. This patch gradient requires
  the computation of the adjacent cell gradient, which is done on the
  fly, on a per patch basis. To compute this patch adjacent gradient
  however, the field under the grad sign is interpolated on the entire
  mesh. If many outlets exist, this leads to a huge computational
  overhead. Solved by caching the interpolated field to the database and
  re-using it, in a way similar to the caching of gradient fields (see
  fvc::grad).

WIP: functions returning references to primal and adjoint boundary
fields within boundaryAdjointContributions seem to have a non-negligible
overhead for cases with many patches. No easy work-around here since
these are virtual and cannot be inlined.

WIP: introduced the code structure for caching the contributions to
the adjoint boundary conditions that depend only on the primal fields
and reusing. The process needs to be completed and evaluated, to make
sure that the extra code complexity is justified by gains in
performance.

src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -259,6 +259,12 @@ bool ATCModel::writeData(Ostream&) const
 }
 
 
+void ATCModel::updatePrimalBasedQuantities()
+{
+    // Does nothing in base
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -201,6 +201,9 @@ public:
 
         //- Dummy writeData function required from regIOobject
         virtual bool writeData(Ostream&) const;
+
+        //- Update quantities related with the primal fields
+        virtual void updatePrimalBasedQuantities();
 };
 
 

src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCUaGradU/ATCUaGradU.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -97,7 +97,7 @@ void ATCUaGradU::addATC(fvVectorMatrix& UaEqn)
     smoothATC();
 
     // Actual ATC term
-    UaEqn += fvm::Su(ATC_, Ua);
+    UaEqn += ATC_.internalField();
 }
 
 

src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,20 @@ ATCstandard::ATCstandard
     const dictionary& dict
 )
 :
-    ATCModel(mesh, primalVars, adjointVars, dict)
+    ATCModel(mesh, primalVars, adjointVars, dict),
+    gradU_
+    (
+        IOobject
+        (
+            "gradUATC",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedTensor(dimless/dimTime, Zero)
+    )
 {}
 
 
@@ -65,23 +78,14 @@ ATCstandard::ATCstandard
 
 void ATCstandard::addATC(fvVectorMatrix& UaEqn)
 {
+    addProfiling(ATCstandard, "ATCstandard::addATC");
     const volVectorField& U = primalVars_.U();
     const volVectorField& Ua = adjointVars_.UaInst();
     const surfaceScalarField& phi = primalVars_.phi();
 
-    // Build U to go into the ATC term, based on whether to smooth field or not
-    autoPtr<volVectorField> UForATC(nullptr);
-    if (reconstructGradients_)
-    {
-        UForATC.reset(new volVectorField(fvc::reconstruct(phi)));
-    }
-    else
-    {
-        UForATC.reset(new volVectorField(U));
-    }
 
     // Main ATC term
-    ATC_ = (fvc::grad(UForATC(), "gradUATC") & Ua);
+    ATC_ = gradU_ & Ua;
 
     if (extraConvection_ > 0)
     {
@@ -97,55 +101,62 @@ void ATCstandard::addATC(fvVectorMatrix& UaEqn)
     smoothATC();
 
     // actual ATC term
-    UaEqn += fvm::Su(ATC_, Ua);
+    UaEqn += ATC_.internalField();
 }
 
 
 tmp<volTensorField> ATCstandard::getFISensitivityTerm() const
 {
-    tmp<volTensorField> tvolSDTerm
-    (
-        new volTensorField
-        (
-            IOobject
-            (
-                "ATCFISensitivityTerm" + type(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedTensor(sqr(dimLength)/pow(dimTime, 3), Zero)
-        )
-    );
-
-    volTensorField& volSDTerm = tvolSDTerm.ref();
-
     const volVectorField& U = primalVars_.U();
     const volVectorField& Ua = adjointVars_.Ua();
 
-    volTensorField gradU(fvc::grad(U));
+    // We only need to modify the boundaryField of gradU locally.
+    // If grad(U) is cached then
+    // a. The .ref() call fails since the tmp is initialised from a
+    //    const ref
+    // b. we would be changing grad(U) for all other places in the code
+    //    that need it
+    // So, always allocate new memory and avoid registering the new field
+    tmp<volTensorField> tgradU(volTensorField::New("gradULocal", fvc::grad(U)));
+    volTensorField::Boundary& gradUbf = tgradU.ref().boundaryFieldRef();
 
-    // Explicitly correct the boundary gradient to get rid of the
+    // Explicitly correct the boundary gradient to get rid of
-    // tangential component
+    // 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();
+            gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad();
-            gradU.boundaryFieldRef()[patchI] =
-                nf*U.boundaryField()[patchI].snGrad();
         }
     }
 
     const volScalarField& mask = getLimiter();
 
-    volSDTerm = -(gradU & Ua)*U*mask;
+    return
+        tmp<volTensorField>::New
+        (
+            "ATCFISensitivityTerm" + type(),
+          - (tgradU & Ua)*U*mask
+        );
+}
 
-    return tvolSDTerm;
+
+void ATCstandard::updatePrimalBasedQuantities()
+{
+    const volVectorField& U = primalVars_.U();
+    const surfaceScalarField& phi = primalVars_.phi();
+    // Build U to go into the ATC term, based on whether to smooth field or not
+    autoPtr<volVectorField> UForATC(nullptr);
+    if (reconstructGradients_)
+    {
+        gradU_ = fvc::grad(fvc::reconstruct(phi), "gradUATC");
+    }
+    else
+    {
+        gradU_ = fvc::grad(U, "gradUATC");
+    }
 }
 
 

src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class ATCstandard Declaration
+                         Class ATCstandard Declaration
 \*---------------------------------------------------------------------------*/
 
 class ATCstandard
@@ -58,6 +58,12 @@ class ATCstandard
 {
 private:
 
+    // Private data
+
+        //- The gradU used in the computation of the standard ATC
+        //  Cached to avoid costly recomputation in each adjoint iteration
+        volTensorField gradU_;
+
     // Private Member Functions
 
         //- No copy construct
@@ -96,6 +102,9 @@ public:
 
         //- Get the FI sensitivity derivatives term coming from the ATC
         virtual tmp<volTensorField> getFISensitivityTerm() const;
+
+        //- Update quantities related with the primal fields
+        virtual void updatePrimalBasedQuantities();
 };
 
 

src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -73,19 +73,48 @@ adjointBoundaryCondition<Type>::computePatchGrad(word name)
     word schemeData(is);
     */
 
-    tmp<surfaceInterpolationScheme<Type2>> tinterpScheme
+    // If there are many patches calling this function, the computation of
-    (
+    // the surface field might be a significant computational
-        surfaceInterpolationScheme<Type2>::New
+    // burden. Cache the interpolated field and fetch from the registry when
-        (
+    // possible.
-            mesh,
+    const word surfFieldName("interpolated" + name + "ForBoundaryGrad");
-            mesh.interpolationScheme("interpolate(" + name + ")")
+    typedef GeometricField<Type2, fvsPatchField, surfaceMesh> surfFieldType;
-        )
+    surfFieldType* surfFieldPtr =
-    );
+        mesh.objectRegistry::template getObjectPtr<surfFieldType>(surfFieldName);
 
-    GeometricField<Type2, fvsPatchField, surfaceMesh> surfField
+    if (!surfFieldPtr)
-    (
+    {
-        tinterpScheme().interpolate(field)
+        solution::cachePrintMessage("Calculating and caching", name, field);
-    );
+
+        surfFieldPtr =
+            surfaceInterpolationScheme<Type2>::New
+            (
+                mesh, mesh.interpolationScheme("interpolate(" + name + ")")
+            )().interpolate(field).ptr();
+        surfFieldPtr->rename(surfFieldName);
+        regIOobject::store(surfFieldPtr);
+    }
+    else
+    {
+        if (surfFieldPtr->upToDate(field))
+        {
+            solution::cachePrintMessage("Reusing", name, field);
+        }
+        else
+        {
+            solution::cachePrintMessage("Updating", name, field);
+            delete surfFieldPtr;
+
+            surfFieldPtr =
+                surfaceInterpolationScheme<Type2>::New
+                (
+                    mesh, mesh.interpolationScheme("interpolate(" + name + ")")
+                )().interpolate(field).ptr();
+            surfFieldPtr->rename(surfFieldName);
+            regIOobject::store(surfFieldPtr);
+        }
+    }
+    surfFieldType& surfField = *surfFieldPtr;
 
     // Auxiliary fields
     const surfaceVectorField& Sf = mesh.Sf();
@@ -271,6 +300,13 @@ const ATCModel& adjointBoundaryCondition<Type>::getATC() const
 }
 
 
+template<class Type>
+void adjointBoundaryCondition<Type>::updatePrimalBasedQuantities()
+{
+    // Does nothing in base
+}
+
+
 template<class Type>
 tmp
 <

src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -152,6 +152,10 @@ public:
             <
                 Field<typename Foam::outerProduct<Foam::vector, Type>::type>
             > dxdbMult() const;
+
+        //- Update the primal based quantities related to the adjoint boundary
+        //- conditions
+        virtual void updatePrimalBasedQuantities();
 };
 
 

src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -97,6 +97,11 @@ void Foam::adjointOutletVelocityFluxFvPatchVectorField::manipulateMatrix
     fvMatrix<vector>& matrix
 )
 {
+    addProfiling
+    (
+        adjointOutletVelocityFluxFvPatchVectorField,
+        "adjointOutletVelocityFluxFvPatchVectorField::manipulateMatrix"
+    );
     vectorField& source = matrix.source();
     const vectorField& Sf = patch().Sf();
     const labelList& faceCells = patch().faceCells();

src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -63,10 +63,6 @@ class adjointOutletVelocityFluxFvPatchVectorField
     public fixedValueFvPatchVectorField,
     public adjointVectorBoundaryCondition
 {
-    // Private Member Functions
-
-        tmp<tensorField> computeLocalGrad();
-
 
 public:
 

src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointWallVelocity/adjointWallVelocityFvPatchVectorField.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,11 @@ void Foam::adjointWallVelocityFvPatchVectorField::manipulateMatrix
     fvMatrix<vector>& matrix
 )
 {
+    addProfiling
+    (
+        adjointWallVelocityFvPatchVectorField,
+        "adjointWallVelocityFvPatchVectorField::manipulateMatrix"
+    );
     // Grab ref to the diagonal matrix
     vectorField& source = matrix.source();
 

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -89,11 +89,10 @@ boundaryAdjointContributionIncompressible
 tmp<vectorField> boundaryAdjointContributionIncompressible::velocitySource()
 {
     // Objective function contribution
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
     tmp<vectorField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdv
         );
     vectorField& source = tsource.ref();
@@ -110,11 +109,10 @@ tmp<vectorField> boundaryAdjointContributionIncompressible::velocitySource()
 tmp<scalarField> boundaryAdjointContributionIncompressible::pressureSource()
 {
     // Objective function contribution
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
     tmp<scalarField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdvn
         );
 
@@ -126,10 +124,7 @@ tmp<scalarField> boundaryAdjointContributionIncompressible::pressureSource()
     const vectorField& adjointTurbulenceContr =
         adjointRAS().adjointMomentumBCSource()[patch_.index()];
 
-    tmp<vectorField> tnf = patch_.nf();
+    source += adjointTurbulenceContr & patch_.nf();
-    const vectorField& nf = tnf();
-
-    source += adjointTurbulenceContr & nf;
 
     return (tsource);
 }
@@ -139,11 +134,10 @@ tmp<vectorField>
 boundaryAdjointContributionIncompressible::tangentVelocitySource()
 {
     // Objective function contribution
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
     tmp<vectorField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdvt
         );
 
@@ -167,59 +161,51 @@ boundaryAdjointContributionIncompressible::tangentVelocitySource()
 tmp<vectorField>
 boundaryAdjointContributionIncompressible::normalVelocitySource()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
+    return
-    tmp<vectorField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdp
         );
 
-    return (tsource);
 }
 
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::energySource()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
+    return
-    tmp<scalarField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdT
         );
 
-    return (tsource);
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::adjointTMVariable1Source()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
+    return
-    tmp<scalarField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdTMvar1
         );
 
-    return (tsource);
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::adjointTMVariable2Source()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
+    return
-    tmp<scalarField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdTMvar2
         );
 
-    return (tsource);
 }
 
 
@@ -249,15 +235,8 @@ boundaryAdjointContributionIncompressible::dJdGradU()
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::momentumDiffusion()
 {
-    tmp<scalarField> tnuEff(new scalarField(patch_.size(), Zero));
-    scalarField& nuEff = tnuEff.ref();
 
-    const autoPtr<incompressibleAdjoint::adjointRASModel>&
+    return adjointVars().adjointTurbulence()().nuEff(patch_.index());
-        adjointTurbulenceModel = adjointVars().adjointTurbulence();
-
-    nuEff = adjointTurbulenceModel().nuEff()().boundaryField()[patch_.index()];
-
-    return tnuEff;
 }
 
 
@@ -296,64 +275,40 @@ tmp<scalarField> boundaryAdjointContributionIncompressible::thermalDiffusion()
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::wallDistance()
 {
-    tmp<scalarField> twallDist(new scalarField(patch_.size(), Zero));
+    return primalVars_.turbulence()->y()[patch_.index()];
-    scalarField& wallDist = twallDist.ref();
-
-    wallDist = primalVars_.turbulence()->y()[patch_.index()];
-
-    return twallDist;
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::TMVariable1Diffusion()
 {
-    const autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS =
+    return
-        adjointVars().adjointTurbulence();
+        adjointVars().adjointTurbulence()->diffusionCoeffVar1(patch_.index());
 
-    tmp<scalarField> tdiffCoeff =
-        adjointRAS().diffusionCoeffVar1(patch_.index());
-
-    return tdiffCoeff;
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::TMVariable2Diffusion()
 {
-    const autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS =
+    return
-        adjointVars().adjointTurbulence();
+        adjointVars().adjointTurbulence()->diffusionCoeffVar2(patch_.index());
-
-    tmp<scalarField> tdiffCoeff =
-        adjointRAS().diffusionCoeffVar2(patch_.index());
-
-    return tdiffCoeff;
 }
 
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::TMVariable1()
 {
-    const autoPtr<incompressible::RASModelVariables>& RASVariables =
+    return
-        primalVars_.RASModelVariables();
+        primalVars_.RASModelVariables()->TMVar1().
-    tmp<scalarField> tboundField(new scalarField(patch_.size(), Zero));
+            boundaryField()[patch_.index()];
-    scalarField& boundField = tboundField.ref();
-
-    boundField = RASVariables().TMVar1().boundaryField()[patch_.index()];
-
-    return tboundField;
 }
 
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::TMVariable2()
 {
-    const autoPtr<incompressible::RASModelVariables>& RASVariables =
+    return
-        primalVars_.RASModelVariables();
+        primalVars_.RASModelVariables()->TMVar2().
-    tmp<scalarField> tboundField(new scalarField(patch_.size(), Zero));
+            boundaryField()[patch_.index()];
-    scalarField& boundField = tboundField.ref();
-
-    boundField = RASVariables().TMVar2().boundaryField()[patch_.index()];
-
-    return tboundField;
 }
 
 

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2019, 2022 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2019, 2022 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -70,28 +70,7 @@ objectiveForce::objectiveForce
     forceDirection_(dict.get<vector>("direction")),
     Aref_(dict.get<scalar>("Aref")),
     rhoInf_(dict.get<scalar>("rhoInf")),
-    UInf_(dict.get<scalar>("UInf")),
+    UInf_(dict.get<scalar>("UInf"))
-    stressXPtr_
-    (
-        Foam::createZeroFieldPtr<vector>
-        (
-            mesh_, "stressX", dimLength/sqr(dimTime)
-        )
-    ),
-    stressYPtr_
-    (
-        Foam::createZeroFieldPtr<vector>
-        (
-            mesh_, "stressY", dimLength/sqr(dimTime)
-        )
-    ),
-    stressZPtr_
-    (
-        Foam::createZeroFieldPtr<vector>
-        (
-            mesh_, "stressZ", dimLength/sqr(dimTime)
-        )
-    )
 {
     // Sanity check and print info
     if (forcePatches_.empty())
@@ -135,16 +114,9 @@ scalar objectiveForce::J()
 
     for (const label patchI : forcePatches_)
     {
-        pressureForce += gSum
+        const vectorField& Sf = mesh_.Sf().boundaryField()[patchI];
-        (
+        pressureForce += gSum(Sf*p.boundaryField()[patchI]);
-            mesh_.Sf().boundaryField()[patchI] * p.boundaryField()[patchI]
+        viscousForce += gSum(devReff.boundaryField()[patchI] & Sf);
-        );
-        // Viscous term calculated using the full tensor derivative
-        viscousForce += gSum
-        (
-            devReff.boundaryField()[patchI]
-          & mesh_.Sf().boundaryField()[patchI]
-        );
     }
 
     cumulativeForce = pressureForce + viscousForce;
@@ -178,8 +150,8 @@ void objectiveForce::update_dSdbMultiplier()
     // Compute contributions with mean fields, if present
     const volScalarField& p = vars_.p();
     const volVectorField& U = vars_.U();
-    const autoPtr<incompressible::RASModelVariables>&
+    const autoPtr<incompressible::RASModelVariables>& turbVars =
-        turbVars = vars_.RASModelVariables();
+        vars_.RASModelVariables();
     const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
 
     tmp<volSymmTensorField> tdevReff = turbVars->devReff(lamTransp, U);
@@ -208,11 +180,16 @@ void objectiveForce::update_dxdbMultiplier()
         turbVars = vars_.RASModelVariables();
     const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
 
-    //tmp<volSymmTensorField> tdevReff = turbVars->devReff(lamTransp, U);
+    // We only need to modify the boundaryField of gradU locally.
-    //const volSymmTensorField& devReff = tdevReff();
+    // If grad(U) is cached then
-
+    // a. The .ref() call fails since the tmp is initialised from a
-    volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
+    //    const ref
-    volTensorField gradU(fvc::grad(U));
+    // b. we would be changing grad(U) for all other places in the code
+    //    that need it
+    // So, always allocate new memory and avoid registering the new field
+    tmp<volTensorField> tgradU =
+        volTensorField::New("gradULocal", fvc::grad(U));
+    volTensorField& gradU = tgradU.ref();
     volTensorField::Boundary& gradUbf = gradU.boundaryFieldRef();
 
     // Explicitly correct the boundary gradient to get rid of
@@ -222,49 +199,42 @@ void objectiveForce::update_dxdbMultiplier()
         const fvPatch& patch = mesh_.boundary()[patchI];
         if (isA<wallFvPatch>(patch))
         {
-            tmp<vectorField> nf = patch.nf();
+            tmp<vectorField> tnf = patch.nf();
-            gradUbf[patchI] = nf*U.boundaryField()[patchI].snGrad();
+            gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad();
         }
     }
 
-    volTensorField stress(nuEff*(gradU + T(gradU)));
+    // Term coming from gradp
-
+    tmp<volVectorField> tgradp(fvc::grad(p));
-    stressXPtr_().replace(0, stress.component(0));
+    const volVectorField& gradp = tgradp.cref();
-    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_()));
-
-    // the notorious second-order derivative at the wall. Use with caution!
-    volVectorField gradp(fvc::grad(p));
-
     for (const label patchI : forcePatches_)
     {
-        const fvPatch& patch = mesh_.boundary()[patchI];
-        tmp<vectorField> tnf = patch.nf();
-        const vectorField& nf = tnf();
         bdxdbMultPtr_()[patchI] =
-        (
+            (forceDirection_ & mesh_.boundary()[patchI].nf())
-            (
+           *gradp.boundaryField()[patchI]/denom();
-                (
+    }
-                   -(forceDirection_.x() * gradStressX.boundaryField()[patchI])
+    tgradp.clear();
-                   -(forceDirection_.y() * gradStressY.boundaryField()[patchI])
+
-                   -(forceDirection_.z() * gradStressZ.boundaryField()[patchI])
+    // Term coming from stresses
-                ) & nf
+    tmp<volScalarField> tnuEff = lamTransp.nu() + turbVars->nutRef();
-            )
+    tmp<volSymmTensorField> tstress = tnuEff*twoSymm(tgradU);
-            + (forceDirection_ & nf)*gradp.boundaryField()[patchI]
+    const volSymmTensorField& stress = tstress.cref();
-        )
+    autoPtr<volVectorField> ptemp
-        /denom();
+        (Foam::createZeroFieldPtr<vector>( mesh_, "temp", sqr(dimVelocity)));
+    volVectorField& temp = ptemp.ref();
+
+    for (label idir = 0; idir < pTraits<vector>::nComponents; ++idir)
+    {
+        unzipRow(stress, idir, temp);
+        volTensorField gradStressDir(fvc::grad(temp));
+        for (const label patchI : forcePatches_)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            tmp<vectorField> tnf = patch.nf();
+            bdxdbMultPtr_()[patchI] -=
+                forceDirection_.component(idir)
+               *(gradStressDir.boundaryField()[patchI] & tnf)/denom();
+        }
     }
 }
 
@@ -278,17 +248,17 @@ void objectiveForce::update_boundarydJdnut()
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
         tmp<vectorField> tnf = patch.nf();
-        const vectorField& nf = tnf();
         bdJdnutPtr_()[patchI] =
-           -((devGradU.boundaryField()[patchI] & forceDirection_) & nf)/denom();
+          - ((devGradU.boundaryField()[patchI] & forceDirection_) & tnf)
+           /denom();
     }
 }
 
 
 void objectiveForce::update_boundarydJdGradU()
 {
-    const autoPtr<incompressible::RASModelVariables>&
+    const autoPtr<incompressible::RASModelVariables>& turbVars =
-        turbVars = vars_.RASModelVariables();
+        vars_.RASModelVariables();
     const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
     volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
     for (const label patchI : forcePatches_)

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2020, 2022 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2020, 2022 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -70,10 +70,6 @@ protected:
         scalar rhoInf_;
         scalar UInf_;
 
-        autoPtr<volVectorField> stressXPtr_;
-        autoPtr<volVectorField> stressYPtr_;
-        autoPtr<volVectorField> stressZPtr_;
-
 
 public:
 

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2019, 2022 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2019, 2022 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -76,27 +76,6 @@ objectiveMoment::objectiveMoment
     rhoInf_(dict.get<scalar>("rhoInf")),
     UInf_(dict.get<scalar>("UInf")),
     invDenom_(2./(rhoInf_*UInf_*UInf_*Aref_*lRef_)),
-    stressXPtr_
-    (
-        Foam::createZeroFieldPtr<vector>
-        (
-            mesh_, "stressX", dimLength/sqr(dimTime)
-        )
-    ),
-    stressYPtr_
-    (
-        Foam::createZeroFieldPtr<vector>
-        (
-            mesh_, "stressY", dimLength/sqr(dimTime)
-        )
-    ),
-    stressZPtr_
-    (
-        Foam::createZeroFieldPtr<vector>
-        (
-            mesh_, "stressZ", dimLength/sqr(dimTime)
-        )
-    ),
     devReff_(vars_.turbulence()->devReff()())
 {
     // Sanity check and print info
@@ -170,8 +149,8 @@ void objectiveMoment::update_meanValues()
     if (computeMeanFields_)
     {
         const volVectorField& U = vars_.U();
-        const autoPtr<incompressible::RASModelVariables>&
+        const autoPtr<incompressible::RASModelVariables>& turbVars =
-           turbVars = vars_.RASModelVariables();
+            vars_.RASModelVariables();
         const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
 
         devReff_ = turbVars->devReff(lamTransp, U)();
@@ -184,8 +163,8 @@ void objectiveMoment::update_boundarydJdp()
     for (const label patchI : momentPatches_)
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
-        vectorField dx(patch.Cf() - rotationCentre_);
+        tmp<vectorField> tdx = patch.Cf() - rotationCentre_;
-        bdJdpPtr_()[patchI] = (momentDirection_ ^ dx)*invDenom_*rhoInf_;
+        bdJdpPtr_()[patchI] = (momentDirection_ ^ tdx)*invDenom_*rhoInf_;
     }
 }
 
@@ -197,19 +176,19 @@ void objectiveMoment::update_dSdbMultiplier()
     for (const label patchI : momentPatches_)
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
-        const vectorField dx(patch.Cf() - rotationCentre_);
+        tmp<vectorField> tdx = patch.Cf() - rotationCentre_;
         bdSdbMultPtr_()[patchI] =
             (
                 (
                     rhoInf_*
                     (
-                        (momentDirection_^dx) &
+                        (momentDirection_ ^ tdx()) &
                         (
                             devReff_.boundaryField()[patchI]
                         )
                     )
                 )
-              + rhoInf_ * (momentDirection_^dx) * p.boundaryField()[patchI]
+              + rhoInf_*(momentDirection_ ^ tdx())*p.boundaryField()[patchI]
             )
            *invDenom_;
     }
@@ -221,11 +200,20 @@ void objectiveMoment::update_dxdbMultiplier()
     const volScalarField& p = vars_.p();
     const volVectorField& U = vars_.U();
 
-    const autoPtr<incompressible::RASModelVariables>&
+    const autoPtr<incompressible::RASModelVariables>& turbVars =
-       turbVars = vars_.RASModelVariables();
+        vars_.RASModelVariables();
     const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
-    volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
+
-    volTensorField gradU(fvc::grad(U));
+    // We only need to modify the boundaryField of gradU locally.
+    // If grad(U) is cached then
+    // a. The .ref() call fails since the tmp is initialised from a
+    //    const ref
+    // b. we would be changing grad(U) for all other places in the code
+    //    that need it
+    // So, always allocate new memory and avoid registering the new field
+    tmp<volTensorField> tgradU =
+        volTensorField::New("gradULocal", fvc::grad(U));
+    volTensorField::Boundary& gradUbf = tgradU.ref().boundaryFieldRef();
 
     // Explicitly correct the boundary gradient to get rid of the
     // tangential component
@@ -235,51 +223,46 @@ void objectiveMoment::update_dxdbMultiplier()
         if (isA<wallFvPatch>(patch))
         {
             tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
-            const vectorField& nf = tnf();
+            gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad();
-            gradU.boundaryFieldRef()[patchI] =
-                nf * U.boundaryField()[patchI].snGrad();
         }
     }
 
-    volTensorField stress(nuEff*(gradU + T(gradU)));
+    // Term coming from gradp
-
+    tmp<volVectorField> tgradp = fvc::grad(p);
-    stressXPtr_().replace(0, stress.component(0));
+    const volVectorField& gradp = tgradp.cref();
-    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_()));
-
-    volVectorField gradp(fvc::grad(p));
-
     for (const label patchI : momentPatches_)
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
         tmp<vectorField> tnf = patch.nf();
-        const vectorField& nf = tnf();
+        tmp<vectorField> tdx = patch.Cf() - rotationCentre_;
-        vectorField dx(patch.Cf() - rotationCentre_);
-        vectorField aux(momentDirection_^dx);
         bdxdbMultPtr_()[patchI] =
-        (
+            (momentDirection_ & (tdx ^ tnf))*gradp.boundaryField()[patchI]
-            (
+           *invDenom_*rhoInf_;
-                (
+    }
-                   -(aux.component(0) * gradStressX.boundaryField()[patchI])
+    tgradp.clear();
-                   -(aux.component(1) * gradStressY.boundaryField()[patchI])
+
-                   -(aux.component(2) * gradStressZ.boundaryField()[patchI])
+    // Term coming from stresses
-                ) & nf
+    tmp<volScalarField> tnuEff = lamTransp.nu() + turbVars->nutRef();
-            )
+    tmp<volSymmTensorField> tstress = tnuEff*twoSymm(tgradU);
-            + (momentDirection_ & (dx^nf))*gradp.boundaryField()[patchI]
+    const volSymmTensorField& stress = tstress.cref();
-        )
+    autoPtr<volVectorField> ptemp
-        *invDenom_*rhoInf_;
+        (Foam::createZeroFieldPtr<vector>( mesh_, "temp", sqr(dimVelocity)));
+    volVectorField& temp = ptemp.ref();
+
+    for (label idir = 0; idir < pTraits<vector>::nComponents; ++idir)
+    {
+        unzipRow(stress, idir, temp);
+        volTensorField gradStressDir(fvc::grad(temp));
+        for (const label patchI : momentPatches_)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            tmp<vectorField> tnf = patch.nf();
+            tmp<vectorField> tdx = patch.Cf() - rotationCentre_;
+            tmp<scalarField> taux = (momentDirection_ ^ tdx)().component(idir);
+            bdxdbMultPtr_()[patchI] -=
+                taux*(gradStressDir.boundaryField()[patchI] & tnf)
+               *invDenom_*rhoInf_;
+        }
     }
 }
 
@@ -316,13 +299,12 @@ void objectiveMoment::update_boundarydJdnut()
     for (const label patchI : momentPatches_)
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
-        tmp<vectorField> nf(patch.nf());
+        tmp<vectorField> tnf = patch.nf();
-        const vectorField dx(patch.Cf() - rotationCentre_);
+        tmp<vectorField> tdx = patch.Cf() - rotationCentre_;
+        const fvPatchSymmTensorField& bdevGradU =
+            devGradU.boundaryField()[patchI];
         bdJdnutPtr_()[patchI] =
-           -rhoInf_
+          - rhoInf_*((tdx ^ (bdevGradU & tnf)) & momentDirection_)*invDenom_;
-           *(
-                (dx^(devGradU.boundaryField()[patchI] & nf)) & momentDirection_
-            )*invDenom_;
     }
 }
 

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2019, 2022 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2019, 2022 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -68,10 +68,6 @@ class objectiveMoment
         scalar UInf_;
         scalar invDenom_;
 
-        autoPtr<volVectorField> stressXPtr_;
-        autoPtr<volVectorField> stressYPtr_;
-        autoPtr<volVectorField> stressZPtr_;
-
         // Store this in order to computed only once per objective call
         volSymmTensorField devReff_;
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/adjointMeshMovementSolver/adjointMeshMovementSolverIncompressible.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2007-2022 PCOpt/NTUA
-    Copyright (C) 2013-2021 FOSS GP
+    Copyright (C) 2013-2022 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -76,9 +76,9 @@ adjointMeshMovementSolver::adjointMeshMovementSolver
         (
             word
             (
-                adjointSensitivity.adjointVars().useSolverNameForFields() ?
+                adjointSensitivity.adjointVars().useSolverNameForFields()
-                "ma" + adjointSensitivity.adjointSolver().solverName() :
+              ? "ma" + adjointSensitivity.adjointSolver().solverName()
-                "ma"
+              : "ma"
             ),
             mesh.time().timeName(),
             mesh,

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2007-2022 PCOpt/NTUA
-    Copyright (C) 2013-2021 FOSS GP
+    Copyright (C) 2013-2022 FOSS GP
     Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -100,24 +100,6 @@ autoPtr<adjointSensitivity> adjointSensitivity::New
 
 // * * * * * * * * * * * * * * *  Member Functions   * * * * * * * * * * * * //
 
-const incompressibleVars& adjointSensitivity::primalVars() const
-{
-    return primalVars_;
-}
-
-
-const incompressibleAdjointVars& adjointSensitivity::adjointVars() const
-{
-    return adjointVars_;
-}
-
-
-const incompressibleAdjointSolver& adjointSensitivity::adjointSolver() const
-{
-    return adjointSolver_;
-}
-
-
 const scalarField& adjointSensitivity::calculateSensitivities()
 {
     assembleSensitivities();

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2007-2022 PCOpt/NTUA
-    Copyright (C) 2013-2021 FOSS GP
+    Copyright (C) 2013-2022 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -156,13 +156,22 @@ public:
     // Member Functions
 
         //- Get primal variables
-        const incompressibleVars& primalVars() const;
+        inline const incompressibleVars& primalVars() const
+        {
+            return primalVars_;
+        }
 
         //- Get adjoint variables
-        const incompressibleAdjointVars& adjointVars() const;
+        inline const incompressibleAdjointVars& adjointVars() const
+        {
+            return adjointVars_;
+        }
 
         //- Get adjoint solver
-        const incompressibleAdjointSolver& adjointSolver() const;
+        inline const incompressibleAdjointSolver& adjointSolver() const
+        {
+            return adjointSolver_;
+        }
 
         //- Accumulate sensitivity integrands
         //  Corresponds to the flow and adjoint part of the sensitivities

src/optimisation/adjointOptimisation/adjoint/parameterization/Bezier/Bezier.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2019, 2022 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2019, 2022 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -236,9 +236,7 @@ tmp<vectorField> Bezier::dndbBasedSensitivities
     const label patchStart = ppatch.start();
     const tensorField& dxdbInt = dxidXj_[cpI].primitiveField();
     vectorField dxdbDir(dxdbInt.size(), Zero);
-    dxdbDir.replace(0, dxdbInt.component(3*idir));
+    unzipRow(dxdbInt, vector::components(idir), dxdbDir);
-    dxdbDir.replace(1, dxdbInt.component(3*idir + 1));
-    dxdbDir.replace(2, dxdbInt.component(3*idir + 2));
 
     // Loop over patch faces
     forAll(patch, fI)

src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2007-2022 PCOpt/NTUA
-    Copyright (C) 2013-2021 FOSS GP
+    Copyright (C) 2013-2022 FOSS GP
     Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -653,17 +653,7 @@ Foam::NURBS3DVolume::NURBS3DVolume
     maxIter_(dict.getOrDefault<label>("maxIterations", 10)),
     tolerance_(dict.getOrDefault<scalar>("tolerance", 1.e-10)),
     nMaxBound_(dict.getOrDefault<scalar>("nMaxBoundIterations", 4)),
-    cps_
+    cps_(),
-    (
-        found("controlPoints") ?
-        vectorField
-        (
-            "controlPoints",
-            *this,
-            basisU_.nCPs()*basisV_.nCPs()*basisW_.nCPs()
-        ) :
-        vectorField(0)
-    ),
     mapPtr_(nullptr),
     reverseMapPtr_(nullptr),
     parametricCoordinatesPtr_(nullptr),
@@ -761,7 +751,17 @@ Foam::NURBS3DVolume::NURBS3DVolume
     }
 
     // Construct control points, if not already read from file
-    if (cps_.empty())
+    if (found("controlPoints"))
+    {
+        cps_ =
+            vectorField
+            (
+                "controlPoints",
+                *this,
+                basisU_.nCPs()*basisV_.nCPs()*basisW_.nCPs()
+            );
+    }
+    else
     {
         controlPointsDefinition::New(*this);
     }
@@ -927,17 +927,7 @@ Foam::tensor Foam::NURBS3DVolume::JacobianUVW
     vector vDeriv = volumeDerivativeV(u, v, w);
     vector wDeriv = volumeDerivativeW(u, v, w);
 
-    tensor Jacobian(Zero);
+    tensor Jacobian(uDeriv, vDeriv, wDeriv, true);
-
-    Jacobian[0] = uDeriv.component(0);
-    Jacobian[1] = vDeriv.component(0);
-    Jacobian[2] = wDeriv.component(0);
-    Jacobian[3] = uDeriv.component(1);
-    Jacobian[4] = vDeriv.component(1);
-    Jacobian[5] = wDeriv.component(1);
-    Jacobian[6] = uDeriv.component(2);
-    Jacobian[7] = vDeriv.component(2);
-    Jacobian[8] = wDeriv.component(2);
 
     return Jacobian;
 }
@@ -1768,7 +1758,7 @@ Foam::tmp<Foam::volTensorField> Foam::NURBS3DVolume::getDxCellsDb
 Foam::label Foam::NURBS3DVolume::nUSymmetry() const
 {
     label nU(basisU_.nCPs());
-    return label(nU % 2 == 0 ? 0.5*nU : (nU - 1)/2 + 1);
+    return label(nU % 2 == 0 ? 0.5*nU : 0.5*(nU - 1) + 1);
 }
 
 

src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/controlPointsDefinition/controlPointsDefinition/controlPointsDefinition.C

@@ -62,9 +62,7 @@ void Foam::controlPointsDefinition::transformControlPoints
     vector scale(dict.get<vector>("scale"));
 
     // Scale box
-    cps_.replace(0, cps_.component(0)*scale.x());
+    cps_ = cmptMultiply(cps_, scale);
-    cps_.replace(1, cps_.component(1)*scale.y());
-    cps_.replace(2, cps_.component(2)*scale.z());
 
     // Rotation matrices
     tensor Rx

src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/adjointSolver/adjointSolver.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -153,6 +153,7 @@ Foam::objectiveManager& Foam::adjointSolver::getObjectiveManager()
 
 void Foam::adjointSolver::postLoop()
 {
+    addProfiling(adjointSolver, "adjointSolver::postLoop");
     computeObjectiveSensitivities();
     // The solver dictionary has been already written after the termination
     // of the adjoint loop. Force re-writing it to include the sensitivities

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -206,6 +206,7 @@ void Foam::adjointSimple::preIter()
 
 void Foam::adjointSimple::mainIter()
 {
+    addProfiling(adjointSimple, "adjointSimple::mainIter");
     // Grab primal references
     const surfaceScalarField& phi = primalVars_.phi();
     // Grab adjoint references
@@ -340,6 +341,7 @@ void Foam::adjointSimple::postIter()
 
 void Foam::adjointSimple::solve()
 {
+    addProfiling(adjointSimple, "adjointSimple::solve");
     if (active_)
     {
         preLoop();

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -163,6 +163,8 @@ void Foam::incompressibleAdjointSolver::updatePrimalBasedQuantities()
     if (vars_)
     {
         getAdjointVars().adjointTurbulence()->setChangedPrimalSolution();
+        ATCModel_().updatePrimalBasedQuantities();
+        getAdjointVars().updatePrimalBasedQuantities();
     }
 }
 

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -379,81 +379,6 @@ surfaceScalarField& incompressibleVars::phi()
 }
 
 
-const volScalarField& incompressibleVars::pInst() const
-{
-    return pPtr_();
-}
-
-
-volScalarField& incompressibleVars::pInst()
-{
-    return pPtr_();
-}
-
-
-const volVectorField& incompressibleVars::UInst() const
-{
-    return UPtr_();
-}
-
-
-volVectorField& incompressibleVars::UInst()
-{
-    return UPtr_();
-}
-
-
-const surfaceScalarField& incompressibleVars::phiInst() const
-{
-    return phiPtr_();
-}
-
-
-surfaceScalarField& incompressibleVars::phiInst()
-{
-    return phiPtr_();
-}
-
-
-const singlePhaseTransportModel& incompressibleVars::laminarTransport() const
-{
-    return laminarTransportPtr_();
-}
-
-
-singlePhaseTransportModel& incompressibleVars::laminarTransport()
-{
-    return laminarTransportPtr_();
-}
-
-
-const autoPtr<incompressible::turbulenceModel>&
-incompressibleVars::turbulence() const
-{
-    return turbulence_;
-}
-
-
-autoPtr<incompressible::turbulenceModel>& incompressibleVars::turbulence()
-{
-    return turbulence_;
-}
-
-
-const autoPtr<incompressible::RASModelVariables>&
-incompressibleVars::RASModelVariables() const
-{
-    return RASModelVariables_;
-}
-
-
-autoPtr<incompressible::RASModelVariables>&
-incompressibleVars::RASModelVariables()
-{
-    return RASModelVariables_;
-}
-
-
 void incompressibleVars::restoreInitValues()
 {
     if (solverControl_.storeInitValues())

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -182,42 +182,44 @@ public:
             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
                 //- Return const reference to pressure
-                const volScalarField& pInst() const;
+                inline const volScalarField& pInst() const;
 
                 //- Return reference to pressure
-                volScalarField& pInst();
+                inline volScalarField& pInst();
 
                 //- Return const reference to velocity
-                const volVectorField& UInst() const;
+                inline const volVectorField& UInst() const;
 
                 //- Return reference to velocity
-                volVectorField& UInst();
+                inline volVectorField& UInst();
 
                 //- Return const reference to volume flux
-                const surfaceScalarField& phiInst() const;
+                inline const surfaceScalarField& phiInst() const;
 
                 //- Return reference to volume flux
-                surfaceScalarField& phiInst();
+                inline surfaceScalarField& phiInst();
 
 
             //- Return const reference to transport model
-            const singlePhaseTransportModel& laminarTransport() const;
+            inline const singlePhaseTransportModel& laminarTransport() const;
 
             //- Return reference to transport model
-            singlePhaseTransportModel& laminarTransport();
+            inline singlePhaseTransportModel& laminarTransport();
 
             //- Return const reference to the turbulence model
-            const autoPtr<incompressible::turbulenceModel>& turbulence() const;
+            inline const autoPtr<incompressible::turbulenceModel>&
+                turbulence() const;
 
             //- Return  reference to the turbulence model
-            autoPtr<incompressible::turbulenceModel>& turbulence();
+            inline autoPtr<incompressible::turbulenceModel>& turbulence();
 
             //- Return const reference to the turbulence model variables
-            const autoPtr<incompressible::RASModelVariables>&
+            inline const autoPtr<incompressible::RASModelVariables>&
                 RASModelVariables() const;
 
             //- Return reference to the turbulence model variables
-            autoPtr<incompressible::RASModelVariables>& RASModelVariables();
+            inline autoPtr<incompressible::RASModelVariables>&
+                RASModelVariables();
 
             //- Restore field values to the initial ones
             void restoreInitValues();
@@ -260,6 +262,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "incompressibleVarsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVarsI.H

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 PCOpt/NTUA
+    Copyright (C) 2021 FOSS GP
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::volScalarField& Foam::incompressibleVars::pInst() const
+{
+    return pPtr_();
+}
+
+
+inline Foam::volScalarField& Foam::incompressibleVars::pInst()
+{
+    return pPtr_();
+}
+
+
+inline const Foam::volVectorField& Foam::incompressibleVars::UInst() const
+{
+    return UPtr_();
+}
+
+
+inline Foam::volVectorField& Foam::incompressibleVars::UInst()
+{
+    return UPtr_();
+}
+
+
+inline const Foam::surfaceScalarField&
+Foam::incompressibleVars::phiInst() const
+{
+    return phiPtr_();
+}
+
+
+inline Foam::surfaceScalarField& Foam::incompressibleVars::phiInst()
+{
+    return phiPtr_();
+}
+
+
+inline const Foam::singlePhaseTransportModel&
+Foam::incompressibleVars::laminarTransport() const
+{
+    return laminarTransportPtr_();
+}
+
+
+inline Foam::singlePhaseTransportModel&
+Foam::incompressibleVars::laminarTransport()
+{
+    return laminarTransportPtr_();
+}
+
+
+inline const Foam::autoPtr<Foam::incompressible::turbulenceModel>&
+Foam::incompressibleVars::turbulence() const
+{
+    return turbulence_;
+}
+
+
+inline Foam::autoPtr<Foam::incompressible::turbulenceModel>&
+Foam::incompressibleVars::turbulence()
+{
+    return turbulence_;
+}
+
+
+inline const Foam::autoPtr<Foam::incompressible::RASModelVariables>&
+Foam::incompressibleVars::RASModelVariables() const
+{
+    return RASModelVariables_;
+}
+
+
+inline Foam::autoPtr<Foam::incompressible::RASModelVariables>&
+Foam::incompressibleVars::RASModelVariables()
+{
+    return RASModelVariables_;
+}
+
+
+// ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "incompressibleAdjointVars.H"
+#include "adjointBoundaryCondition.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,20 +67,6 @@ incompressibleAdjointVars::incompressibleAdjointVars
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const autoPtr<incompressibleAdjoint::adjointRASModel>&
-incompressibleAdjointVars::adjointTurbulence() const
-{
-    return adjointTurbulence_;
-}
-
-
-autoPtr<incompressibleAdjoint::adjointRASModel>&
-incompressibleAdjointVars::adjointTurbulence()
-{
-    return adjointTurbulence_;
-}
-
-
 void incompressibleAdjointVars::resetMeanFields()
 {
     if (solverControl_.average())
@@ -123,6 +110,33 @@ void incompressibleAdjointVars::nullify()
 }
 
 
+void incompressibleAdjointVars::updatePrimalBasedQuantities()
+{
+    /*
+    // WIP
+    for (fvPatchVectorField& pf : UaInst().boundaryFieldRef())
+    {
+        if (isA<adjointBoundaryCondition<vector>>(pf))
+        {
+            adjointBoundaryCondition<vector>& adjointBC =
+                refCast<adjointBoundaryCondition<vector>>(pf);
+            adjointBC.updatePrimalBasedQuantities();
+        }
+    }
+
+    for (fvPatchScalarField& pf : paInst().boundaryFieldRef())
+    {
+        if (isA<adjointBoundaryCondition<scalar>>(pf))
+        {
+            adjointBoundaryCondition<scalar>& adjointBC =
+                refCast<adjointBoundaryCondition<scalar>>(pf);
+            adjointBC.updatePrimalBasedQuantities();
+        }
+    }
+    */
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -104,11 +104,11 @@ public:
         // Access
 
             //- Return const reference to the adjointRASModel
-            const autoPtr<incompressibleAdjoint::adjointRASModel>&
+            inline const autoPtr<incompressibleAdjoint::adjointRASModel>&
                 adjointTurbulence() const;
 
             //- Return non-const reference to the adjointRASModel
-            autoPtr<incompressibleAdjoint::adjointRASModel>&
+            inline autoPtr<incompressibleAdjoint::adjointRASModel>&
                 adjointTurbulence();
 
             //- Reset mean fields to zero
@@ -119,6 +119,10 @@ public:
 
             //- Nullify all adjoint fields
             virtual void nullify();
+
+            //- Update primal based quantities of the adjoint boundary
+            //  conditions
+            virtual void updatePrimalBasedQuantities();
 };
 
 
@@ -128,6 +132,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "incompressibleAdjointVarsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVarsI.H

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 PCOpt/NTUA
+    Copyright (C) 2021 FOSS GP
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::autoPtr<Foam::incompressibleAdjoint::adjointRASModel>&
+Foam::incompressibleAdjointVars::adjointTurbulence() const
+{
+    return adjointTurbulence_;
+}
+
+
+inline Foam::autoPtr<Foam::incompressibleAdjoint::adjointRASModel>&
+Foam::incompressibleAdjointVars::adjointTurbulence()
+{
+    return adjointTurbulence_;
+}
+
+
+// ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -222,42 +222,6 @@ surfaceScalarField& incompressibleAdjointMeanFlowVars::phia()
 }
 
 
-const volScalarField& incompressibleAdjointMeanFlowVars::paInst() const
-{
-    return paPtr_();
-}
-
-
-volScalarField& incompressibleAdjointMeanFlowVars::paInst()
-{
-    return paPtr_();
-}
-
-
-const volVectorField& incompressibleAdjointMeanFlowVars::UaInst() const
-{
-    return UaPtr_();
-}
-
-
-volVectorField& incompressibleAdjointMeanFlowVars::UaInst()
-{
-    return UaPtr_();
-}
-
-
-const surfaceScalarField& incompressibleAdjointMeanFlowVars::phiaInst() const
-{
-    return phiaPtr_();
-}
-
-
-surfaceScalarField& incompressibleAdjointMeanFlowVars::phiaInst()
-{
-    return phiaPtr_();
-}
-
-
 bool incompressibleAdjointMeanFlowVars::computeMeanFields() const
 {
     return solverControl_.average();

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -158,22 +158,22 @@ public:
             // conditions should use these fields to execute a solver iteration
 
                 //- Return const reference to pressure
-                const volScalarField& paInst() const;
+                inline const volScalarField& paInst() const;
 
                 //- Return reference to pressure
-                volScalarField& paInst();
+                inline volScalarField& paInst();
 
                 //- Return const reference to velocity
-                const volVectorField& UaInst() const;
+                inline const volVectorField& UaInst() const;
 
                 //- Return reference to velocity
-                volVectorField& UaInst();
+                inline volVectorField& UaInst();
 
                 //- Return const reference to volume flux
-                const surfaceScalarField& phiaInst() const;
+                inline const surfaceScalarField& phiaInst() const;
 
                 //- Return reference to volume flux
-                surfaceScalarField& phiaInst();
+                inline surfaceScalarField& phiaInst();
 
             //- Return computeMeanFields bool
             bool computeMeanFields() const;
@@ -192,6 +192,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "incompressibleAdjointMeanFlowVarsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVarsI.H

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 PCOpt/NTUA
+    Copyright (C) 2021 FOSS GP
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::volScalarField&
+Foam::incompressibleAdjointMeanFlowVars::paInst() const
+{
+    return paPtr_();
+}
+
+
+inline Foam::volScalarField& Foam::incompressibleAdjointMeanFlowVars::paInst()
+{
+    return paPtr_();
+}
+
+
+inline const Foam::volVectorField&
+Foam::incompressibleAdjointMeanFlowVars::UaInst() const
+{
+    return UaPtr_();
+}
+
+
+inline Foam::volVectorField& Foam::incompressibleAdjointMeanFlowVars::UaInst()
+{
+    return UaPtr_();
+}
+
+
+inline const Foam::surfaceScalarField&
+Foam::incompressibleAdjointMeanFlowVars::phiaInst() const
+{
+    return phiaPtr_();
+}
+
+
+inline Foam::surfaceScalarField&
+Foam::incompressibleAdjointMeanFlowVars::phiaInst()
+{
+    return phiaPtr_();
+}
+
+
+
+// ************************************************************************* //

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -108,17 +108,15 @@ tmp<volScalarField> adjointSpalartAllmaras::r
     const volScalarField& Stilda
 ) const
 {
-    tmp<volScalarField> tr
+    auto tr =
-    (
+        tmp<volScalarField>::New
-        new volScalarField
         (
             min
             (
                 nuTilda()/(max(Stilda, minStilda_)*sqr(kappa_*y_)),
                 scalar(10)
             )
-        )
+        );
-    );
     tr.ref().boundaryFieldRef() == Zero;
 
     return tr;
@@ -138,10 +136,9 @@ tmp<volScalarField> adjointSpalartAllmaras::fw
 
 tmp<volScalarField> adjointSpalartAllmaras::DnuTildaEff() const
 {
-    return tmp<volScalarField>
+    return
-    (
+        tmp<volScalarField>::New
-        new volScalarField("DnuTildaEff", (nuTilda() + this->nu())/sigmaNut_)
+            ("DnuTildaEff", (nuTilda() + this->nu())/sigmaNut_);
-    );
 }
 
 
@@ -166,7 +163,7 @@ tmp<volScalarField> adjointSpalartAllmaras::dFv1_dChi
 {
     volScalarField chi3(pow3(chi));
 
-    return 3.0*pow3(Cv1_)*sqr(chi/(chi3+pow3(Cv1_)));
+    return 3.0*pow3(Cv1_)*sqr(chi/(chi3 + pow3(Cv1_)));
 }
 
 
@@ -246,8 +243,8 @@ tmp<volScalarField> adjointSpalartAllmaras::dr_dStilda
 {
     tmp<volScalarField> tdrdStilda
     (
-        - nuTilda()/sqr(max(Stilda, minStilda_)*kappa_*y_)
+      - nuTilda()/sqr(max(Stilda, minStilda_)*kappa_*y_)
-        *(scalar(10) - r_)/(scalar(10) - r_ + SMALL)
+       *(scalar(10) - r_)/(scalar(10) - r_ + SMALL)
     );
     tdrdStilda.ref().boundaryFieldRef() == Zero;
 
@@ -262,8 +259,8 @@ tmp<volScalarField> adjointSpalartAllmaras::dr_dDelta
 {
     tmp<volScalarField> tdrdDelta
     (
-        -2.*nuTilda()/(max(Stilda, minStilda_)*sqr(kappa_*y_)*y_)
+      - 2.*nuTilda()/(max(Stilda, minStilda_)*sqr(kappa_*y_)*y_)
-        *(scalar(10) - r_)/(scalar(10) - r_ + SMALL)
+       *(scalar(10) - r_)/(scalar(10) - r_ + SMALL)
     );
     tdrdDelta.ref().boundaryFieldRef() == Zero;
 
@@ -281,9 +278,10 @@ tmp<volScalarField> adjointSpalartAllmaras::dfw_dr
     dimensionedScalar pow6Cw3 = pow6(Cw3_);
     volScalarField pow6g(pow6(g));
 
-    return  pow6Cw3/(pow6g + pow6Cw3)
+    return
-           *pow((1.0 + pow6Cw3)/(pow6g + pow6Cw3), 1.0/6.0)
+        pow6Cw3/(pow6g + pow6Cw3)
-           *(1.0 + Cw2_*(6.0*pow5(r_) - 1.0));
+       *pow((1.0 + pow6Cw3)/(pow6g + pow6Cw3), 1.0/6.0)
+       *(1.0 + Cw2_*(6.0*pow5(r_) - 1.0));
 }
 
 
@@ -361,9 +359,9 @@ tmp<volVectorField> adjointSpalartAllmaras::conservativeMomentumSource()
         const fvPatch& patch = mesh_.boundary()[pI];
         if(!isA<coupledFvPatch>(patch))
         {
-            vectorField nf(patch.nf());
+            tmp<vectorField> tnf = patch.nf();
             adjMomentumBCSourcePtr_()[pI] =
-                (nf & momentumSourceMult_.boundaryField()[pI])
+                (tnf & momentumSourceMult_.boundaryField()[pI])
                *nuaTilda().boundaryField()[pI];
         }
     }
@@ -745,6 +743,8 @@ tmp<fvVectorMatrix> adjointSpalartAllmaras::divDevReff(volVectorField& Ua) const
 
 tmp<volVectorField> adjointSpalartAllmaras::adjointMeanFlowSource()
 {
+    addProfiling
+        (adjointSpalartAllmaras, "adjointSpalartAllmaras::addMomentumSource");
     // cm formulation
     //return (- nuTilda()*fvc::grad(nuaTilda() - conservativeMomentumSource());
 
@@ -774,7 +774,7 @@ tmp<scalarField> adjointSpalartAllmaras::diffusionCoeffVar1(label patchI) const
 
     diffCoeff =
         (nuTilda().boundaryField()[patchI] + nu()().boundaryField()[patchI])
-        /sigmaNut_.value();
+       /sigmaNut_.value();
 
     return tdiffCoeff;
 }
@@ -796,14 +796,13 @@ const boundaryVectorField& adjointSpalartAllmaras::wallShapeSensitivities()
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
 
-        tmp<vectorField> tnf(patch.nf());
+        tmp<vectorField> tnf = patch.nf();
-        const vectorField& nf = tnf();
         if (isA<wallFvPatch>(patch) && patch.size() != 0)
         {
             wallShapeSens[patchI] =
               - nuaTilda().boundaryField()[patchI].snGrad()
-              * diffusionCoeffVar1(patchI)()
+               *diffusionCoeffVar1(patchI)
-              * nuTilda().boundaryField()[patchI].snGrad() * nf;
+               *nuTilda().boundaryField()[patchI].snGrad()*tnf;
         }
     }
 
@@ -818,11 +817,10 @@ const boundaryVectorField& adjointSpalartAllmaras::wallFloCoSensitivities()
     forAll(mesh_.boundary(), patchI)
     {
         tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
-        const vectorField& nf = tnf();
 
         wallFloCoSens[patchI] =
             nuaTilda().boundaryField()[patchI]
-          * nuTilda().boundaryField()[patchI] * nf;
+           *nuTilda().boundaryField()[patchI]*tnf;
     }
 
     return wallFloCoSens;
@@ -848,18 +846,15 @@ tmp<volScalarField> adjointSpalartAllmaras::distanceSensitivities()
     volScalarField dfw_dDelta
         (this->dfw_dDelta(Stilda_, dfw_dr, dStilda_dDelta));
 
-
+    auto tadjointEikonalSource =
-    tmp<volScalarField> tadjointEikonalSource
+        tmp<volScalarField>::New
-    (
-        new volScalarField
         (
             "adjointEikonalSource" + type(),
             (
-                - Cb1_*nuTilda()*dStilda_dDelta
+              - Cb1_*nuTilda()*dStilda_dDelta
-                + Cw1_*sqr(nuTilda()/y_)*(dfw_dDelta - 2.*fw_/y_)
+              + Cw1_*sqr(nuTilda()/y_)*(dfw_dDelta - 2.*fw_/y_)
             )*nuaTilda()
-        )
+        );
-    );
     volScalarField& adjointEikonalSource = tadjointEikonalSource.ref();
 
     // if wall functions are used, add appropriate source terms
@@ -883,8 +878,8 @@ tmp<volScalarField> adjointSpalartAllmaras::distanceSensitivities()
         {
             const scalar kappa_(0.41);
             const scalar E_(9.8);
-            const tmp<vectorField> tnf(patch.nf());
+            tmp<vectorField> tnf = patch.nf();
-            const vectorField& nf = tnf();
+            const vectorField& nf = tnf.ref();
             const scalarField& magSf = patch.magSf();
 
             const fvPatchVectorField& Up = U.boundaryField()[patchi];
@@ -936,8 +931,8 @@ tmp<volScalarField> adjointSpalartAllmaras::distanceSensitivities()
                 label cellI = faceCells[faceI];
                 adjointEikonalSource[cellI] -=
                     2.*( rt[faceI] + Uap_t[faceI] )
-                  * vtau[faceI]*Cwf_d[faceI]*magSf[faceI]
+                   *vtau[faceI]*Cwf_d[faceI]*magSf[faceI]
-                  / V[cellI]; // Divide with cell volume since the term
+                   /V[cellI]; // Divide with cell volume since the term
                               // will be used as a source term in the
                               // adjoint eikonal equation
             }
@@ -952,9 +947,13 @@ tmp<volTensorField> adjointSpalartAllmaras::FISensitivityTerm()
 {
     const volVectorField& U  = primalVars_.U();
 
-    volTensorField gradU(fvc::grad(U));
+    tmp<volTensorField> tgradU = fvc::grad(U);
-    volVectorField gradNuTilda(fvc::grad(nuTilda()));
+    const volTensorField& gradU = tgradU.cref();
-    volVectorField gradNuaTilda(fvc::grad(nuaTilda()));
+    tmp<volVectorField> tgradNuTilda = fvc::grad(nuTilda());
+    volVectorField& gradNuTilda = tgradNuTilda.ref();
+    tmp<volVectorField> tgradNuaTilda = fvc::grad(nuaTilda());
+    volVectorField::Boundary& gradNuaTildab =
+        tgradNuaTilda.ref().boundaryFieldRef();
 
     // Explicitly correct the boundary gradient to get rid of
     // the tangential component
@@ -963,7 +962,7 @@ tmp<volTensorField> adjointSpalartAllmaras::FISensitivityTerm()
         const fvPatch& patch = mesh_.boundary()[patchI];
         if (isA<wallFvPatch>(patch))
         {
-            tmp<vectorField> tnf(patch.nf());
+            tmp<vectorField> tnf = patch.nf();
             const vectorField& nf = tnf();
             // gradU:: can cause problems in zeroGradient patches for U
             // and zero fixedValue for nuTilda.
@@ -971,9 +970,9 @@ tmp<volTensorField> adjointSpalartAllmaras::FISensitivityTerm()
             //gradU.boundaryField()[patchI] =
             //  nf * U_.boundaryField()[patchI].snGrad();
             gradNuTilda.boundaryFieldRef()[patchI]  =
-                nf * nuTilda().boundaryField()[patchI].snGrad();
+                nf*nuTilda().boundaryField()[patchI].snGrad();
-            gradNuaTilda.boundaryFieldRef()[patchI] =
+            gradNuaTildab[patchI] =
-                nf * nuaTilda().boundaryField()[patchI].snGrad();
+                nf*nuaTilda().boundaryField()[patchI].snGrad();
         }
     }
 
@@ -987,6 +986,7 @@ tmp<volTensorField> adjointSpalartAllmaras::FISensitivityTerm()
         )
        /(Omega + dimensionedScalar("SMALL", Omega.dimensions(), SMALL))
     );
+    tgradU.clear();
 
     volScalarField chi(this->chi());
     volScalarField fv1(this->fv1(chi));
@@ -997,29 +997,24 @@ tmp<volTensorField> adjointSpalartAllmaras::FISensitivityTerm()
     volScalarField dfw_dOmega
         (this->dfw_dOmega(Stilda_, dfw_dr, dStilda_dOmega));
 
-    // Assemply of the return field
+    return
-    tmp<volTensorField> tvolSensTerm
+        tmp<volTensorField>::New
-    (
-        new volTensorField
         (
             "volSensTerm",
             // jk, cm formulation for the TM model convection
-            - (nuaTilda() * (U * gradNuTilda))
+          - (nuaTilda()*(U*gradNuTilda))
             // jk, symmetric in theory
-            + nuaTilda()*fvc::grad(DnuTildaEff() * gradNuTilda)().T()
+          + nuaTilda()*fvc::grad(DnuTildaEff()*gradNuTilda)().T()
             // jk
-            - DnuTildaEff() * (gradNuaTilda * gradNuTilda)
+          - DnuTildaEff()*(tgradNuaTilda*gradNuTilda)
             // symmetric
-            + 2.*nuaTilda()*Cb2_/sigmaNut_ * (gradNuTilda * gradNuTilda)
+          + 2.*nuaTilda()*Cb2_/sigmaNut_*(gradNuTilda*gradNuTilda)
-            + (
+          + (
-                - Cb1_*nuTilda()*dStilda_dOmega
+              - Cb1_*nuTilda()*dStilda_dOmega
-                + Cw1_*sqr(nuTilda()/y_)*dfw_dOmega
+              + Cw1_*sqr(nuTilda()/y_)*dfw_dOmega
-              )
+            )
-            * nuaTilda() * deltaOmega // jk
+           *nuaTilda()*deltaOmega // jk
-         )
+        );
-    );
-
-    return tvolSensTerm;
 }
 
 
@@ -1031,6 +1026,8 @@ void adjointSpalartAllmaras::nullify()
 
 void adjointSpalartAllmaras::correct()
 {
+    addProfiling
+        (adjointSpalartAllmaras, "adjointSpalartAllmaras::correct");
     if (!adjointTurbulence_)
     {
         return;
@@ -1046,7 +1043,7 @@ void adjointSpalartAllmaras::correct()
     volScalarField gradNua(gradNuTilda_ & fvc::grad(nuaTilda()));
     volScalarField gradUaR
     (
-        2.0*fvc::grad(Ua,"adjointProductionUa") && symmAdjointProductionU_
+        2.0*fvc::grad(Ua) && symmAdjointProductionU_
     );
 
     dimensionedScalar oneOverSigmaNut = 1./sigmaNut_;

src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointTurbulenceModel/adjointTurbulenceModel.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -183,6 +183,23 @@ public:
             //return primalVars_.turbulence()().nuEff();
         }
 
+        //- Return the effective viscosity on a given patch
+        virtual tmp<scalarField> nuEff(const label patchI) const
+        {
+            // Go through RASModelVariables::nutRef in order to obtain
+            // the mean field, if present
+            const singlePhaseTransportModel& lamTrans =
+                primalVars_.laminarTransport();
+            const autoPtr<incompressible::RASModelVariables>&
+                turbVars = primalVars_.RASModelVariables();
+
+            return
+                (
+                    lamTrans.nu()().boundaryField()[patchI]
+                  + turbVars().nutRef().boundaryField()[patchI]
+                );
+        }
+
         //- Return the effective stress tensor including the laminar stress
         virtual tmp<volSymmTensorField> devReff() const = 0;