ENH: changes reducing the peak memory consumption of shape sensitivities

The multiplier of grad(dxdb) is a volTensorField which, by itself, is
memory consuming. The function computing it though was sloppy in terms
of memory management, constituting the peak memory consumption during an
adjoint optimisation. Initial changes to remedy the problem include the
deallocation of some of the volTensorFields included in the computation
of grad(dxdb) once unneeded, the utilisation of volSymmTensorFields
instead of volTensorFields where possible and avoiding allocating some
unnecessary intermediate fields.

Actions to further reduce memory consumption:

- For historical reasons, the code computes/stores the transpose of
  grad(dxdb), which is then transposed when used in the computation of
  the FI or the ESI sensitivity derivatives.  This redundant
  transposition can be avoid, saving the allocation of an additional
  volTensorField, but the changes need to permeate a number of places in
  the code that contribute to grad(dxdb) (e.g.  ATC, adjoint turbulence
  models, adjoint MRF, etc).
- Allocation of unnecessary pointers in the objective class should be
  avoided.

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.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
@@ -133,8 +133,11 @@ void FIBase::accumulateIntegrand(const scalar dt)
     PtrList<objective>& functions(objectiveManager_.getObjectiveFunctions());
     for (objective& func : functions)
     {
-        divDxDbMult_ +=
-            func.weight()*func.divDxDbMultiplier().primitiveField()*dt;
+        if (func.hasDivDxDbMult())
+        {
+            divDxDbMult_ +=
+                func.weight()*func.divDxDbMultiplier().primitiveField()*dt;
+        }
     }
 
     // Terms from fvOptions
@@ -154,7 +157,6 @@ void FIBase::accumulateIntegrand(const scalar dt)
 
     // Accumulate sensitivities due to boundary conditions
     accumulateBCSensitivityIntegrand(dt);
-
 }
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C

@@ -602,61 +602,6 @@ void sensitivitySurface::accumulateIntegrand(const scalar dt)
     autoPtr<incompressibleAdjoint::adjointRASModel>& adjointTurbulence =
         adjointVars_.adjointTurbulence();
 
