diff --git a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.C b/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.C
index 2be1febb3505e72a515f712a073a96dfb7e458b3..057996e18ab4f7e18f812bad73794cbf1f1da071 100644
--- a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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
diff --git a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.H b/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.H
index d80c75b12eb43bb4cb62649d895b5fbf10dbbb0d..7738ef491c5d63f1e530d6b2e6b8fcf407ea905b 100644
--- a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCModel/ATCModel.H
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
 };
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCUaGradU/ATCUaGradU.C b/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCUaGradU/ATCUaGradU.C
index 99d73932ff4a06fb78dddae2919c059cfbbcba18..660a697ee2f0808471cfffc550973a7e50938ec4 100644
--- a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCUaGradU/ATCUaGradU.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.C b/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.C
index 214c3aefbcfbfce67bba88d9557e4b971d74d6fd..ec389d20d75c8ba292708883afff79d7bd624772 100644
--- a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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,7 +101,7 @@ void ATCstandard::addATC(fvVectorMatrix& UaEqn)
     smoothATC();
 
     // actual ATC term
-    UaEqn += fvm::Su(ATC_, Ua);
+    UaEqn += ATC_.internalField();
 }
 
 
@@ -149,6 +153,24 @@ tmp<volTensorField> ATCstandard::getFISensitivityTerm() const
 }
 
 
+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_)
+    {
+        UForATC.reset(new volVectorField(fvc::reconstruct(phi)));
+    }
+    else
+    {
+        UForATC.reset(new volVectorField(U));
+    }
+    gradU_ = fvc::grad(UForATC(), "gradUATC");
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.H b/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.H
index 4c65a9c280331162c5947d89ab8fbee6a9c187f0..0d7f4eb99f41d55c1a338ce4e5682806e2c3a36d 100644
--- a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.H
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
 };
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C
index 1d87524b235e64cfa3a9f335e9653f401ad1ebae..0e49b2815e3cf3e786e892154bb04b02578cfe72 100644
--- a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.C
+++ b/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) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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
-    (
-        surfaceInterpolationScheme<Type2>::New
-        (
-            mesh,
-            mesh.interpolationScheme("interpolate(" + name + ")")
-        )
-    );
+    // If there are many patches calling this function, the computation of
+    // the surface field might be a significant computational
+    // burden. Cache the interpolated field and fetch from the registry when
+    // possible.
+    const word surfFieldName("interpolated" + name + "ForBoundaryGrad");
+    typedef GeometricField<Type2, fvsPatchField, surfaceMesh> surfFieldType;
+    surfFieldType* surfFieldPtr =
+        mesh.objectRegistry::template getObjectPtr<surfFieldType>(surfFieldName);
 
-    GeometricField<Type2, fvsPatchField, surfaceMesh> surfField
-    (
-        tinterpScheme().interpolate(field)
-    );
+    if (!surfFieldPtr)
+    {
+        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
 <
diff --git a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.H b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.H
index 4ec46c64c96867bc520b24b39cbfe082eb8c8933..dd2f92c98d32ba174fa1a9c51378872da67cdaf2 100644
--- a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryCondition.H
+++ b/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) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
 };
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.C b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.C
index 8abaf44d00c9eadc46ee9e0168c4620201974ef3..8b70903941d87ba536188bfe5146a0fb85e0ae75 100644
--- a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.C
+++ b/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) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
diff --git a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.H b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.H
index 3201b245ab79d4d7c403ef690b1944635b0b4b1a..5dd565f9d730a96168ec5558ac4cd19ca95de75c 100644
--- a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvPatchVectorField.H
+++ b/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) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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:
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointWallVelocity/adjointWallVelocityFvPatchVectorField.C b/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointWallVelocity/adjointWallVelocityFvPatchVectorField.C
index 32372f99affc5c919f480a9d14400b602efd258e..c26a7b863a7a799aa1e735fd1e8ed3a576069e97 100644
--- a/src/optimisation/adjointOptimisation/adjoint/adjointBoundaryConditions/adjointWallVelocity/adjointWallVelocityFvPatchVectorField.C
+++ b/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) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C
index 6f108478a0fd0e157079c4d09687e19ebd727554..f694aa7f1ee0ceeba64d0276d1f01d2ed458c3b3 100644
--- a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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,8 +124,7 @@ tmp<scalarField> boundaryAdjointContributionIncompressible::pressureSource()
     const vectorField& adjointTurbulenceContr =
         adjointRAS().adjointMomentumBCSource()[patch_.index()];
 
