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); }