From b550a23acbfd680ea86b109d1ef3aab67853f0e6 Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsis Date: Mon, 7 Feb 2022 12:01:34 +0200 Subject: [PATCH] ENH: add the infrastructure for computing and utilising the Jacobian of an objective function, defined at the boundary, wrt nut and gradU. Also modified the current objectives that include such contributions --- .../boundaryAdjointContribution.C | 12 ++++ .../boundaryAdjointContribution.H | 2 + ...oundaryAdjointContributionIncompressible.C | 24 +++++++ ...oundaryAdjointContributionIncompressible.H | 2 + .../objectiveForce/objectiveForce.C | 28 +++++++- .../objectiveForce/objectiveForce.H | 7 +- .../objectiveIncompressible.C | 69 ++++++++++++++++++- .../objectiveIncompressible.H | 26 +++++++ .../objectiveIncompressibleI.H | 12 ++++ .../objectiveMoment/objectiveMoment.C | 21 ++++++ .../objectiveMoment/objectiveMoment.H | 8 +++ .../adjoint/objectives/objective/objective.C | 29 -------- .../adjoint/objectives/objective/objective.H | 14 ---- .../adjoint/objectives/objective/objectiveI.H | 6 -- 14 files changed, 204 insertions(+), 56 deletions(-) diff --git a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.C b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.C index d1696041de..004ea4b0ae 100644 --- a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.C +++ b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.C @@ -103,6 +103,18 @@ tmp boundaryAdjointContribution::adjointTMVariable2Source() } +tmp boundaryAdjointContribution::dJdnut() +{ + return tmp::New(patch_.size(), Zero); +} + + +tmp boundaryAdjointContribution::dJdGradU() +{ + return tmp::New(patch_.size(), Zero); +} + + tmp boundaryAdjointContribution::TMVariable1Diffusion() { return tmp::New(patch_.size(), Zero); diff --git a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.H b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.H index d8063d0407..b9e52b8294 100644 --- a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.H +++ b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContribution/boundaryAdjointContribution.H @@ -138,6 +138,8 @@ public: virtual tmp normalVelocitySource() = 0; virtual tmp adjointTMVariable1Source(); virtual tmp adjointTMVariable2Source(); + virtual tmp dJdnut(); + virtual tmp dJdGradU(); virtual tmp energySource() = 0; virtual tmp momentumDiffusion() = 0; diff --git a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C index 6f108478a0..3b01bb971e 100644 --- a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C +++ b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.C @@ -223,6 +223,30 @@ boundaryAdjointContributionIncompressible::adjointTMVariable2Source() } +tmp +boundaryAdjointContributionIncompressible::dJdnut() +{ + return + sumContributions + ( + objectiveManager_.getObjectiveFunctions(), + &objectiveIncompressible::boundarydJdnut + ); +} + + +tmp +boundaryAdjointContributionIncompressible::dJdGradU() +{ + return + sumContributions + ( + objectiveManager_.getObjectiveFunctions(), + &objectiveIncompressible::boundarydJdGradU + ); +} + + tmp boundaryAdjointContributionIncompressible::momentumDiffusion() { tmp tnuEff(new scalarField(patch_.size(), Zero)); diff --git a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.H index 356fcac8ac..d77cec9b25 100644 --- a/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.H +++ b/src/optimisation/adjointOptimisation/adjoint/boundaryAdjointContributions/boundaryAdjointContributionIncompressible/boundaryAdjointContributionIncompressible.H @@ -139,6 +139,8 @@ public: tmp energySource(); tmp adjointTMVariable1Source(); tmp adjointTMVariable2Source(); + tmp dJdnut(); + tmp dJdGradU(); tmp momentumDiffusion(); tmp laminarDiffusivity(); diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.C b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.C index 2ab4bfe870..18e694f1d8 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.C +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.C @@ -113,7 +113,8 @@ objectiveForce::objectiveForce bdJdpPtr_.reset(createZeroBoundaryPtr(mesh_)); bdSdbMultPtr_.reset(createZeroBoundaryPtr(mesh_)); bdxdbMultPtr_.reset(createZeroBoundaryPtr(mesh_)); - bdJdStressPtr_.reset(createZeroBoundaryPtr(mesh_)); + bdJdnutPtr_.reset(createZeroBoundaryPtr(mesh_)); + bdJdGradUPtr_.reset(createZeroBoundaryPtr(mesh_)); } @@ -268,14 +269,35 @@ void objectiveForce::update_dxdbMultiplier() } -void objectiveForce::update_dJdStressMultiplier() +void objectiveForce::update_boundarydJdnut() { + const volVectorField& U = vars_.U(); + volSymmTensorField devGradU(dev(twoSymm(fvc::grad(U)))); + for (const label patchI : forcePatches_) { const fvPatch& patch = mesh_.boundary()[patchI]; tmp tnf = patch.nf(); const vectorField& nf = tnf(); - bdJdStressPtr_()[patchI] = (forceDirection_ * nf)/denom(); + bdJdnutPtr_()[patchI] = + -((devGradU.boundaryField()[patchI] & forceDirection_) & nf)/denom(); + } +} + + +void objectiveForce::update_boundarydJdGradU() +{ + const autoPtr& + turbVars = vars_.RASModelVariables(); + const singlePhaseTransportModel& lamTransp = vars_.laminarTransport(); + volScalarField nuEff(lamTransp.nu() + turbVars->nutRef()); + for (const label patchI : forcePatches_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + const vectorField& Sf = patch.Sf(); + bdJdGradUPtr_()[patchI] = + - nuEff.boundaryField()[patchI] + *dev(forceDirection_*Sf + Sf*forceDirection_); } } diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.H b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.H index 4cff0bf3fe..01bf362ae2 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.H +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveForce/objectiveForce.H @@ -111,8 +111,11 @@ public: //- Update delta(x)/delta b multiplier void update_dxdbMultiplier(); - //- Update dJ/dStress multiplier - void update_dJdStressMultiplier(); + //- Update dJ/dnut multiplier + void update_boundarydJdnut(); + + //- Update dJ/dGradU multiplier + void update_boundarydJdGradU(); //- Return denominator, without density virtual scalar denom() const; diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C index af2ae64ee9..6c71f1e66a 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C @@ -84,7 +84,9 @@ objectiveIncompressible::objectiveIncompressible bdJdpPtr_(nullptr), bdJdTPtr_(nullptr), bdJdTMvar1Ptr_(nullptr), - bdJdTMvar2Ptr_(nullptr) + bdJdTMvar2Ptr_(nullptr), + bdJdnutPtr_(nullptr), + bdJdGradUPtr_(nullptr) { weight_ = dict.get("weight"); computeMeanFields_ = vars_.computeMeanFields(); @@ -182,6 +184,14 @@ void objectiveIncompressible::doNormalization() { bdJdTMvar2Ptr_() *= oneOverNorm; } + if (hasBoundarydJdnut()) + { + bdJdnutPtr_() *= oneOverNorm; + } + if (hasBoundarydJdGradU()) + { + bdJdGradUPtr_() *= oneOverNorm; + } // Normalize geometric fields objective::doNormalization(); @@ -375,6 +385,32 @@ const fvPatchScalarField& objectiveIncompressible::boundarydJdTMvar2 } +const fvPatchScalarField& objectiveIncompressible::boundarydJdnut +( + const label patchI +) +{ + if (!bdJdnutPtr_) + { + bdJdnutPtr_.reset(createZeroBoundaryPtr(mesh_)); + } + return bdJdnutPtr_()[patchI]; +} + + +const fvPatchTensorField& objectiveIncompressible::boundarydJdGradU +( + const label patchI +) +{ + if (!bdJdGradUPtr_) + { + bdJdGradUPtr_.reset(createZeroBoundaryPtr(mesh_)); + } + return bdJdGradUPtr_()[patchI]; +} + + const boundaryVectorField& objectiveIncompressible::boundarydJdv() { if (!bdJdvPtr_) @@ -445,6 +481,26 @@ const boundaryScalarField& objectiveIncompressible::boundarydJdTMvar2() } +const boundaryScalarField& objectiveIncompressible::boundarydJdnut() +{ + if (!bdJdnutPtr_) + { + bdJdnutPtr_.reset(createZeroBoundaryPtr(mesh_)); + } + return bdJdnutPtr_(); +} + + +const boundaryTensorField& objectiveIncompressible::boundarydJdGradU() +{ + if (!bdJdGradUPtr_) + { + bdJdGradUPtr_.reset(createZeroBoundaryPtr(mesh_)); + } + return *bdJdGradUPtr_; +} + + void objectiveIncompressible::update() { // Objective function value @@ -472,13 +528,14 @@ void objectiveIncompressible::update() update_boundarydJdT(); update_boundarydJdTMvar1(); update_boundarydJdTMvar2(); + update_boundarydJdnut(); + update_boundarydJdGradU(); update_boundarydJdb(); update_dSdbMultiplier(); update_dndbMultiplier(); update_dxdbMultiplier(); update_dxdbDirectMultiplier(); update_boundaryEdgeContribution(); - update_dJdStressMultiplier(); // Divide everything with normalization factor doNormalization(); @@ -539,6 +596,14 @@ void objectiveIncompressible::nullify() { bdJdTMvar2Ptr_() == scalar(0); } + if (hasBoundarydJdnut()) + { + bdJdnutPtr_() == scalar(0); + } + if (hasBoundarydJdGradU()) + { + bdJdGradUPtr_() == tensor::zero; + } // Nullify geometric fields and sets nullified_ to true objective::nullify(); diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H index a96a07b003..b7a01154c1 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.H @@ -97,6 +97,12 @@ protected: //- Adjoint outlet turbulence model var 2 autoPtr bdJdTMvar2Ptr_; + //- Jacobian wrt to nut + autoPtr bdJdnutPtr_; + + //- Term multiplying gradU variations + autoPtr bdJdGradUPtr_; + private: @@ -207,6 +213,12 @@ public: //- patch const fvPatchScalarField& boundarydJdTMvar2(const label); + //- Objective partial deriv wrt nut for a specific patch + const fvPatchScalarField& boundarydJdnut(const label); + + //- Objective partial deriv wrt stress tensor + const fvPatchTensorField& boundarydJdGradU(const label); + //- Objective partial deriv wrt velocity for all patches const boundaryVectorField& boundarydJdv(); @@ -228,6 +240,12 @@ public: //- Objective partial deriv wrt turbulence model var 2 for all patches const boundaryScalarField& boundarydJdTMvar2(); + //- Objective partial deriv wrt nut for all patches + const boundaryScalarField& boundarydJdnut(); + + //- Objective partial deriv wrt gradU + const boundaryTensorField& boundarydJdGradU(); + //- Update objective function derivatives virtual void update(); @@ -282,6 +300,12 @@ public: virtual void update_boundarydJdTMvar2() {} + virtual void update_boundarydJdnut() + {} + + virtual void update_boundarydJdGradU() + {} + virtual void update_boundarydJdb() {} @@ -321,6 +345,8 @@ public: inline bool hasBoundarydJdT() const; inline bool hasBoundarydJdTMVar1() const; inline bool hasBoundarydJdTMVar2() const; + inline bool hasBoundarydJdnut() const; + inline bool hasBoundarydJdGradU() const; }; diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressibleI.H b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressibleI.H index 4b1bb5f442..a1f72dafb3 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressibleI.H +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressibleI.H @@ -102,4 +102,16 @@ inline bool Foam::objectiveIncompressible::hasBoundarydJdTMVar2() const } +inline bool Foam::objectiveIncompressible::hasBoundarydJdnut() const +{ + return bool(bdJdnutPtr_); +} + + +inline bool Foam::objectiveIncompressible::hasBoundarydJdGradU() const +{ + return bool(bdJdGradUPtr_); +} + + // ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.C b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.C index d71fcfc21a..b6c3378193 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.C +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.C @@ -120,6 +120,8 @@ objectiveMoment::objectiveMoment bdSdbMultPtr_.reset(createZeroBoundaryPtr(mesh_)); bdxdbMultPtr_.reset(createZeroBoundaryPtr(mesh_)); bdxdbDirectMultPtr_.reset(createZeroBoundaryPtr(mesh_)); + bdJdnutPtr_.reset(createZeroBoundaryPtr(mesh_)); + //bdJdGradUPtr_.reset(createZeroBoundaryPtr(mesh_)); } @@ -306,6 +308,25 @@ void objectiveMoment::update_dxdbDirectMultiplier() } +void objectiveMoment::update_boundarydJdnut() +{ + const volVectorField& U = vars_.U(); + volSymmTensorField devGradU(dev(twoSymm(fvc::grad(U)))); + + for (const label patchI : momentPatches_) + { + const fvPatch& patch = mesh_.boundary()[patchI]; + tmp nf(patch.nf()); + const vectorField dx(patch.Cf() - rotationCentre_); + bdJdnutPtr_()[patchI] = + -rhoInf_ + *( + (dx^(devGradU.boundaryField()[patchI] & nf)) & momentDirection_ + )*invDenom_; + } +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace objectives diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.H b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.H index 9a369e9d10..c5a54afb9b 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.H +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveMoment/objectiveMoment.H @@ -118,6 +118,14 @@ public: //- Update delta(x)/delta b multiplier coming directly from the //- objective void update_dxdbDirectMultiplier(); + + //- Update dJ/dnut multiplier + void update_boundarydJdnut(); + + //- Update dJ/dGradU multiplier + /* WIP + void update_boundarydJdGradU(); + */ }; diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.C b/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.C index 0eb776d1d3..86bb2b8fb2 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.C +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.C @@ -156,7 +156,6 @@ objective::objective bdxdbMultPtr_(nullptr), bdxdbDirectMultPtr_(nullptr), bEdgeContribution_(nullptr), - bdJdStressPtr_(nullptr), divDxDbMultPtr_(nullptr), gradDxDbMultPtr_(nullptr), @@ -376,10 +375,6 @@ void objective::doNormalization() { gradDxDbMultPtr_() *= oneOverNorm; } - if (hasBoundarydJdStress()) - { - bdJdStressPtr_() *= oneOverNorm; - } } } @@ -507,16 +502,6 @@ const vectorField& objective::boundaryEdgeMultiplier } -const fvPatchTensorField& objective::boundarydJdStress(const label patchI) -{ - if (!bdJdStressPtr_) - { - bdJdStressPtr_.reset(createZeroBoundaryPtr(mesh_)); - } - return bdJdStressPtr_()[patchI]; -} - - const boundaryVectorField& objective::boundarydJdb() { if (!bdJdbPtr_) @@ -580,16 +565,6 @@ const vectorField3& objective::boundaryEdgeMultiplier() } -const boundaryTensorField& objective::boundarydJdStress() -{ - if (!bdJdStressPtr_) - { - bdJdStressPtr_.reset(createZeroBoundaryPtr(mesh_)); - } - return *bdJdStressPtr_; -} - - const volScalarField& objective::divDxDbMultiplier() { if (!divDxDbMultPtr_) @@ -666,10 +641,6 @@ void objective::nullify() field = vector::zero; } } - if (hasBoundarydJdStress()) - { - bdJdStressPtr_() == tensor::zero; - } if (hasDivDxDbMult()) { divDxDbMultPtr_() == diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.H b/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.H index f95242076b..b290784f0e 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.H +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.H @@ -127,9 +127,6 @@ protected: //- NURBS surfaces G1 continuity autoPtr bEdgeContribution_; - //- For use with discrete-like sensitivities - autoPtr bdJdStressPtr_; - // Contribution to volume-based sensitivities from volume-based // objective functions //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -298,9 +295,6 @@ public: const label edgeI ); - //- Objective partial deriv wrt stress tensor - const fvPatchTensorField& boundarydJdStress(const label); - //- Contribution to surface sensitivities for all patches const boundaryVectorField& boundarydJdb(); @@ -319,9 +313,6 @@ public: //- Multiplier located at patch boundary edges const vectorField3& boundaryEdgeMultiplier(); - //- Objective partial deriv wrt stress tensor - const boundaryTensorField& boundarydJdStress(); - //- Multiplier of grad( delta(x)/delta b) for volume-based sensitivities const volScalarField& divDxDbMultiplier(); @@ -363,10 +354,6 @@ public: virtual void update_boundaryEdgeContribution() {} - //- Update dJ/dStress field - virtual void update_dJdStressMultiplier() - {} - //- Update div( dx/db multiplier). Volume-based sensitivity term virtual void update_divDxDbMultiplier() {} @@ -415,7 +402,6 @@ public: inline bool hasdxdbMult() const; inline bool hasdxdbDirectMult() const; inline bool hasBoundaryEdgeContribution() const; - inline bool hasBoundarydJdStress() const; inline bool hasDivDxDbMult() const; inline bool hasGradDxDbMult() const; diff --git a/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objectiveI.H b/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objectiveI.H index 9609184627..4c11fa047d 100644 --- a/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objectiveI.H +++ b/src/optimisation/adjointOptimisation/adjoint/objectives/objective/objectiveI.H @@ -90,12 +90,6 @@ inline bool Foam::objective::hasGradDxDbMult() const } -inline bool Foam::objective::hasBoundarydJdStress() const -{ - return bool(bdJdStressPtr_); -} - - inline bool Foam::objective::hasIntegrationStartTime() const { return bool(integrationStartTimePtr_);