-    tmp<vectorField> tnf = patch_.nf();
-    const vectorField& nf = tnf();
+    tmp<vectorField> nf = patch_.nf();
 
     source += adjointTurbulenceContr & nf;
 
@@ -139,11 +136,10 @@ tmp<vectorField>
 boundaryAdjointContributionIncompressible::tangentVelocitySource()
 {
     // Objective function contribution
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
     tmp<vectorField> tsource =
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdvt
         );
 
@@ -167,73 +163,58 @@ boundaryAdjointContributionIncompressible::tangentVelocitySource()
 tmp<vectorField>
 boundaryAdjointContributionIncompressible::normalVelocitySource()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
-    tmp<vectorField> tsource =
+    return
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdp
         );
 
-    return (tsource);
 }
 
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::energySource()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
-    tmp<scalarField> tsource =
+    return
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdT
         );
 
-    return (tsource);
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::adjointTMVariable1Source()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
-    tmp<scalarField> tsource =
+    return
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdTMvar1
         );
 
-    return (tsource);
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::adjointTMVariable2Source()
 {
-    PtrList<objective>& objectives = objectiveManager_.getObjectiveFunctions();
-    tmp<scalarField> tsource =
+    return
         sumContributions
         (
-            objectives,
+            objectiveManager_.getObjectiveFunctions(),
             &objectiveIncompressible::boundarydJdTMvar2
         );
 
-    return (tsource);
 }
 
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::momentumDiffusion()
 {
-    tmp<scalarField> tnuEff(new scalarField(patch_.size(), Zero));
-    scalarField& nuEff = tnuEff.ref();
-
-    const autoPtr<incompressibleAdjoint::adjointRASModel>&
-        adjointTurbulenceModel = adjointVars().adjointTurbulence();
 
-    nuEff = adjointTurbulenceModel().nuEff()().boundaryField()[patch_.index()];
-
-    return tnuEff;
+    return adjointVars().adjointTurbulence()().nuEff(patch_.index());
 }
 
 
