From 36ca117192d24afc0b127ebd8a161f7f994d59ac Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsis Date: Thu, 9 Dec 2021 18:53:14 +0200 Subject: [PATCH] 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. --- .../FIBase/FIBaseIncompressible.C | 12 +- .../sensitivitySurfaceIncompressible.C | 191 ++++++++--------- .../sensitivitySurfacePointsIncompressible.C | 189 ++++++++--------- .../sensitivityVolBSplinesFIIncompressible.C | 125 +++++------- .../shapeSensitivitiesIncompressible.C | 11 +- .../incompressibleAdjointSolver.C | 193 ++++++++---------- 6 files changed, 346 insertions(+), 375 deletions(-) diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C index 5f787c33fd..c41eda0013 100644 --- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C +++ b/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& 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); - } diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C index 64c49e5f5f..f6eb845d09 100644 --- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C +++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C @@ -602,61 +602,6 @@ void sensitivitySurface::accumulateIntegrand(const scalar dt) autoPtr& adjointTurbulence = adjointVars_.adjointTurbulence(); - Info<< " Calculating auxiliary quantities " << endl; - // Fields needed to calculate adjoint sensitivities - const autoPtr& - 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(patch)) - { - tmp 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 stressXPtr - ( - createZeroFieldPtr(mesh_, "stressX", stress.dimensions()) - ); - autoPtr stressYPtr - ( - createZeroFieldPtr(mesh_, "stressY", stress.dimensions()) - ); - autoPtr stressZPtr - ( - createZeroFieldPtr(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 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 tgradp = fvc::grad(p); + const volVectorField& gradp = tgradp.ref(); + for (const label patchI : sensitivityPatchIDs_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + tmp 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 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(patch)) + { + tmp tnf = mesh_.boundary()[patchI].nf(); + gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad(); + } + } + + tmp tstress = nuEff*twoSymm(tgradU); + const volSymmTensorField& stress = tstress.cref(); + autoPtr ptemp + (Foam::createZeroFieldPtr(mesh_, "temp", sqr(dimVelocity))); + volVectorField& temp = ptemp.ref(); + for (label idir = 0; idir < pTraits::nComponents; ++idir) + { + unzipRow(stress, idir, temp); + volTensorField gradStressDir(fvc::grad(temp)); + for (const label patchI : sensitivityPatchIDs_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + tmp 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 tgradUa = fvc::grad(Ua); + const volTensorField::Boundary& gradUabf = + tgradUa.ref().boundaryField(); + + for (const label patchI : sensitivityPatchIDs_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + tmp 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 diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C index 355c3dca0f..3c9bb45ffc 100644 --- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C +++ b/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& 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(patch)) - { - tmp 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 stressXPtr - ( - createZeroFieldPtr(mesh_, "stressX", stress.dimensions()) - ); - autoPtr stressYPtr - ( - createZeroFieldPtr(mesh_, "stressY", stress.dimensions()) - ); - autoPtr stressZPtr - ( - createZeroFieldPtr(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 meshMovementSensPtr(nullptr); if (includeMeshMovement_) { meshMovementSolver_->accumulateIntegrand(dt); @@ -491,6 +436,101 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt) DebugInfo << " Calculating adjoint sensitivity. " << endl; + tmp 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 tgradp = fvc::grad(p); + const volVectorField& gradp = tgradp.ref(); + for (const label patchI : sensitivityPatchIDs_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + tmp 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 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(patch)) + { + tmp tnf = mesh_.boundary()[patchI].nf(); + gradUbf[patchI] = tnf*U.boundaryField()[patchI].snGrad(); + } + } + + tmp tstress = nuEff*twoSymm(tgradU); + const volSymmTensorField& stress = tstress.cref(); + autoPtr ptemp + (Foam::createZeroFieldPtr(mesh_, "temp", sqr(dimVelocity))); + volVectorField& temp = ptemp.ref(); + for (label idir = 0; idir < pTraits::nComponents; ++idir) + { + unzipRow(stress, idir, temp); + volTensorField gradStressDir(fvc::grad(temp)); + for (const label patchI : sensitivityPatchIDs_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + tmp 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 tgradUa = fvc::grad(Ua); + const volTensorField::Boundary& gradUabf = + tgradUa.ref().boundaryField(); + for (const label patchI : sensitivityPatchIDs_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + tmp 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 diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C index 42b0af9e5c..abd779ab9e 100644 --- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C +++ b/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 dxdbI(boxes[iNURB].getDxDb(cpI)); - tmp tvolDxDbI(volPointInter.interpolate(dxdbI)); + tmp 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::nComponents; ++idir) + { + // Gradient of parameterization info + auto ttemp = + tmp::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(); } diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C index 9b8ddc8b7a..be4ae28bb5 100644 --- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C +++ b/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& functions(objectiveManager_.getObjectiveFunctions()); + PtrList& 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 DvDbMult(dvdbMult()); // Accumulate sensitivities due to boundary conditions @@ -102,7 +104,8 @@ tmp shapeSensitivities::dvdbMult() const turbVars = primalVars_.RASModelVariables(); const singlePhaseTransportModel& lamTransp = primalVars_.laminarTransport(); volScalarField nuEff(lamTransp.nu() + turbVars->nutRef()); - volTensorField gradUa(fvc::grad(Ua)); + tmp tgradUa = fvc::grad(Ua); + const volTensorField::Boundary& gradUabf = tgradUa.ref().boundaryField(); for (const label patchI : sensitivityPatchIDs_) { @@ -115,7 +118,7 @@ tmp shapeSensitivities::dvdbMult() const nuEff.boundaryField()[patchI] * ( Ua.boundaryField()[patchI].snGrad() - + (gradUa.boundaryField()[patchI] & nf) + + (gradUabf[patchI] & nf) ) ) - (nf*pa.boundaryField()[patchI]) 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 2ffa2d1b03..c5691c4477 100644 --- a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.C +++ b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/incompressibleAdjointSolver/incompressibleAdjointSolver.C @@ -172,42 +172,34 @@ void Foam::incompressibleAdjointSolver::updatePrimalBasedQuantities() Foam::tmp Foam::incompressibleAdjointSolver::computeGradDxDbMultiplier() { - // Term depending on the adjoint turbulence model + /* + addProfiling + ( + incompressibleAdjointSolver, + "incompressibleAdjointSolver::computeGradDxDbMultiplier" + ); + */ autoPtr& adjointRAS ( getAdjointVars().adjointTurbulence() ); - tmp tturbulenceTerm(adjointRAS->FISensitivityTerm()); - volTensorField& turbulenceTerm = tturbulenceTerm.ref(); - - // nu effective - tmp tnuEff(adjointRAS->nuEff()); - const volScalarField& nuEff = tnuEff(); - - tmp 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 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(patch)) { tmp 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 stressXPtr - ( - createZeroFieldPtr(mesh_, "stressX", stress.dimensions()) - ); - autoPtr stressYPtr - ( - createZeroFieldPtr(mesh_, "stressY", stress.dimensions()) - ); - autoPtr stressZPtr - ( - createZeroFieldPtr(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& functions - (objectiveManagerPtr_->getObjectiveFunctions()); - forAll(functions, funcI) - { - objectiveContributions += - functions[funcI].weight() - *functions[funcI].gradDxDbMultiplier(); - } - + tmp tnuEff = adjointRAS->nuEff(); + tmp 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 tgradUa = fvc::grad(Ua); + auto tflowTerm = + tmp::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::nComponents; ++idir) { - if (!isA(mesh_.boundary()[pI])) + autoPtr stressDirPtr + ( + createZeroFieldPtr + (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(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& functions + (objectiveManagerPtr_->getObjectiveFunctions()); + + for (objective& objI : functions) + { + if (objI.hasGradDxDbMult()) + { + flowTerm += objI.weight()*objI.gradDxDbMultiplier(); + } + } flowTerm.correctBoundaryConditions(); + //profiling::writeNow(); + return (tflowTerm); }