-    Info<< "    Calculating auxiliary quantities " << endl;
-    // Fields needed to calculate adjoint sensitivities
-    const autoPtr<incompressible::RASModelVariables>&
-       turbVars = primalVars_.RASModelVariables();
-    const singlePhaseTransportModel& lamTransp = primalVars_.laminarTransport();
-    volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
-    volTensorField gradUa(fvc::grad(Ua));
-    volTensorField gradU(fvc::grad(U));
-
-    // Explicitly correct the boundary gradient to get rid of 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();
-            gradU.boundaryFieldRef()[patchI] =
-                nf*U.boundaryField()[patchI].snGrad();
-        }
-    }
-
-    // Auxiliary terms
-    volVectorField gradp(fvc::grad(p));
-    volTensorField stress(nuEff*(gradU + T(gradU)));
-    autoPtr<volVectorField> stressXPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressX", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressYPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressY", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressZPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressZ", stress.dimensions())
-    );
-
-    stressXPtr().replace(0, stress.component(0));
-    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()));
-
     // Accumulate source for additional post-processing PDEs, if necessary
     if (includeDistance_)
     {
@@ -675,8 +620,105 @@ void sensitivitySurface::accumulateIntegrand(const scalar dt)
     DebugInfo
         << "    Calculating adjoint sensitivity. " << endl;
 
+    tmp<volScalarField> tnuEff = adjointTurbulence->nuEff();
+    const volScalarField& nuEff = tnuEff.ref();
+
     // Sensitivities do not include locale surface area by default.
-    // Part of the sensitivities that multiplies dxFace/db
+    // The part of the sensitivities that multiplies dxFace/db follows
+
+    // Deal with the stress part first since it's the most awkward in terms
+    // of memory managment
+    if (includeGradStressTerm_)
+    {
+        // Terms corresponding to contributions from converting delta
+        // to thetas are added through the corresponding adjoint
+        // boundary conditions instead of grabbing contributions from
+        // the objective function.  Useful to have a unified
+        // formulation for low- and high-re meshes
+
+        tmp<volVectorField> tgradp = fvc::grad(p);
+        const volVectorField& gradp = tgradp.ref();
+        for (const label patchI : sensitivityPatchIDs_)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            tmp<vectorField> tnf = patch.nf();
+            const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
+            wallFaceSensVecPtr_()[patchI] -=
+                (Uab & tnf)*gradp.boundaryField()[patchI]*dt;
+        }
+        tgradp.clear();
+
+        // 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
+        forAll(mesh_.boundary(), patchI)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            if (isA<wallFvPatch>(patch))
+            {
+                tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
+                gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad();
+            }
+        }
+
+        tmp<volSymmTensorField> tstress = nuEff*twoSymm(tgradU);
+        const volSymmTensorField& stress = tstress.cref();
+        autoPtr<volVectorField> ptemp
+            (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 : sensitivityPatchIDs_)
+            {
+                const fvPatch& patch = mesh_.boundary()[patchI];
+                tmp<vectorField> tnf = patch.nf();
+                const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
+                wallFaceSensVecPtr_()[patchI] +=
+                    (
+                        Uab.component(idir)
+                       *(gradStressDir.boundaryField()[patchI] & tnf)
+                    )*dt;
+            }
+        }
+    }
+
+    // Transpose part of the adjoint stresses
+    // Dealt with separately to deallocate gradUa as soon as possible
+    if (includeTransposeStresses_)
+    {
+        tmp<volTensorField> tgradUa = fvc::grad(Ua);
+        const volTensorField::Boundary& gradUabf =
+            tgradUa.ref().boundaryField();
+
+        for (const label patchI : sensitivityPatchIDs_)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            tmp<vectorField> tnf = patch.nf();
+            const vectorField& nf = tnf();
+            vectorField gradUaNf
+                (
+                    useSnGradInTranposeStresses_
+                  ? (Ua.boundaryField()[patchI].snGrad() & nf)*nf
+                  : (gradUabf[patchI] & nf)
+                );
+            wallFaceSensVecPtr_()[patchI] -=
+                nuEff.boundaryField()[patchI]
+               *(gradUaNf & U.boundaryField()[patchI].snGrad())*nf;
+        }
+    }
+
     for (const label patchI : sensitivityPatchIDs_)
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
@@ -706,51 +748,17 @@ void sensitivitySurface::accumulateIntegrand(const scalar dt)
           * nf
         );
 
-
-        if (includeTransposeStresses_)
-        {
-            vectorField gradUaNf
-                (
-                    useSnGradInTranposeStresses_ ?
-                    (Ua.boundaryField()[patchI].snGrad() & nf)*nf :
-                    (gradUa.boundaryField()[patchI] & nf)
-                );
-
-            stressTerm -=
-                nuEff.boundaryField()[patchI]
-               *(gradUaNf & U.boundaryField()[patchI].snGrad())
-               *nf;
-        }
-
         if (includeDivTerm_)
         {
             stressTerm +=
                 scalar(1./3.)*nuEff.boundaryField()[patchI]
               * (
                     ((Ua.boundaryField()[patchI].snGrad() &nf)*nf)
-                    & U.boundaryField()[patchI].snGrad()
+                  & U.boundaryField()[patchI].snGrad()
                 )
               * nf;
         }
 
-        vectorField gradStressTerm(patch.size(), Zero);
-        if (includeGradStressTerm_)
-        {
-            // Terms corresponding to contributions from converting delta to
-            // thetas are added through the corresponding adjoint boundary
-            // conditions instead of grabbing contributions from the objective
-            // function.  Useful to have a unified formulation for low- and
-            // high-re meshes
-            const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
-            gradStressTerm = - ((Uab & nf)*gradp.boundaryField()[patchI]);
-            gradStressTerm +=
-            (
-                Uab.component(0) * gradStressX.boundaryField()[patchI]
-              + Uab.component(1) * gradStressY.boundaryField()[patchI]
-              + Uab.component(2) * gradStressZ.boundaryField()[patchI]
-            ) & nf;
-        }
-
         // Adjoint pressure terms
         vectorField pressureTerm(patch.size(), Zero);
         if (includePressureTerm_)
