From 2ad6c9f5d2a2ca03d7257c3dfb6268822bef9b2e Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsis Date: Thu, 19 Oct 2023 14:44:31 +0300 Subject: [PATCH] ENH: adjointRASModel will not crash if set to laminar and the Jacobian of the objective function wrt the turbulence variables is called (rare/unorthodox case). Additionally, objectivePowerDissipation dissipation can now be used in topology optimisation, adding the necessary blockage dependency to it. --- .../objectiveIncompressible.C | 2 +- .../objectiveNutSqr/objectiveNutSqr.C | 2 + .../objectivePowerDissipation.C | 60 +++++++++++++++++++ .../objectivePowerDissipation.H | 3 + .../adjointRASModel/adjointRASModel.C | 26 ++++---- 5 files changed, 82 insertions(+), 11 deletions(-) diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C index 8fb87d88e6..68ba25e6c7 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C @@ -373,7 +373,7 @@ void objectiveIncompressible::update_dJdTMvar const labelList& zones ) { - if (dJdTMvarPtr) + if (dJdTMvarPtr.good()) { // nut Jacobians are currently computed in the adjoint turbulence // models, though they would be better placed within the primal diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveNutSqr/objectiveNutSqr.C b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveNutSqr/objectiveNutSqr.C index 96e868c086..972f8bdef0 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveNutSqr/objectiveNutSqr.C +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveNutSqr/objectiveNutSqr.C @@ -51,6 +51,8 @@ addToRunTimeSelectionTable ); +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + void objectiveNutSqr::populateFieldNames() { if (adjointTurbulenceNames_.empty()) diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.C b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.C index 2a910db9a2..aee1298e08 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.C +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.C @@ -29,6 +29,7 @@ License #include "objectivePowerDissipation.H" #include "incompressibleAdjointSolver.H" #include "createZeroField.H" +#include "topOVariablesBase.H" #include "IOmanip.H" #include "addToRunTimeSelectionTable.H" @@ -51,6 +52,8 @@ addToRunTimeSelectionTable ); +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + void objectivePowerDissipation::populateFieldNames() { if (fieldNames_.size() == 1) @@ -134,6 +137,9 @@ objectivePowerDissipation::objectivePowerDissipation sqr(dimLength)/pow3(dimTime) ) ); + + // Allocate direct sensitivity contributions for topology optimisation + dJdbPtr_.reset(createZeroFieldPtr(mesh_, "dJdb", dimless)); } @@ -157,6 +163,19 @@ scalar objectivePowerDissipation::J() scalarField integrandZone(integrand.primitiveField(), zoneI); J_ += 0.5*gSum(integrandZone*VZone); + if (mesh_.foundObject("topoVars")) + { + const topOVariablesBase& vars = + mesh_.lookupObject("topoVars"); + const volScalarField& beta = vars.beta(); + scalar porosityContr = Zero; + for (const label cellI : zoneI) + { + porosityContr += beta[cellI]*magSqr(U[cellI])*V[cellI]; + } + porosityContr *= vars.getBetaMax(); + J_ += returnReduce(porosityContr, sumOp()); + } } return J_; @@ -204,6 +223,23 @@ void objectivePowerDissipation::update_dJdv() dJdvPtr_()[cellI] += integrand[cellI]; } } + + // Add source from porosity dependencies + if (mesh_.foundObject("topoVars")) + { + const topOVariablesBase& vars = + mesh_.lookupObject("topoVars"); + const volScalarField& beta = vars.beta(); + const scalar betaMax = vars.getBetaMax(); + for (const label zI : zones_) + { + const cellZone& zoneI = mesh_.cellZones()[zI]; + for (const label cellI : zoneI) + { + dJdvPtr_()[cellI] += 2*betaMax*beta[cellI]*U[cellI]; + } + } + } } @@ -276,6 +312,30 @@ void objectivePowerDissipation::update_gradDxDbMultiplier() } +void objectivePowerDissipation::update_dJdb() +{ + if (mesh_.foundObject("topoVars")) + { + scalarField& dJdb = dJdbPtr_().primitiveFieldRef(); + dJdb = Zero; + const volVectorField& U = vars_.UInst(); + const topOVariablesBase& vars = + mesh_.lookupObject("topoVars"); + const scalar betaMax = vars.getBetaMax(); + for (const label zI : zones_) + { + const cellZone& zoneI = mesh_.cellZones()[zI]; + for (const label cellI : zoneI) + { + dJdb[cellI] += betaMax*magSqr(U[cellI]); + } + } + } +} + + + + void objectivePowerDissipation::addSource(fvScalarMatrix& matrix) { populateFieldNames(); diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.H b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.H index 7f2bb93a48..70d10f44af 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.H +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectivePowerDissipation/objectivePowerDissipation.H @@ -114,6 +114,9 @@ public: //- Update grad(dx/db multiplier). Volume-based sensitivity term virtual void update_gradDxDbMultiplier(); + //- Contribution to field sensitivities + virtual void update_dJdb(); + //- Add source terms to the adjoint turbulence model equations virtual void addSource(fvScalarMatrix& matrix); }; diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C index 50f7c226da..c98b5d5e18 100644 --- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C @@ -380,6 +380,13 @@ const wordList& adjointRASModel::getAdjointTMVariablesBaseNames() const tmp adjointRASModel::nutJacobianTMVar1() const { + dimensionSet dims + ( + bool(adjointTMVariable1Ptr_) + ? dimViscosity/adjointTMVariable1Ptr_().dimensions() + : dimless + ); + return tmp::New ( @@ -392,17 +399,20 @@ tmp adjointRASModel::nutJacobianTMVar1() const IOobject::NO_WRITE ), mesh_, - dimensionedScalar - ( - dimViscosity/adjointTMVariable1Ptr_().dimensions(), - Zero - ) + dimensionedScalar(dims, Zero) ); } tmp adjointRASModel::nutJacobianTMVar2() const { + dimensionSet dims + ( + bool(adjointTMVariable2Ptr_) + ? dimViscosity/adjointTMVariable2Ptr_().dimensions() + : dimless + ); + return tmp::New ( @@ -415,11 +425,7 @@ tmp adjointRASModel::nutJacobianTMVar2() const IOobject::NO_WRITE ), mesh_, - dimensionedScalar - ( - dimViscosity/adjointTMVariable2Ptr_().dimensions(), - Zero - ) + dimensionedScalar(dims, Zero) ); }