@@ -272,64 +253,40 @@ tmp<scalarField> boundaryAdjointContributionIncompressible::thermalDiffusion()
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::wallDistance()
 {
-    tmp<scalarField> twallDist(new scalarField(patch_.size(), Zero));
-    scalarField& wallDist = twallDist.ref();
-
-    wallDist = primalVars_.turbulence()->y()[patch_.index()];
-
-    return twallDist;
+    return primalVars_.turbulence()->y()[patch_.index()];
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::TMVariable1Diffusion()
 {
-    const autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS =
-        adjointVars().adjointTurbulence();
-
-    tmp<scalarField> tdiffCoeff =
-        adjointRAS().diffusionCoeffVar1(patch_.index());
+    return
+        adjointVars().adjointTurbulence()->diffusionCoeffVar1(patch_.index());
 
-    return tdiffCoeff;
 }
 
 
 tmp<scalarField>
 boundaryAdjointContributionIncompressible::TMVariable2Diffusion()
 {
-    const autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS =
-        adjointVars().adjointTurbulence();
-
-    tmp<scalarField> tdiffCoeff =
-        adjointRAS().diffusionCoeffVar2(patch_.index());
-
-    return tdiffCoeff;
+    return
+        adjointVars().adjointTurbulence()->diffusionCoeffVar2(patch_.index());
 }
 
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::TMVariable1()
 {
-    const autoPtr<incompressible::RASModelVariables>& RASVariables =
-        primalVars_.RASModelVariables();
-    tmp<scalarField> tboundField(new scalarField(patch_.size(), Zero));
-    scalarField& boundField = tboundField.ref();
-
-    boundField = RASVariables().TMVar1().boundaryField()[patch_.index()];
-
-    return tboundField;
+    return
+        primalVars_.RASModelVariables()->TMVar1().
+            boundaryField()[patch_.index()];
 }
 
 
 tmp<scalarField> boundaryAdjointContributionIncompressible::TMVariable2()
 {
-    const autoPtr<incompressible::RASModelVariables>& RASVariables =
-        primalVars_.RASModelVariables();
-    tmp<scalarField> tboundField(new scalarField(patch_.size(), Zero));
-    scalarField& boundField = tboundField.ref();
-
-    boundField = RASVariables().TMVar2().boundaryField()[patch_.index()];
-
-    return tboundField;
+    return
+        primalVars_.RASModelVariables()->TMVar2().
+            boundaryField()[patch_.index()];
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/adjointSolver/adjointSolver.C b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/adjointSolver/adjointSolver.C
index aff18976f3ea6c9621abcd31f8933a78c35e3ea0..e212b596dc8975f9de1af1eaebb9add2fceb6447 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/adjointSolver/adjointSolver.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C
index 064c176ee7ffd5d1278bfc3df5f8146b095ff9ea..5fa4603e8c97e8213ed868ca775c885beae3f803 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C
+++ b/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) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.C b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.C
index b158a1310bce49c288983917f493ac4579b46652..2ffa2d1b0371aa40d8e920d0d68944e4676b5611 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.C
+++ b/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) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
     }
 }
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C
index 45fb8baffe6ac61ad8cd11c9c6973c21175c660c..0fe83537d5ab48643484bb68763264b38923fe38 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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())
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.H b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.H
index 1ddc919cb93190756bbf84b07cc19259a2199314..38c34a8484f1987fc9d221d5b1c8af6fc231d380 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.H
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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
 
 // ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVarsI.H b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVarsI.H
new file mode 100644
index 0000000000000000000000000000000000000000..c1777b45796ea6b89da6bc65c5358c51b9568afb
--- /dev/null
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVarsI.H
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.C b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.C
index c544f74d9005b611189859386c8cba2e5ec1cf94..6bf4cbc2dfc4a9001b97ceb1ef10f6c57dea7b41 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.H b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.H
index 815edfb43b833c2ede1ae277c20dec59d8951ff0..8bb60a23ebfcf47d96dded9aa4f11ffad774dbc1 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.H
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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
 
 // ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVarsI.H b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVarsI.H
new file mode 100644
index 0000000000000000000000000000000000000000..75576a43c3c00fdf7d3005fde356410e61cdea07
--- /dev/null
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVarsI.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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_;
+}
+
+
+
+
+
+// ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C
index eab4d49e5e7aa9730f85e97a76688531248df5ff..7a23aba9254ffa9f3bbc863822bea61e6b39398c 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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();
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.H b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.H
index 5082cf6f525bb183ee1892df2518f693cf906519..9e523c61e2e86ed9af73114aeae26de363bdba04 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVars.H
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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
 
 // ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVarsI.H b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVarsI.H
new file mode 100644
index 0000000000000000000000000000000000000000..e70f29030f7c2e6a6ad6b6e1c30d58d73e951a6a
--- /dev/null
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjointMeanFlow/incompressibleAdjointMeanFlowVarsI.H
@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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_();
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C
index 7878349822708ebf00a5f0326e72597cb5e38974..edb83dded5ac0fccf81fe94b0ed42a417c57b2e0 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -743,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());
 
@@ -1029,6 +1031,8 @@ void adjointSpalartAllmaras::nullify()
 
 void adjointSpalartAllmaras::correct()
 {
+    addProfiling
+        (adjointSpalartAllmaras, "adjointSpalartAllmaras::correct");
     if (!adjointTurbulence_)
     {
         return;
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointTurbulenceModel/adjointTurbulenceModel.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointTurbulenceModel/adjointTurbulenceModel.H
index e7039465e53fee1835ff04b9a2c9d1490455ddb8..04adc50c37e51e887e3e06dbb9ab3c7715b6a5a6 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointTurbulenceModel/adjointTurbulenceModel.H
+++ b/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) 2013-2019 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    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;