@@ -783,7 +791,6 @@ void sensitivitySurface::accumulateIntegrand(const scalar dt)
         wallFaceSensVecPtr_()[patchI] +=
         (
             stressTerm
-          + gradStressTerm
           + pressureTerm
           + adjointTMsensitivities[patchI]
           + dxdbMultiplierTot

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.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-2020, 2022 PCOpt/NTUA
+    Copyright (C) 2013-2020, 2022 FOSS GP
     Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -415,67 +415,12 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt)
     autoPtr<incompressibleAdjoint::adjointRASModel>& adjointTurbulence =
         adjointVars_.adjointTurbulence();
 
-    DebugInfo
-        << "    Calculating auxilary quantities " << endl;
-
-    // Fields needed to calculate adjoint sensitivities
-    volScalarField nuEff(adjointTurbulence->nuEff());
-    volTensorField gradUa(fvc::grad(Ua));
-    volTensorField gradU(fvc::grad(U));
-
-    // Explicitly correct the boundary gradient to get rid of 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();
-            gradU.boundaryFieldRef()[patchI] =
-                nf*U.boundaryField()[patchI].snGrad();
-        }
-    }
-
-    // Auxiliary terms
-    volVectorField gradp(fvc::grad(p));
-    volTensorField stress(nuEff*(gradU + T(gradU)));
-    autoPtr<volVectorField> stressXPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressX", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressYPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressY", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressZPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressZ", stress.dimensions())
-    );
-
-    stressXPtr().replace(0, stress.component(0));
-    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()));
-
     // Solve extra equations if necessary
     if (includeDistance_)
     {
         eikonalSolver_->accumulateIntegrand(dt);
     }
 
-    autoPtr<boundaryVectorField> meshMovementSensPtr(nullptr);
     if (includeMeshMovement_)
     {
         meshMovementSolver_->accumulateIntegrand(dt);
@@ -491,6 +436,101 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt)
     DebugInfo
         << "    Calculating adjoint sensitivity. " << endl;
 
+    tmp<volScalarField> tnuEff = adjointTurbulence->nuEff();
+    const volScalarField& nuEff = tnuEff.ref();
+
+    // Deal with the stress part first since it's the most awkward in terms
+    // of memory managment
+    if (includeGradStressTerm_)
+    {
+        // Terms corresponding to contributions from converting delta
+        // to thetas are added through the corresponding adjoint
+        // boundary conditions instead of grabbing contributions from
+        // the objective function.  Useful to have a unified
+        // formulation for low- and high-re meshes
+
+        tmp<volVectorField> tgradp = fvc::grad(p);
+        const volVectorField& gradp = tgradp.ref();
+        for (const label patchI : sensitivityPatchIDs_)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            tmp<vectorField> tnf = patch.nf();
+            const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
+            wallFaceSens_()[patchI] -=
+                (Uab & tnf)*gradp.boundaryField()[patchI]*dt;
+        }
+        tgradp.clear();
+
+        // 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
+        forAll(mesh_.boundary(), patchI)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            if (isA<wallFvPatch>(patch))
+            {
+                tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
+                gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad();
+            }
+        }
+
+        tmp<volSymmTensorField> tstress = nuEff*twoSymm(tgradU);
+        const volSymmTensorField& stress = tstress.cref();
+        autoPtr<volVectorField> ptemp
+            (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 : sensitivityPatchIDs_)
+            {
+                const fvPatch& patch = mesh_.boundary()[patchI];
+                tmp<vectorField> tnf = patch.nf();
+                const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
+                wallFaceSens_()[patchI] +=
+                    (
+                        Uab.component(idir)
+                       *(gradStressDir.boundaryField()[patchI] & tnf)
+                    )*dt;
+            }
+        }
+    }
+
+    // Transpose part of the adjoint stresses
+    // Dealt with separately to deallocate gradUa as soon as possible
+    if (includeTransposeStresses_)
+    {
+        tmp<volTensorField> tgradUa = fvc::grad(Ua);
+        const volTensorField::Boundary& gradUabf =
+            tgradUa.ref().boundaryField();
+        for (const label patchI : sensitivityPatchIDs_)
+        {
+            const fvPatch& patch = mesh_.boundary()[patchI];
+            tmp<vectorField> tnf = patch.nf();
+            const vectorField& nf = tnf();
+            vectorField gradUaNf
+                (
+                    useSnGradInTranposeStresses_
+                  ? (Ua.boundaryField()[patchI].snGrad() & nf)*nf
+                  : (gradUabf[patchI] & nf)
+                );
+            wallFaceSens_()[patchI] -=
+                nuEff.boundaryField()[patchI]
+               *(gradUaNf & U.boundaryField()[patchI].snGrad())*tnf;
+        }
+    }
+
     // The face-based part of the sensitivities, i.e. terms that multiply
     // dxFace/dxPoint.
     for (const label patchI : sensitivityPatchIDs_)
@@ -517,38 +557,6 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt)
           * nf
         );
 
-        vectorField gradStressTerm(patch.size(), Zero);
-        if (includeGradStressTerm_)
-        {
-            // Terms corresponding to contributions from converting delta to
-            // thetas are added through the corresponding adjoint boundary
-            // conditions instead of grabing contributions from the objective
-            // function. Useful to have a unified formulation for low- and
-            // high-re meshes
-            const fvPatchVectorField& Uab = Ua.boundaryField()[patchI];
-            gradStressTerm = (-((Uab & nf)*gradp.boundaryField()[patchI]));
-            gradStressTerm +=
-            (
-                Uab.component(0)*gradStressX.boundaryField()[patchI]
-              + Uab.component(1)*gradStressY.boundaryField()[patchI]
-              + Uab.component(2)*gradStressZ.boundaryField()[patchI]
-            ) & nf;
-        }
-
-        if (includeTransposeStresses_)
-        {
-            vectorField gradUaNf
-                (
-                    useSnGradInTranposeStresses_ ?
-                    (Ua.boundaryField()[patchI].snGrad() & nf)*nf :
-                    (gradUa.boundaryField()[patchI] & nf)
-                );
-            stressTerm -=
-                nuEff.boundaryField()[patchI]
-               *(gradUaNf & U.boundaryField()[patchI].snGrad())
-               *nf;
-        }
-
         if (includeDivTerm_)
         {
             stressTerm +=
@@ -566,7 +574,7 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt)
         {
             pressureTerm =
             (
-                (nf * pa.boundaryField()[patchI])
+                (nf*pa.boundaryField()[patchI])
               & U.boundaryField()[patchI].snGrad()
             )
            *nf;
@@ -597,7 +605,6 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt)
         wallFaceSens_()[patchI] +=
         (
             stressTerm
-          + gradStressTerm
           + pressureTerm
           + adjointTMsensitivities[patchI]
           + dxdbMultiplierTot

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.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,13 @@ sensitivityVolBSplinesFI::sensitivityVolBSplinesFI
 
 void sensitivityVolBSplinesFI::assembleSensitivities()
 {
+    /*
+    addProfiling
+    (
+        sensitivityVolBSplinesFI,
+        "sensitivityVolBSplinesFI::assembleSensitivities"
+    );
+    */
     read();
 
     // Interpolation engine
@@ -147,8 +154,10 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
             label globalCP = passedCPs + cpI;
 
             // Parameterization info
-            tmp<pointTensorField> dxdbI(boxes[iNURB].getDxDb(cpI));
-            tmp<volTensorField> tvolDxDbI(volPointInter.interpolate(dxdbI));
+            tmp<volTensorField> tvolDxDbI
+            (
+                volPointInter.interpolate(boxes[iNURB].getDxDb(cpI))
+            );
             const volTensorField& volDxDbI = tvolDxDbI();
 
             // Chain rule used to get dx/db at cells
@@ -158,53 +167,47 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
             volTensorField& volDxDbI = tvolDxDbI.ref();
             */
 
-            // Gradient of parameterization info
-            volVectorField temp
-            (
-                IOobject
+            const tensorField& gradDxDbMultInt = gradDxDbMult_.primitiveField();
+            for (label idir = 0; idir < pTraits<vector>::nComponents; ++idir)
+            {
+                // Gradient of parameterization info
+                auto ttemp =
+                    tmp<volVectorField>::New
+                    (
+                        IOobject
+                        (
+                            "dxdb",
+                            mesh_.time().timeName(),
+                            mesh_,
+                            IOobject::NO_READ,
+                            IOobject::NO_WRITE
+                        ),
+                        mesh_,
+                        dimensionedVector(dimless, Zero)
+                    );
+                volVectorField& temp = ttemp.ref();
+                unzipCol(volDxDbI, vector::components(idir), temp);
+
+                volTensorField gradDxDb(fvc::grad(temp));
+                // Volume integral terms
+                flowSens_[globalCP].component(idir) = gSum
                 (
-                    "dxdb",
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh_,
-                vector::zero
-            );
+                    (gradDxDbMultInt && gradDxDb.primitiveField())
+                   *mesh_.V()
+                );
 
-            temp.replace(0, volDxDbI.component(0));
-            temp.replace(1, volDxDbI.component(3));
-            temp.replace(2, volDxDbI.component(6));
-            volTensorField gradDxDb1(fvc::grad(temp));
-
-            temp.replace(0, volDxDbI.component(1));
-            temp.replace(1, volDxDbI.component(4));
-            temp.replace(2, volDxDbI.component(7));
-            volTensorField gradDxDb2(fvc::grad(temp));
-
-            temp.replace(0, volDxDbI.component(2));
-            temp.replace(1, volDxDbI.component(5));
-            temp.replace(2, volDxDbI.component(8));
-            volTensorField gradDxDb3(fvc::grad(temp));
-
-
-            // Volume integral terms
-            flowSens_[globalCP].x() = gSum
-            (
-                (gradDxDbMult_.primitiveField() && gradDxDb1.primitiveField())
-               *mesh_.V()
-            );
-            flowSens_[globalCP].y() = gSum
-            (
-                (gradDxDbMult_.primitiveField() && gradDxDb2.primitiveField())
-               *mesh_.V()
-            );
-            flowSens_[globalCP].z() = gSum
-            (
-                (gradDxDbMult_.primitiveField() && gradDxDb3.primitiveField())
-               *mesh_.V()
-            );
+                if (includeDistance_)
+                {
+                    const tensorField& distSensInt =
+                        distanceSensPtr().primitiveField();
+                    distanceSens_[globalCP].component(idir) =
+                        gSum
+                        (
+                            (distSensInt && gradDxDb.primitiveField())
+                            *mesh_.V()
+                        );
+                }
+            }
 
             // Contribution from objective function term from
             // delta( n dS ) / delta b and
@@ -234,28 +237,6 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
                    *mesh_.V()
                 );
 
-            // Distance dependent term
-            if (includeDistance_)
-            {
-                const tensorField& distSensInt =
-                    distanceSensPtr().primitiveField();
-                distanceSens_[globalCP].x() =
-                    gSum
-                    (
-                        (distSensInt && gradDxDb1.primitiveField())*mesh_.V()
-                    );
-                distanceSens_[globalCP].y() =
-                    gSum
-                    (
-                        (distSensInt && gradDxDb2.primitiveField())*mesh_.V()
-                    );
-                distanceSens_[globalCP].z() =
-                    gSum
-                    (
-                        (distSensInt && gradDxDb3.primitiveField()) *mesh_.V()
-                    );
-            }
-
             // Terms from fvOptions
             optionsSens_[globalCP] +=
                 gSum((optionsDxDbMult_ & volDxDbI.primitiveField())*mesh_.V());
@@ -303,6 +284,8 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
     volBSplinesBase_.boundControlPointMovement(dxdbDirectSens_);
     volBSplinesBase_.boundControlPointMovement(optionsSens_);
     volBSplinesBase_.boundControlPointMovement(bcSens_);
+
+  //profiling::writeNow();
 }
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C

@@ -47,7 +47,7 @@ defineTypeNameAndDebug(shapeSensitivities, 0);
 void shapeSensitivities::accumulateDirectSensitivityIntegrand(const scalar dt)
 {
     // Accumulate direct sensitivities
-    PtrList<objective>& functions(objectiveManager_.getObjectiveFunctions());
+    PtrList<objective>& functions = objectiveManager_.getObjectiveFunctions();
     for (const label patchI : sensitivityPatchIDs_)
     {
         const scalarField magSfDt(mesh_.boundary()[patchI].magSf()*dt);
@@ -65,7 +65,9 @@ void shapeSensitivities::accumulateDirectSensitivityIntegrand(const scalar dt)
 
 void shapeSensitivities::accumulateBCSensitivityIntegrand(const scalar dt)
 {
-    auto& UaBoundary = adjointVars_.Ua().boundaryFieldRef();
+    // Avoid updating the event number to keep consistency with cases caching
+    // gradUa
+    auto& UaBoundary = adjointVars_.Ua().boundaryFieldRef(false);
     tmp<boundaryVectorField> DvDbMult(dvdbMult());
 
     // Accumulate sensitivities due to boundary conditions
@@ -102,7 +104,8 @@ tmp<boundaryVectorField> shapeSensitivities::dvdbMult() const
        turbVars = primalVars_.RASModelVariables();
     const singlePhaseTransportModel& lamTransp = primalVars_.laminarTransport();
     volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
-    volTensorField gradUa(fvc::grad(Ua));
+    tmp<volTensorField> tgradUa = fvc::grad(Ua);
+    const volTensorField::Boundary& gradUabf = tgradUa.ref().boundaryField();
 
     for (const label patchI : sensitivityPatchIDs_)
     {
@@ -115,7 +118,7 @@ tmp<boundaryVectorField> shapeSensitivities::dvdbMult() const
                 nuEff.boundaryField()[patchI]
               * (
                     Ua.boundaryField()[patchI].snGrad()
-                  + (gradUa.boundaryField()[patchI] & nf)
+                  + (gradUabf[patchI] & nf)
                 )
             )
           - (nf*pa.boundaryField()[patchI])

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

@@ -172,42 +172,34 @@ void Foam::incompressibleAdjointSolver::updatePrimalBasedQuantities()
 Foam::tmp<Foam::volTensorField>
 Foam::incompressibleAdjointSolver::computeGradDxDbMultiplier()
 {
-    // Term depending on the adjoint turbulence model
+    /*
+    addProfiling
+    (
+        incompressibleAdjointSolver,
+        "incompressibleAdjointSolver::computeGradDxDbMultiplier"
+    );
+    */
     autoPtr<incompressibleAdjoint::adjointRASModel>& adjointRAS
     (
         getAdjointVars().adjointTurbulence()
     );
-    tmp<volTensorField> tturbulenceTerm(adjointRAS->FISensitivityTerm());
-    volTensorField& turbulenceTerm = tturbulenceTerm.ref();
-
-    // nu effective
-    tmp<volScalarField> tnuEff(adjointRAS->nuEff());
-    const volScalarField& nuEff = tnuEff();
-
-    tmp<volTensorField> tflowTerm
-    (
-        new volTensorField
-        (
-            IOobject
-            (
-               "flowTerm",
-               mesh_.time().timeName(),
-               mesh_,
-               IOobject::NO_READ,
-               IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedTensor(sqr(dimLength)/pow3(dimTime), Zero)
-        )
-    );
-    volTensorField& flowTerm = tflowTerm.ref();
 
     const volScalarField& p = primalVars_.p();
     const volVectorField& U = primalVars_.U();
     const volScalarField& pa = getAdjointVars().pa();
     const volVectorField& Ua = getAdjointVars().Ua();
-    volTensorField gradU(fvc::grad(U));
-    volTensorField gradUa(fvc::grad(Ua));
+
+    // 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& gradU = tgradU.ref();
+    volTensorField::Boundary& gradUbf = gradU.boundaryFieldRef();
 
     // Explicitly correct the boundary gradient to get rid of
     // the tangential component
@@ -217,114 +209,91 @@ Foam::incompressibleAdjointSolver::computeGradDxDbMultiplier()
         if (isA<wallFvPatch>(patch))
         {
             tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
-            const vectorField& nf = tnf();
-            gradU.boundaryFieldRef()[patchI] =
-                nf*U.boundaryField()[patchI].snGrad();
-            //gradUa.boundaryField()[patchI] =
-            //    nf*Ua.boundaryField()[patchI].snGrad();
+            gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad();
         }
     }
 
-    volTensorField stress(nuEff*(gradU + T(gradU)));
-    autoPtr<volVectorField> stressXPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressX", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressYPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressY", stress.dimensions())
-    );
-    autoPtr<volVectorField> stressZPtr
-    (
-        createZeroFieldPtr<vector>(mesh_, "stressZ", stress.dimensions())
-    );
-
-    stressXPtr().replace(0, stress.component(0));
-    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()));
-
-    // Contribution from objective functions and constraints
-    volTensorField objectiveContributions
-    (
-        IOobject
-        (
-            "objectiveContributions",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedTensor(sqr(dimLength)/pow3(dimTime), Zero)
-    );
-    PtrList<objective>& functions
-        (objectiveManagerPtr_->getObjectiveFunctions());
-    forAll(functions, funcI)
-    {
-        objectiveContributions +=
-            functions[funcI].weight()
-           *functions[funcI].gradDxDbMultiplier();
-    }
-
+    tmp<volScalarField> tnuEff = adjointRAS->nuEff();
+    tmp<volSymmTensorField> stress = tnuEff()*twoSymm(gradU);
     // Note:
     // term4 (Ua & grad(stress)) is numerically tricky.  Its div leads to third
-    // order spatial derivs in E-SI based computations Applying the product
+    // order spatial derivs in E-SI based computations. Applying the product
     // derivative rule (putting Ua inside the grad) gives better results in
     // NACA0012, SA, WF.  However, the original formulation should be kept at
     // the boundary in order to respect the Ua boundary conditions (necessary
     // for E-SI to give the same sens as FI).  A mixed approach is hence
     // followed
-    volTensorField term4
-    (
-      - nuEff*(gradUa & (gradU + T(gradU)))
-      + fvc::grad(nuEff * Ua & (gradU + T(gradU)))
-    );
 
-    forAll(mesh_.boundary(), pI)
+    // Term 3, used also to allocated the return field
+    tmp<volTensorField> tgradUa = fvc::grad(Ua);
+    auto tflowTerm =
+        tmp<volTensorField>::New
+        (
+            "flowTerm",
+          - tnuEff*(gradU & twoSymm(tgradUa()))
+        );
+    volTensorField& flowTerm = tflowTerm.ref();
+    // Term 4, only for the internal field
+    flowTerm.ref() +=
+        (
+          - (tgradUa & stress())
+          + fvc::grad(Ua & stress())
+        )().internalField();
+
+    // Boundary conditions from term 4
+    for (label idir = 0; idir < pTraits<vector>::nComponents; ++idir)
     {
-        if (!isA<coupledFvPatch>(mesh_.boundary()[pI]))
+        autoPtr<volVectorField> stressDirPtr
+        (
+            createZeroFieldPtr<vector>
+                (mesh_, "stressDir", stress().dimensions())
+        );
+        // Components need to be in the [0-5] range since stress is a
+        // volSymmTensorField
+        unzipRow(stress(), idir, stressDirPtr());
+        volTensorField gradStressDir(fvc::grad(stressDirPtr()));
+        forAll(mesh_.boundary(), pI)
         {
-            term4.boundaryFieldRef()[pI] =
-                Ua.component(0)().boundaryField()[pI]
-               *gradStressX.boundaryField()[pI]
-              + Ua.component(1)().boundaryField()[pI]
-               *gradStressY.boundaryField()[pI]
-              + Ua.component(2)().boundaryField()[pI]
-               *gradStressZ.boundaryField()[pI];
+            if (!isA<coupledFvPatch>(mesh_.boundary()[pI]))
+            {
+                flowTerm.boundaryFieldRef()[pI] +=
+                    Ua.component(idir)().boundaryField()[pI]
+                   *gradStressDir.boundaryField()[pI];
+            }
         }
     }
+    // Release memory
+    stress.clear();
 
     // Compute dxdb multiplier
-    flowTerm =
+    flowTerm +=
         // Term 1, ATC
         ATCModel_->getFISensitivityTerm()
         // Term 2
-      - fvc::grad(p) * Ua
-        // Term 3
-      - nuEff*(gradU & (gradUa + T(gradUa)))
-        // Term 4
-      + term4
-        // Term 5
-      + (pa * gradU)
-        // Term 6, from the adjoint turbulence model
-      + turbulenceTerm.T()
-        // Term 7, term from objective functions
-      + objectiveContributions;
+      - fvc::grad(p)*Ua;
+
+    // Term 5
+    flowTerm += pa*tgradU;
+
+    // Term 6, from the adjoint turbulence model
+    flowTerm += T(adjointRAS->FISensitivityTerm());
+
+    // Term 7, term from objective functions
+    PtrList<objective>& functions
+        (objectiveManagerPtr_->getObjectiveFunctions());
+
+    for (objective& objI : functions)
+    {
+        if (objI.hasGradDxDbMult())
+        {
+            flowTerm += objI.weight()*objI.gradDxDbMultiplier();
+        }
+    }
 
     flowTerm.correctBoundaryConditions();
 
+  //profiling::writeNow();
+
     return (tflowTerm);
 }