From 8f7fc54c27d5672b3cf2a64de096da04b2a65a72 Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsis Date: Mon, 7 Feb 2022 17:31:35 +0200 Subject: [PATCH] ENH: added the adjoint to the kOmega SST turbulence model MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit for incompressible flows. The typical convention of appending the primal field name with 'a' to form the adjoint field is followed for the adjoint turbulent kinetic energy (i.e. 'ka') but since this would produce an ugly variable name for the adjoint to omega (i.e. omegaa), the latter is abbreviated to 'wa'. The work is based on \verbatim Kavvadias, I., Papoutsis-Kiachagias, E., Dimitrakopoulos, G., & Giannakoglou, K. (2014). The continuous adjoint approach to the k–$omega$ SST turbulence model with applications in shape optimization Engineering Optimization, 47(11), 1523-1542. https://doi.org/10.1080/0305215X.2014.979816 \endverbatim with changes in the discretisation of a number of differential operators and the formulation of the adjoint to the wall functions employed by the primal model. Regarding the latter, the code assumes (and differentiates) the default behaviour of nutkWallFunction (i.e. nutWallFunction::blendingType::STEPWISE) and omegaWallFunction (i.e. omegaWallFunction::blendingType::BINOMIAL2). Due to the availability of a number of terms required for the formulation of the wall function for ka, the latter is implemented within adjointkOmegaSST itself, with contributions from objective functions implemented within kaqRWallFunction. Wall functions for wa are implemented within waWallFunction. The initial implementation of the above-mentioned reference was performed by Dr. Ioannis Kavvadias --- .../adjointOptimisation/adjoint/Make/files | 9 +- .../adjointRASModel/adjointRASModel.C | 7 + .../adjointRASModel/adjointRASModel.H | 7 + .../adjointSpalartAllmaras.C | 2 + .../adjointkOmegaSST/adjointkOmegaSST.C | 2276 +++++++++++++++++ .../adjointkOmegaSST/adjointkOmegaSST.H | 620 +++++ .../adjointkOmegaSSTTemplates.C | 69 + .../adjointFarFieldTMVar1FvPatchScalarField.C | 178 ++ .../adjointFarFieldTMVar1FvPatchScalarField.H | 154 ++ .../adjointFarFieldTMVar2FvPatchScalarField.C | 178 ++ .../adjointFarFieldTMVar2FvPatchScalarField.H | 154 ++ .../adjointOutletFluxFvPatchField.C | 150 ++ .../adjointOutletFluxFvPatchField.H | 169 ++ .../adjointOutletFluxFvPatchFields.C | 46 + .../adjointOutletFluxFvPatchFields.H | 53 + .../adjointOutletFluxFvPatchFieldsFwd.H | 53 + .../adjointOutletKaFvPatchScalarField.C | 142 + .../adjointOutletKaFvPatchScalarField.H | 138 + .../adjointOutletWaFvPatchScalarField.C | 142 + .../adjointOutletWaFvPatchScalarField.H | 138 + .../kaqRWallFunctionFvPatchScalarField.C | 176 ++ .../kaqRWallFunctionFvPatchScalarField.H | 144 ++ .../waWallFunctionFvPatchScalarField.C | 318 +++ .../waWallFunctionFvPatchScalarField.H | 184 ++ .../RAS/RASModelVariables/RASModelVariables.H | 13 +- .../RAS/kOmegaSST/kOmegaSST.C | 92 + .../RAS/kOmegaSST/kOmegaSST.H | 21 + 27 files changed, 5629 insertions(+), 4 deletions(-) create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFieldsFwd.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.H create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.C create mode 100644 src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.H diff --git a/src/optimisation/adjointOptimisation/adjoint/Make/files b/src/optimisation/adjointOptimisation/adjoint/Make/files index b18b9abec6..7dbddc342d 100644 --- a/src/optimisation/adjointOptimisation/adjoint/Make/files +++ b/src/optimisation/adjointOptimisation/adjoint/Make/files @@ -64,11 +64,18 @@ turbulenceModels/incompressibleAdjoint/adjointTurbulenceModel/adjointTurbulenceM turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C turbulenceModels/incompressibleAdjoint/adjointRAS/adjointLaminar/adjointLaminar.C turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C +turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointInletNuaTilda/adjointInletNuaTildaFvPatchScalarField.C turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletNuaTilda/adjointOutletNuaTildaFvPatchScalarField.C turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldNuaTilda/adjointFarFieldNuaTildaFvPatchScalarField.C +turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.C +turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.C turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletNuaTildaFlux/adjointOutletNuaTildaFluxFvPatchScalarField.C -turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchScalarField.C +turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.C +turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.C +turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.C +turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.C +turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.C /* ADJOINT BOUNDARY CONDITIONS */ adjointBoundaryConditions/adjointBoundaryCondition/adjointBoundaryConditions.C 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 15ea65c8b1..8dea7abbaf 100644 --- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C @@ -147,6 +147,7 @@ adjointRASModel::adjointRASModel adjointTMVariable1Ptr_(nullptr), adjointTMVariable2Ptr_(nullptr), + adjointTMVariablesBaseNames_(0), adjointTMVariable1MeanPtr_(nullptr), adjointTMVariable2MeanPtr_(nullptr), adjMomentumBCSourcePtr_( createZeroBoundaryPtr(mesh_) ), @@ -351,6 +352,12 @@ autoPtr& adjointRASModel::getAdjointTMVariable2InstPtr() } +const wordList& adjointRASModel::getAdjointTMVariablesBaseNames() +{ + return adjointTMVariablesBaseNames_; +} + + tmp adjointRASModel::nutJacobianTMVar1() const { return diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.H index 6abd793144..6aa126a34a 100644 --- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.H +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.H @@ -116,6 +116,9 @@ protected: //- Adjoint turbulence model variable 2 autoPtr adjointTMVariable2Ptr_; + //- Base names of the adjoint fields + wordList adjointTMVariablesBaseNames_; + //- Adjoint turbulence model variable 1, mean value autoPtr adjointTMVariable1MeanPtr_; @@ -258,6 +261,10 @@ public: //- Return non-constant autoPtr to adjoint turbulence model variable 2 autoPtr& getAdjointTMVariable2InstPtr(); + //- Return reference to the adjoint turbulence model variables base + //- names + const wordList& getAdjointTMVariablesBaseNames(); + //- Return the effective stress tensor including the laminar stress virtual tmp devReff() const = 0; diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C index 7878349822..c7fffd0a29 100644 --- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointSpalartAllmaras/adjointSpalartAllmaras.C @@ -675,6 +675,8 @@ adjointSpalartAllmaras::adjointSpalartAllmaras gradNuTilda_(fvc::grad(nuTilda())), minStilda_("SMALL", Stilda_.dimensions(), SMALL) { + adjointTMVariablesBaseNames_.setSize(1); + adjointTMVariablesBaseNames_[0] = "nuaTilda"; // Read nuaTilda field and reset pointer to the first // adjoint turbulence model variable variablesSet::setField diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C new file mode 100644 index 0000000000..06bf3a9ea8 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.C @@ -0,0 +1,2276 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "adjointkOmegaSST.H" +#include "wallFvPatch.H" +#include "fixedValueFvPatchFields.H" +#include "zeroGradientFvPatchFields.H" +#include "linear.H" +#include "reverseLinear.H" +#include "nutkWallFunctionFvPatchScalarField.H" +#include "omegaWallFunctionFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressibleAdjoint +{ +namespace adjointRASModels +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(adjointkOmegaSST, 0); +addToRunTimeSelectionTable(adjointRASModel, adjointkOmegaSST, dictionary); + + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +tmp adjointkOmegaSST::F1() const +{ + tmp arg1 = min + ( + min + ( + max + ( + (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + (4*alphaOmega2_)*k()/(CDkOmegaPlus_*sqr(y_)) + ), + scalar(10) + ); + + return tanh(pow4(arg1)); +} + + +tmp adjointkOmegaSST::F2() const +{ + tmp arg2 = min + ( + max + ( + (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + scalar(100) + ); + + return tanh(sqr(arg2)); +} + + +tmp adjointkOmegaSST::GbyNu +( + const volScalarField& GbyNu0, + const volScalarField& F2, + const volScalarField& S2 +) const +{ + return min + ( + GbyNu0, + (c1_/a1_)*betaStar_*omega() + *max(a1_*omega(), b1_*F2*sqrt(S2)) + ); +} + + +tmp adjointkOmegaSST::GbyNu +( + const volScalarField::Internal& GbyNu0, + const volScalarField::Internal& F2, + const volScalarField::Internal& S2 +) const +{ + return min + ( + GbyNu0, + (c1_/a1_)*betaStar_*omega()() + *max(a1_*omega()(), b1_*F2*sqrt(S2)) + ); +} + + +tmp adjointkOmegaSST::zeroFirstCell() +{ + auto tzeroFirstCell = + tmp::New + ( + IOobject + ( + "zeroFirstCell", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Foam::one{}) + ); + volScalarField& zeroFirstCell = tzeroFirstCell.ref(); + + firstCellIDs_.resize(mesh_.nCells(), -1); + label counter(0); + for (const fvPatchScalarField& omegab : omega().boundaryField()) + { + const fvPatch& patch = omegab.patch(); + if (isA(omegab)) + { + const label patchi = patch.index(); + const labelList& faceCells = patch.faceCells(); + fvPatchScalarField& bf = zeroFirstCell.boundaryFieldRef()[patchi]; + forAll(faceCells, faceI) + { + const label cellI = faceCells[faceI]; + + zeroFirstCell[cellI] = 0.; + bf[faceI] = 0.; + firstCellIDs_[counter++] = cellI; + } + } + } + firstCellIDs_.setSize(counter); + + zeroFirstCell.correctBoundaryConditions(); + + return tzeroFirstCell; +} + + +tmp adjointkOmegaSST::dR_dnut() +{ + const volVectorField& U = primalVars_.U(); + const volVectorField& Ua = adjointVars_.UaInst(); + word scheme("coeffsDiff"); + tmp tdR_dnut + ( + dNutdbMult(U, Ua, nutRef(), scheme) + + dNutdbMult(k(), ka(), alphaK_, nutRef(), scheme) + + dNutdbMult(omega(), zeroFirstCell_*wa(), alphaOmega_, nutRef(), scheme) + - case_1_Pk_*ka()*GbyNu0_*zeroFirstCell_ + ); + volScalarField& dRdnut = tdR_dnut.ref(); + + forAll(mesh_.boundary(), pI) + { + const fvPatch& patch = mesh_.boundary()[pI]; + const fvPatchScalarField& nutb = nutRef().boundaryField()[pI]; + const labelList& faceCells = patch.faceCells(); + if (isA(nutb)) + { + fvPatchScalarField& bf = dRdnut.boundaryFieldRef()[pI]; + const scalarField kSnGrad(k().boundaryField()[pI].snGrad()); + const scalarField omegaSnGrad + ( + omega().boundaryField()[pI].snGrad() + ); + const fvPatchScalarField& akb = alphaK_.boundaryField()[pI]; + const fvPatchScalarField& aOmegab = alphaOmega_.boundaryField()[pI]; + const vectorField USnGrad(U.boundaryField()[pI].snGrad()); + const fvPatchTensorField& gradUb = gradU_.boundaryField()[pI]; + const vectorField nf(mesh_.boundary()[pI].nf()); + forAll(faceCells, fI) + { + const label cI(faceCells[fI]); + bf[fI] = + - wa()[cI]*zeroFirstCell_[cI]*aOmegab[fI]*omegaSnGrad[fI] + - ka()[cI]*akb[fI]*kSnGrad[fI] + - (Ua[cI] & (USnGrad[fI] + (dev2(gradUb[fI]) & nf[fI]))); + } + } + } + + return tdR_dnut; +} + + +tmp adjointkOmegaSST::dnut_domega() const +{ + return + ( + - case_1_nut_*k()/sqr(omega()) + - a1_*k()/(b1_*S_*F2_*F2_)*dF2_domega() + ); +} + + +tmp adjointkOmegaSST::dnut_dk() const +{ + return + ( + a1_/max(a1_*omega(), b1_*F2_*S_) + - a1_*k()/(b1_*S_*F2_*F2_)*dF2_dk() + ); +} + + +tmp adjointkOmegaSST::dF2_domega() const +{ + tmp arg2 = min + ( + max + ( + (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + scalar(100) + ); + + return + - scalar(2)*arg2*(1 - F2_*F2_)* + ( + case_2_nut_*scalar(2)*sqrt(k())/(betaStar_*sqr(omega())*y_) + + case_3_nut_*scalar(500)*nu()/sqr(omega()*y_) + ); +} + + +tmp adjointkOmegaSST::dF2_dk() const +{ + tmp arg2 = min + ( + max + ( + (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + scalar(100) + ); + + return + case_2_nut_*scalar(2)*arg2*(1 - F2_*F2_)/(betaStar_*omega()*y_ + *sqrt(k())); +} + + +tmp adjointkOmegaSST::dGPrime_domega() const +{ + return + ( + case_2_GPrime_*c1_*betaStar_/a1_ + *( + max(a1_*omega(), b1_*F2_*S_) + + case_1_nut_*a1_*omega() + + (scalar(1) - case_1_nut_)*omega()*b1_*S_*dF2_domega() + ) + ); +} + + +tmp adjointkOmegaSST::dGPrime_dk() const +{ + return case_2_GPrime_*c1_*betaStar_/a1_*omega()*b1_*S_*dF2_dk(); +} + + +tmp adjointkOmegaSST::dR_dF1() const +{ + const volVectorField& U = primalVars_.U(); + const surfaceScalarField& phi = primalVars_.phi(); + tmp tGPrime = GbyNu(GbyNu0_, F2_, S2_); + tmp tdivU = fvc::div(fvc::absolute(phi, U)); + word scheme("coeffsDiff"); + + tmp tdRdF1 + ( + nutRef() + *( + coeffsDifferentiation(k(), ka(), scheme)*(alphaK1_ - alphaK2_) + + coeffsDifferentiation(omega(), zeroFirstCell_*wa(), scheme) + *(alphaOmega1_ - alphaOmega2_) + ) + ); + volScalarField& dRdF1 = tdRdF1.ref(); + + typedef fixedValueFvPatchScalarField fixedValue; + typedef zeroGradientFvPatchScalarField zeroGrad; + typedef omegaWallFunctionFvPatchScalarField omegaWF; + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& kb = k().boundaryField()[pI]; + const fvPatchScalarField& omegab = omega().boundaryField()[pI]; + fvPatchScalarField& dRdF1b = dRdF1.boundaryFieldRef()[pI]; + if + ( + isA(kb) + && (isA(omegab) || isA(omegab)) + ) + { + dRdF1b = dRdF1b.patchInternalField(); + } + else if (isA(kb) && isA(omegab) + ) + { + // Note: might need revisiting + dRdF1b = dRdF1b.patchInternalField(); + } + } + + volScalarField dR_dF1_coeffs + ( + zeroFirstCell_*wa() + *( + (gamma1_ - gamma2_)*((2.0/3.0)*omega()*tdivU - tGPrime) + + omega()*omega()*(beta1_ - beta2_) + CDkOmega_ + ) + ); + + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& kb = k().boundaryField()[pI]; + const fvPatchScalarField& omegab = omega().boundaryField()[pI]; + fvPatchScalarField& dRdF1b = dR_dF1_coeffs.boundaryFieldRef()[pI]; + if + ( + isA(kb) + && (isA(omegab) || isA(omegab)) + ) + { + dRdF1b = Zero; + } + else if (isA(kb) && isA(omegab)) + { + dRdF1b = dRdF1b.patchInternalField(); + } + } + + dRdF1 += dR_dF1_coeffs; + return tdRdF1; +} + + +tmp adjointkOmegaSST::dF1_domega +( + const volScalarField& arg1 +) const +{ + return + ( + 4*pow3(arg1)*(scalar(1) - F1_*F1_) + *( + - case_1_F1_*sqrt(k())/(betaStar_*omega()*omega()*y_) + - case_2_F1_*scalar(500)*nu()/sqr(omega()*y_) + + case_3_F1_*scalar(4)*alphaOmega2_*k()/sqr(CDkOmegaPlus_*y_) + *CDkOmega_/omega() + ) + ); +} + + +tmp adjointkOmegaSST::dF1_dGradOmega +( + const volScalarField& arg1 +) const +{ + return + ( + - case_3_F1_*scalar(4)*pow3(arg1)*(scalar(1) - F1_*F1_) + *scalar(8)*k()*sqr(alphaOmega2_/(CDkOmegaPlus_*y_))/omega()*gradK_ + ); +} + + +tmp adjointkOmegaSST::waEqnSourceFromF1() const +{ + tmp arg1 = min + ( + min + ( + max + ( + (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + (4*alphaOmega2_)*k()/(CDkOmegaPlus_*sqr(y_)) + ), + scalar(10) + ); + + volScalarField dR_dF1(this->dR_dF1()); + volScalarField dF1_domega(this->dF1_domega(arg1)); + volVectorField dF1_dGradOmega(this->dF1_dGradOmega(arg1)); + + surfaceScalarField dR_dGradOmega + ( + interpolationScheme("div(dR_dGradOmega)")(). + interpolate(dR_dF1*dF1_dGradOmega) & mesh_.Sf() + ); + + return + ( + dR_dF1*dF1_domega + - fvc::div(dR_dGradOmega) + ); +} + + +tmp adjointkOmegaSST::waEqnSourceFromCDkOmega() const +{ + tmp tsource + ( + 2*zeroFirstCell_*(1 - F1_)*alphaOmega2_/omega()*wa()*gradK_ + ); + volVectorField& source = tsource.ref(); + + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& omegab = omega().boundaryField()[pI]; + fvPatchVectorField& sourceb = source.boundaryFieldRef()[pI]; + if + ( + isA(omegab) + || isA(omegab) + ) + { + sourceb = Zero; + } + else if (isA(omegab)) + { + sourceb = sourceb.patchInternalField(); + } + } + + surfaceScalarField sourceFaces + ( + interpolationScheme("sourceAdjontkOmegaSST")(). + interpolate(source) & mesh_.Sf() + ); + + return + ( + // Differentiation of omega in CDkOmega + fvm::SuSp(zeroFirstCell_*(1. - F1_)*CDkOmega_/omega(), wa()) + // Differentiation of grad(omega) in CDkOmega + + fvc::div(sourceFaces) + ); +} + + +tmp adjointkOmegaSST::kaEqnSourceFromCDkOmega() const +{ + tmp tsource + ( + 2.*zeroFirstCell_*(1. - F1_)*alphaOmega2_/omega()*wa()*gradOmega_ + ); + volVectorField& source = tsource.ref(); + + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& kb = k().boundaryField()[pI]; + fvPatchVectorField& sourceb = source.boundaryFieldRef()[pI]; + if (isA(kb)) + { + sourceb = Zero; + } + else if (isA(kb)) + { + sourceb = sourceb.patchInternalField(); + } + } + + return + fvc::div + ( + interpolationScheme("sourceAdjontkOmegaSST")(). + interpolate(source) & mesh_.Sf() + ); +} + + +tmp adjointkOmegaSST::dF1_dk +( + const volScalarField& arg1 +) const +{ + return + ( + 4*pow3(arg1)*(scalar(1) - F1_*F1_) + *( + case_1_F1_*0.5/(betaStar_*omega()*sqrt(k())*y_) + + case_4_F1_*scalar(4)*alphaOmega2_/(CDkOmegaPlus_*sqr(y_)) + ) + ); +} + + +tmp adjointkOmegaSST::dF1_dGradK +( + const volScalarField& arg1 +) const +{ + return + ( + - case_3_F1_*scalar(4)*pow3(arg1)*(scalar(1) - F1_*F1_) + *scalar(8)*k()*sqr(alphaOmega2_/(CDkOmegaPlus_*y_))/omega() + *gradOmega_ + ); +} + + +tmp adjointkOmegaSST::kaEqnSourceFromF1() const +{ + tmp arg1 = min + ( + min + ( + max + ( + (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + (4*alphaOmega2_)*k()/(CDkOmegaPlus_*sqr(y_)) + ), + scalar(10) + ); + + volScalarField dR_dF1(this->dR_dF1()); + volScalarField dF1_dk(this->dF1_dk(arg1)); + volVectorField dF1_dGradK(this->dF1_dGradK(arg1)); + + surfaceScalarField dR_dGradK + ( + interpolationScheme("div(dR_dGradK)")(). + interpolate(dR_dF1*dF1_dGradK) & mesh_.Sf() + ); + + return + ( + dR_dF1*dF1_dk + - fvc::div(dR_dGradK) + ); +} + + +tmp adjointkOmegaSST::coeffsDifferentiation +( + const volScalarField& primalField, + const volScalarField& adjointField, + const word& schemeName +) const +{ + tmp> scheme + (interpolationScheme(schemeName)); + surfaceScalarField snGradPrimal(fvc::snGrad(primalField)*mesh_.magSf()); + surfaceScalarField flux(scheme().interpolate(adjointField)*snGradPrimal); + + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& primalbf = primalField.boundaryField()[pI]; + if (isA(primalbf)) + { + // Needless, but just to be safe + snGradPrimal.boundaryFieldRef()[pI] = Zero; + flux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(primalbf)) + { + // Note: to be potentially revisited + snGradPrimal.boundaryFieldRef()[pI] = Zero; + flux.boundaryFieldRef()[pI] = Zero; + } + } + + return (fvc::div(flux) - adjointField*fvc::div(snGradPrimal)); +} + + +tmp adjointkOmegaSST::dNutdbMult +( + const volScalarField& primalField, + const volScalarField& adjointField, + const volScalarField& coeffField, + const volScalarField& bcField, + const word& schemeName +) const +{ + tmp> scheme + (interpolationScheme(schemeName)); + surfaceScalarField snGradPrimal(fvc::snGrad(primalField)*mesh_.magSf()); + surfaceScalarField flux(scheme().interpolate(adjointField)*snGradPrimal); + + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& bc = bcField.boundaryField()[pI]; + if (isA(bc)) + { + const fvPatchScalarField& coeffFieldb = + coeffField.boundaryField()[pI]; + snGradPrimal.boundaryFieldRef()[pI] *= + coeffFieldb/coeffFieldb.patchInternalField(); + flux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(bc)) + { + snGradPrimal.boundaryFieldRef()[pI] = Zero; + flux.boundaryFieldRef()[pI] = Zero; + } + } + + return coeffField*(fvc::div(flux) - adjointField*fvc::div(snGradPrimal)); +} + + +tmp adjointkOmegaSST::dNutdbMult +( + const volVectorField& U, + const volVectorField& Ua, + const volScalarField& nut, + const word& schemeName +) const +{ + tmp> scheme + (interpolationScheme(schemeName)); + surfaceVectorField snGradU(fvc::snGrad(U)*mesh_.magSf()); + surfaceScalarField flux(scheme().interpolate(Ua) & snGradU); + + // Terms form the Laplacian part of the momentum stresses + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& bc = nut.boundaryField()[pI]; + if (isA(bc)) + { + flux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(bc)) + { + snGradU.boundaryFieldRef()[pI] = Zero; + flux.boundaryFieldRef()[pI] = Zero; + } + } + + // Terms form the tranpose gradient in the momentum stresses + const surfaceVectorField& Sf = mesh_.Sf(); + surfaceTensorField fluxTranspose + ( + reverseLinear(mesh_).interpolate(Ua)*Sf + ); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& Ub = U.boundaryField()[pI]; + if (!isA(Ub)) + { + const vectorField Uai(Ua.boundaryField()[pI].patchInternalField()); + const vectorField& Sfb = Sf.boundaryField()[pI]; + fluxTranspose.boundaryFieldRef()[pI] = Uai*Sfb; + } + } + volScalarField M(dev2(gradU_) && fvc::div(fluxTranspose)); + const DimensionedField& V = mesh_.V(); + forAll(mesh_.boundary(), pI) + { + const fvPatchScalarField& bc = nut.boundaryField()[pI]; + if (isA(bc)) + { + const vectorField Uai(Ua.boundaryField()[pI].patchInternalField()); + const tensorField dev2GradU + ( + dev2(gradU_.boundaryField()[pI].patchInternalField()) + ); + const vectorField& Sfb = Sf.boundaryField()[pI]; + const labelList& faceCells = mesh_.boundary()[pI].faceCells(); + forAll(faceCells, fI) + { + const label celli = faceCells[fI]; + M[celli] -= ((Uai[fI] & dev2GradU[fI]) & Sfb[fI])/V[celli]; + } + } + } + M.correctBoundaryConditions(); + + //surfaceScalarField fluxTranspose = + // fvc::interpolate(UaGradU, schemeName) & mesh_.Sf(); + //forAll(mesh_.boundary(), pI) + //{ + // const fvPatchScalarField& bc = nut.boundaryField()[pI]; + // const vectorField& Sf = mesh_.boundary()[pI].Sf(); + // if (isA(bc)) + // { + // fluxTranspose.boundaryFieldRef()[pI] = + // ( + // UaGradU.boundaryField()[pI].patchInternalField() + // - ( + // Ua.boundaryField()[pI].patchInternalField() + // & gradU_.boundaryField()[pI] + // ) + // ) & Sf; + // } + // else if (isA(bc)) + // { + // fluxTranspose.boundaryFieldRef()[pI] = + // UaGradU.boundaryField()[pI].patchInternalField() & Sf; + // } + //} + return + fvc::div(flux) - (Ua & fvc::div(snGradU)) + M; + //fvc::div(flux + fluxTranspose) - (Ua & fvc::div(snGradU)); +} + + +tmp adjointkOmegaSST::convectionMeanFlowSource +( + const volScalarField& primalField, + const volScalarField& adjointField +) const +{ + // Grab the interpolation scheme from the primal convection term + tmp> primalInterpolationScheme + ( + convectionScheme(primalField.name()) + ); + + surfaceVectorField interpolatedPrimal + ( + primalInterpolationScheme().interpolate(primalField)*mesh_.Sf() + ); + surfaceVectorField flux + ( + //reverseLinear(mesh_).interpolate(adjointField) + linear(mesh_).interpolate(adjointField) + *interpolatedPrimal + ); + + const volVectorField& U = primalVars_.U(); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& bc = U.boundaryField()[pI]; + if (isA(bc)) + { + flux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(bc)) + { + interpolatedPrimal.boundaryFieldRef()[pI] = Zero; + flux.boundaryFieldRef()[pI] = Zero; + } + } + + return (-fvc::div(flux) + adjointField*fvc::div(interpolatedPrimal)); +} + + +tmp adjointkOmegaSST::GMeanFlowSource +( + tmp& GbyNuMult +) const +{ + surfaceVectorField flux + ( + mesh_.Sf() & reverseLinear(mesh_).interpolate(GbyNuMult()) + ); + + const volVectorField& U = primalVars_.U(); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& bc = U.boundaryField()[pI]; + if (isA(bc)) + { + flux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(bc)) + { + flux.boundaryFieldRef()[pI] = + mesh_.boundary()[pI].Sf() + & GbyNuMult().boundaryField()[pI].patchInternalField(); + } + } + + return fvc::div(flux); +} + + +tmp adjointkOmegaSST::divUMeanFlowSource +( + tmp& divUMult +) const +{ + surfaceVectorField flux + ( + mesh_.Sf()*reverseLinear(mesh_).interpolate(divUMult()) + ); + + const volVectorField& U = primalVars_.U(); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& bc = U.boundaryField()[pI]; + if (isA(bc)) + { + flux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(bc)) + { + flux.boundaryFieldRef()[pI] = + mesh_.boundary()[pI].Sf() + *divUMult().boundaryField()[pI].patchInternalField(); + } + } + + divUMult.clear(); + + return -fvc::div(flux); +} + + +tmp adjointkOmegaSST::diffusionNutMeanFlowMult +( + const volScalarField& primalField, + const volScalarField& adjointField, + const volScalarField& coeffField +) const +{ + // Note: we could grab the snGrad scheme from the diffusion term directly + surfaceScalarField snGradPrimal(fvc::snGrad(primalField)*mesh_.magSf()); + surfaceScalarField flux + ( + reverseLinear(mesh_).interpolate(adjointField)*snGradPrimal + ); + + const volVectorField& U = primalVars_.U(); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& bc = U.boundaryField()[pI]; + if (!isA(bc)) + { + flux.boundaryFieldRef()[pI] = Zero; + snGradPrimal.boundaryFieldRef()[pI] = Zero; + } + } + return (fvc::div(flux) - adjointField*fvc::div(snGradPrimal))*coeffField; +} + + +tmp adjointkOmegaSST::nutMeanFlowSource +( + tmp& mult +) const +{ + volSymmTensorField M + ( + mult*a1_*k()*(1 - case_1_nut_)/(b1_*F2_*S_*S2_)*twoSymm(gradU_) + ); + M.correctBoundaryConditions(); + mult.clear(); + + surfaceVectorField returnFieldFlux + ( + mesh_.Sf() & reverseLinear(mesh_).interpolate(M) + ); + + const volVectorField& U = primalVars_.U(); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& bc = U.boundaryField()[pI]; + if (isA(bc)) + { + returnFieldFlux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(bc)) + { + returnFieldFlux.boundaryFieldRef()[pI] = + mesh_.boundary()[pI].Sf() + & M.boundaryField()[pI].patchInternalField(); + } + } + + // Note: potentially missing contributions form patches with a calculated + // nut bc + + return fvc::div(returnFieldFlux); +} + + +void adjointkOmegaSST::addWallFunctionTerms +( + fvScalarMatrix& kaEqn, + const volScalarField& dR_dnut +) +{ + // Add source term from the differentiation of nutWallFunction + scalarField& source = kaEqn.source(); + const DimensionedField& V = mesh_.V(); + const volScalarField& ka = this->ka(); + const volScalarField& wa = this->wa(); + const volScalarField& k = this->k(); + const volScalarField& omega = this->omega(); + forAll(nutRef().boundaryFieldRef(), patchi) + { + fvPatchScalarField& nutWall = nutRef().boundaryFieldRef()[patchi]; + if (isA(nutWall)) + { + const fvPatch& patch = mesh_.boundary()[patchi]; + const scalarField& magSf = patch.magSf(); + + const autoPtr& turbModel = + primalVars_.turbulence(); + + const scalarField& y = turbModel().y()[patchi]; + const tmp tnuw = turbModel().nu(patchi); + const scalarField& nuw = tnuw(); + + const nutWallFunctionFvPatchScalarField& nutWF = + refCast(nutWall); + const wallFunctionCoefficients& wallCoeffs = nutWF.wallCoeffs(); + const scalar Cmu = wallCoeffs.Cmu(); + const scalar kappa = wallCoeffs.kappa(); + const scalar E = wallCoeffs.E(); + const scalar yPlusLam = wallCoeffs.yPlusLam(); + + const scalar Cmu25 = pow025(Cmu); + + const labelList& faceCells = patch.faceCells(); + const fvPatchScalarField& dR_dnutw = + dR_dnut.boundaryField()[patchi]; + const fvPatchScalarField& omegaw = omega.boundaryField()[patchi]; + bool addTermsFromOmegaWallFuction + (isA(omegaw)); + + const fvPatchVectorField& Uw = + primalVars_.U().boundaryField()[patchi]; + const scalarField magGradUw(mag(Uw.snGrad())); + forAll(nuw, facei) + { + const label celli = faceCells[facei]; + + const scalar sqrtkCell(sqrt(k[celli])); + const scalar yPlus = Cmu25*y[facei]*sqrtkCell/nuw[facei]; + const scalar logEyPlus = log(E*yPlus); + const scalar dnut_dyPlus = + nuw[facei]*kappa*(logEyPlus - 1)/sqr(logEyPlus); + const scalar dyPlus_dk = + Cmu25*y[facei]/(2*nuw[facei]*sqrtkCell); + const scalar dnut_dk = dnut_dyPlus*dyPlus_dk; + + if (yPlusLam < yPlus) + { + // Terms from nutkWallFunction + source[celli] -= dR_dnutw[facei]*dnut_dk*magSf[facei]; + } + if (addTermsFromOmegaWallFuction) + { + const scalar denom(Cmu25*kappa*y[facei]); + const scalar omegaLog(sqrtkCell/denom); + // Terms from omegaLog in omegaWallFunction + source[celli] += + wa[celli]*omegaLog/omega[celli] + /(2*sqrtkCell*denom); + + // Terms from G at the first cell off the wall + source[celli] += + case_1_Pk_[celli]*ka[celli]*V[celli] + *( + (nutWall[facei] + nuw[facei]) + *magGradUw[facei] + *Cmu25 + /(2.0*sqrtkCell*kappa*y[facei]) + ); + + if (yPlusLam < yPlus) + { + source[celli] += + case_1_Pk_[celli]*ka[celli]*V[celli] + *dnut_dk + *magGradUw[facei] + *Cmu25*sqrtkCell + /(kappa*y[facei]); + } + } + } + } + } +} + + +void adjointkOmegaSST::updatePrimalRelatedFields() +{ + if (changedPrimalSolution_) + { + Info<< "Updating primal-based fields of the adjoint turbulence " + << "model ..." << endl; + + const volVectorField& U = primalVars_.U(); + gradU_ = fvc::grad(U); + gradOmega_ = fvc::grad(omega()); + gradK_ = fvc::grad(k()); + + S2_ = 2*magSqr(symm(gradU_)) + + dimensionedScalar(dimless/sqr(dimTime), 1.e-21); + S_ = sqrt(S2_); + GbyNu0_ = gradU_ && dev(twoSymm(gradU_)); + + // Instead of computing G directly here, delegate to RASModelVariables + // to make sure G is computed consistently when the primal fields are + // averaged. The local value computed here is just used to update + // the switch fields + /* + volScalarField G(primalVars_.turbulence()->GName(), nutRef()*GbyNu0_); + omega().correctBoundaryConditions(); + */ + volScalarField G + ( + IOobject + ( + type() + ":G", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimArea/pow3(dimTime), Zero) + ); + G.ref() = primalVars_.RASModelVariables()->G()(); + + CDkOmega_ = + (2*alphaOmega2_)*(gradK_ & gradOmega_)/omega(); + CDkOmegaPlus_ = max + ( + CDkOmega_, + dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10) + ); + F1_ = F1(); + F2_ = F2(); + case_1_Pk_ = neg(G - c1_*betaStar_*k()*omega()); + case_2_Pk_ = pos0(G - c1_*betaStar_*k()*omega()); + case_3_Pk_ = neg(a1_*omega() - b1_*F2_*S_); + + + alphaK_ = alphaK(F1_); + alphaOmega_ = alphaOmega(F1_); + beta_ = beta(F1_); + gamma_ = gamma(F1_); + + // Switch fields for F1 + { + volScalarField arg1_C3 + ( + (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_) + - scalar(500)*nu()/(sqr(y_)*omega()) + ); + volScalarField arg1_C2 + ( + max + ( + (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ) + - (4*alphaOmega2_)*k()/(CDkOmegaPlus_*sqr(y_)) + ); + volScalarField arg1_C1 + ( + min + ( + max + ( + (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + (4*alphaOmega2_)*k()/(CDkOmegaPlus_*sqr(y_)) + ) - scalar(10) + ); + volScalarField CDkOmegaPlus_C1 + ( + CDkOmega_ + - dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10) + ); + + case_1_F1_ = pos(arg1_C3)*neg(arg1_C2)*neg(arg1_C1); + case_2_F1_ = neg0(arg1_C3)*neg(arg1_C2)*neg(arg1_C1); + case_3_F1_ = pos0(arg1_C2)*neg(arg1_C1)*pos(CDkOmegaPlus_C1); + case_4_F1_ = pos0(arg1_C2)*neg(arg1_C1); + } + + // Switch fields for nut + { + volScalarField nut_C1(a1_*omega() - b1_*F2_*S_); + volScalarField arg2_C2 + ( + (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_) + - scalar(500)*nu()/(sqr(y_)*omega()) + ); + volScalarField arg2_C1 + ( + max + ( + (scalar(2)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ) - scalar(100) + ); + + case_1_nut_ = pos(nut_C1); + case_2_nut_ = neg0(nut_C1)*pos(arg2_C2)*neg(arg2_C1); + case_3_nut_ = neg0(nut_C1)*neg0(arg2_C2)*neg(arg2_C1); + } + + { + volScalarField GPrime_C1 + ( + GbyNu0_ + - (c1_/a1_)*betaStar_*omega()*max(a1_*omega(), b1_*F2_*S_) + ); + case_1_GPrime_ = neg(GPrime_C1); + case_2_GPrime_ = pos0(GPrime_C1); + } + + dnut_domega_ = dnut_domega(); + dnut_dk_ = dnut_dk(); + DOmegaEff_ = DomegaEff(F1_); + DkEff_ = DkEff(F1_); + + changedPrimalSolution_ = false; + } +} + + +tmp> adjointkOmegaSST::convectionScheme +( + const word& varName +) const +{ + const surfaceScalarField& phi = primalVars_.phi(); + const surfaceScalarField& phiInst = primalVars_.phiInst(); + word divEntry("div(" + phiInst.name() + ',' + varName +')'); + ITstream& divScheme = mesh_.divScheme(divEntry); + // Skip the first entry which might be 'bounded' or 'Gauss'. + // If it is 'bounded', skip the second entry as well + word discarded(divScheme); + if (discarded == "bounded") + { + discarded = word(divScheme); + } + return surfaceInterpolationScheme::New(mesh_, phi, divScheme); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +adjointkOmegaSST::adjointkOmegaSST +( + incompressibleVars& primalVars, + incompressibleAdjointMeanFlowVars& adjointVars, + objectiveManager& objManager, + const word& adjointTurbulenceModelName, + const word& modelName +) +: + adjointRASModel + ( + modelName, + primalVars, + adjointVars, + objManager, + adjointTurbulenceModelName + ), + + kappa_ + ( + dimensioned::getOrAddToDict + ( + "kappa", + coeffDict_, + 0.41 + ) + ), + alphaK1_ + ( + dimensioned::getOrAddToDict + ( + "alphaK1", + this->coeffDict_, + 0.85 + ) + ), + alphaK2_ + ( + dimensioned::getOrAddToDict + ( + "alphaK2", + this->coeffDict_, + 1.0 + ) + ), + alphaOmega1_ + ( + dimensioned::getOrAddToDict + ( + "alphaOmega1", + this->coeffDict_, + 0.5 + ) + ), + alphaOmega2_ + ( + dimensioned::getOrAddToDict + ( + "alphaOmega2", + this->coeffDict_, + 0.856 + ) + ), + gamma1_ + ( + dimensioned::getOrAddToDict + ( + "gamma1", + this->coeffDict_, + 5.0/9.0 + ) + ), + gamma2_ + ( + dimensioned::getOrAddToDict + ( + "gamma2", + this->coeffDict_, + 0.44 + ) + ), + beta1_ + ( + dimensioned::getOrAddToDict + ( + "beta1", + this->coeffDict_, + 0.075 + ) + ), + beta2_ + ( + dimensioned::getOrAddToDict + ( + "beta2", + this->coeffDict_, + 0.0828 + ) + ), + betaStar_ + ( + dimensioned::getOrAddToDict + ( + "betaStar", + this->coeffDict_, + 0.09 + ) + ), + a1_ + ( + dimensioned::lookupOrAddToDict + ( + "a1", + this->coeffDict_, + 0.31 + ) + ), + b1_ + ( + dimensioned::lookupOrAddToDict + ( + "b1", + this->coeffDict_, + 1.0 + ) + ), + c1_ + ( + dimensioned::lookupOrAddToDict + ( + "c1", + this->coeffDict_, + 10.0 + ) + ), + F3_ + ( + Switch::lookupOrAddToDict + ( + "F3", + this->coeffDict_, + false + ) + ), + + y_(primalVars_.RASModelVariables()().d()), + + //Primal Gradient Fields + gradU_ + ( + IOobject + ( + "rasModel::gradU", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedTensor(dimless/dimTime, Zero) + ), + gradOmega_ + ( + IOobject + ( + "rasModel::gradOmega", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector(dimless/dimTime/dimLength, Zero) + ), + gradK_ + ( + IOobject + ( + "rasModel::gradK", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector(dimLength/sqr(dimTime), Zero) + ), + + S2_ + ( + IOobject + ( + "S2", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless/sqr(dimTime), Zero) + ), + S_ + ( + IOobject + ( + "kOmegaSST_S", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless/dimTime, Zero) + ), + GbyNu0_ + ( + IOobject + ( + "adjointRASModel::GbyNu0", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless/sqr(dimTime), Zero) + ), + CDkOmega_ + ( + IOobject + ( + "CDkOmega_", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless/sqr(dimTime), Zero) + ), + CDkOmegaPlus_ + ( + IOobject + ( + "CDkOmegaPlus", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless/sqr(dimTime), Zero) + ), + F1_ + ( + IOobject + ( + "F1", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + F2_ + ( + IOobject + ( + "F2", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + // Model Field coefficients + alphaK_ + ( + IOobject + ( + "alphaK", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + alphaOmega_ + ( + IOobject + ( + "alphaOmega", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + beta_ + ( + IOobject + ( + "beta", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + gamma_ + ( + IOobject + ( + "gamma", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + + case_1_F1_ + ( + IOobject + ( + "case_1_F1", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_2_F1_ + ( + IOobject + ( + "case_2_F1", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_3_F1_ + ( + IOobject + ( + "case_3_F1", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_4_F1_ + ( + IOobject + ( + "case_4_F1", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_1_Pk_ + ( + IOobject + ( + "case_1_Pk", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_2_Pk_ + ( + IOobject + ( + "case_2_Pk", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_3_Pk_ + ( + IOobject + ( + "case_3_Pk", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + + case_1_nut_ + ( + IOobject + ( + "case_1_nut", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_2_nut_ + ( + IOobject + ( + "case_2_nut", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_3_nut_ + ( + IOobject + ( + "case_3_nut", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_1_GPrime_ + ( + IOobject + ( + "case_1_GPrime", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + case_2_GPrime_ + ( + IOobject + ( + "case_2_GPrime", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimless, Zero) + ), + + // Zero 1rst cell field + firstCellIDs_(0), + zeroFirstCell_(zeroFirstCell()), + + // Turbulence model multipliers + dnut_domega_ + ( + IOobject + ( + "dnut_domega", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(sqr(dimLength), Zero) + ), + dnut_dk_ + ( + IOobject + ( + "dnut_dk", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimTime, Zero) + ), + DOmegaEff_ + ( + IOobject + ( + "DomegaEff", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(nutRef().dimensions(), Zero) + ), + DkEff_ + ( + IOobject + ( + "DkEff", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(nutRef().dimensions(), Zero) + ) +{ + adjointTMVariablesBaseNames_.setSize(2); + adjointTMVariablesBaseNames_[0] = "ka"; + adjointTMVariablesBaseNames_[1] = "wa"; + // Read in adjoint fields + variablesSet::setField + ( + adjointTMVariable1Ptr_, + mesh_, + "ka", + adjointVars.solverName(), + adjointVars.useSolverNameForFields() + ); + variablesSet::setField + ( + adjointTMVariable2Ptr_, + mesh_, + "wa", + adjointVars.solverName(), + adjointVars.useSolverNameForFields() + ); + + setMeanFields(); + + // No sensitivity contributions from the adjoint to the eikonal equation + // for the moment + includeDistance_ = false; + + // Update the primal related fields here so that functions computing + // sensitivities have the updated fields in case of continuation + updatePrimalRelatedFields(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +tmp adjointkOmegaSST::devReff() const +{ + const volVectorField& Ua = adjointVars_.UaInst(); + return devReff(Ua); +} + + +tmp adjointkOmegaSST::devReff +( + const volVectorField& U +) const +{ + return + tmp::New + ( + IOobject + ( + "devRhoReff", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + -nuEff()*dev(twoSymm(fvc::grad(U))) + ); +} + + +tmp adjointkOmegaSST::divDevReff(volVectorField& Ua) const +{ + tmp tnuEff = nuEff(); + const volScalarField& nuEff = tnuEff(); + + return + ( + - fvm::laplacian(nuEff, Ua) + - fvc::div(nuEff*dev(fvc::grad(Ua)().T())) + ); + + /* WIP + const volVectorField& U = primalVars_.U(); + const surfaceVectorField& Sf = mesh_.Sf(); + tmp tflux = + reverseLinear(mesh_).interpolate(Ua)*Sf; + surfaceTensorField& flux = tflux.ref(); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& Ub = U.boundaryField()[pI]; + if (!isA(Ub)) + { + const vectorField Uai = Ua.boundaryField()[pI].patchInternalField(); + const vectorField& Sfb = Sf.boundaryField()[pI]; + flux.boundaryFieldRef()[pI] = Uai*Sfb; + } + } + volTensorField M(nuEff*dev2(fvc::div(flux))); + const DimensionedField& V = mesh_.V(); + + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& Ub = U.boundaryField()[pI]; + if (!isA(Ub)) + { + const fvPatchScalarField& nuEffb = nuEff.boundaryField()[pI]; + const vectorField nf = mesh_.boundary()[pI].nf(); + const vectorField Uai = Ua.boundaryField()[pI].patchInternalField(); + const labelList& faceCells = mesh_.boundary()[pI].faceCells(); + const vectorField& Sfb = Sf.boundaryField()[pI]; + + forAll(faceCells, fI) + { + const label celli = faceCells[fI]; + const tensor t(dev2(Uai[fI]*Sfb[fI])); + M[celli] -= nuEffb[fI]*(t - nf[fI]*(nf[fI] & t))/V[celli]; + } + } + } + M.correctBoundaryConditions(); + + surfaceVectorField returnFlux = + - (Sf & reverseLinear(mesh_).interpolate(M)); + forAll(mesh_.boundary(), pI) + { + const fvPatchVectorField& Ub = U.boundaryField()[pI]; + if (isA(Ub)) + { + returnFlux.boundaryFieldRef()[pI] = Zero; + } + else if (isA(Ub)) + { + const scalarField& deltaCoeffs = mesh_.boundary()[pI].deltaCoeffs(); + const fvPatchScalarField& nuEffb = nuEff.boundaryField()[pI]; + const vectorField Uai = Ua.boundaryField()[pI].patchInternalField(); + const vectorField nf = mesh_.boundary()[pI].nf(); + const vectorField& Sfb = Sf.boundaryField()[pI]; + + returnFlux.boundaryFieldRef()[pI] = + - (Sfb & M.boundaryField()[pI].patchInternalField()) + + nuEffb*deltaCoeffs*(nf & dev2(Uai*Sfb)); + } + } + + return + ( + - fvm::laplacian(nuEff, Ua) + + fvc::div(returnFlux) + ); + */ +} + + +tmp adjointkOmegaSST::nonConservativeMomentumSource() const +{ + return (ka()*gradK_ + wa()*gradOmega_); +} + + +tmp adjointkOmegaSST::adjointMeanFlowSource() +{ + tmp tmeanFlowSource + ( + tmp::New + ( + IOobject + ( + "adjointMeanFlowSource" + type(), + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector(dimVelocity/dimTime, Zero) + ) + ); + volVectorField& meanFlowSource = tmeanFlowSource.ref(); + + // Contributions from the convection terms of the turbulence model + meanFlowSource += + convectionMeanFlowSource(omega(), zeroFirstCell_*wa()) + + convectionMeanFlowSource(k(), ka()); + + // Contributions from GbyNu, including gradU + tmp twoSymmGradU(twoSymm(gradU_)); + tmp GbyNuMult + ( + // First part of GPrime and G from Pk + 2.*dev(twoSymmGradU())*zeroFirstCell_ + *(wa()*gamma_*case_1_GPrime_ + ka()*nutRef()*case_1_Pk_) + // Second part of GPrime + + twoSymmGradU()*wa()*zeroFirstCell_*gamma_*case_2_GPrime_ + *(1. - case_1_nut_)*c1_/a1_*betaStar_*omega()*b1_*F2_/S_ + ); + twoSymmGradU.clear(); + meanFlowSource += GMeanFlowSource(GbyNuMult); + GbyNuMult.clear(); + + // Contributions from divU + tmp divUMult + ( + (2.0/3.0)*(zeroFirstCell_*wa()*omega()*gamma_ + ka()*k()) + ); + meanFlowSource += divUMeanFlowSource(divUMult); + + // Contributions from S2, existing in nut + const volVectorField& U = primalVars_.U(); + const volVectorField& Ua = adjointVars_.UaInst(); + tmp nutMeanFlowSourceMult + ( + // nut in the diffusion coefficients + diffusionNutMeanFlowMult(k(), ka(), alphaK_) + + diffusionNutMeanFlowMult(omega(), zeroFirstCell_*wa(), alphaOmega_) + + dNutdbMult(U, Ua, nutRef(), "coeffsDiff") + // nut in G + - ka()*case_1_Pk_*zeroFirstCell_*GbyNu0_ + ); + meanFlowSource += nutMeanFlowSource(nutMeanFlowSourceMult); + + // G at the first cell includes mag(U.snGrad()) + // Add term here + forAll(omega().boundaryFieldRef(), patchi) + { + fvPatchScalarField& omegaWall = omega().boundaryFieldRef()[patchi]; + if (isA(omegaWall)) + { + const fvPatch& patch = mesh_.boundary()[patchi]; + + const autoPtr& turbModel = + primalVars_.turbulence(); + + const scalarField& y = turbModel().y()[patchi]; + const tmp tnuw = turbModel().nu(patchi); + const scalarField& nuw = tnuw(); + + const nutWallFunctionFvPatchScalarField& nutw = + refCast + (nutRef().boundaryFieldRef()[patchi]); + const wallFunctionCoefficients& wallCoeffs = nutw.wallCoeffs(); + const scalar Cmu = wallCoeffs.Cmu(); + const scalar kappa = wallCoeffs.kappa(); + const scalar Cmu25 = pow025(Cmu); + + const labelList& faceCells = patch.faceCells(); + + const fvPatchVectorField& Uw = + primalVars_.U().boundaryField()[patchi]; + vectorField snGradUw(Uw.snGrad()); + const scalarField& deltaCoeffs = patch.deltaCoeffs(); + forAll(faceCells, facei) + { + const label celli = faceCells[facei]; + // Volume will be added when meanFlowSource is added to UaEqn + meanFlowSource[celli] += + ka()[celli]*case_1_Pk_[celli] + *(nutw[facei] + nuw[facei]) + *snGradUw[facei].normalise() + *Cmu25*sqrt(k()[celli]) + *deltaCoeffs[facei] + /(kappa*y[facei]); + } + } + } + return tmeanFlowSource; +} + + +tmp adjointkOmegaSST::nutJacobianTMVar1() const +{ + return dnut_dk_; +} + + +tmp adjointkOmegaSST::nutJacobianTMVar2() const +{ + return dnut_domega_; +} + + +tmp adjointkOmegaSST::diffusionCoeffVar1(label patchI) const +{ + return + ( + alphaK_.boundaryField()[patchI]*nutRef().boundaryField()[patchI] + + nu()().boundaryField()[patchI] + ); +} + + +tmp adjointkOmegaSST::diffusionCoeffVar2(label patchI) const +{ + return + ( + alphaOmega_.boundaryField()[patchI]*nutRef().boundaryField()[patchI] + + nu()().boundaryField()[patchI] + ); +} + + +void adjointkOmegaSST::correct() +{ + adjointRASModel::correct(); + + if (!adjointTurbulence_) + { + return; + } + + updatePrimalRelatedFields(); + + // Primal and adjoint fields + const volVectorField& U = primalVars_.U(); + const surfaceScalarField& phi = primalVars_.phi(); + volScalarField dR_dnut(this->dR_dnut()); + volScalarField::Internal divU(fvc::div(fvc::absolute(phi, U))); + + tmp waEqn + ( + fvm::div(-phi, wa()) + + fvm::SuSp(zeroFirstCell_*fvc::div(phi), wa()) + - fvm::laplacian(DOmegaEff_, wa()) + + waEqnSourceFromCDkOmega() + + waEqnSourceFromF1() + + dR_dnut*dnut_domega_ + + fvm::Sp(zeroFirstCell_*scalar(2.)*beta_*omega(), wa()) + + fvm::SuSp + ( + zeroFirstCell_()*gamma_*((2.0/3.0)*divU - dGPrime_domega().ref()()), + wa() + ) + - (case_2_Pk_*c1_ - scalar(1))*betaStar_*k()*ka() + ); + + // Boundary manipulate changes the diagonal component, so relax has to + // come after that + waEqn.ref().boundaryManipulate(wa().boundaryFieldRef()); + waEqn.ref().relax(); + + // Sources from the objective should be added after the boundary + // manipulation + objectiveManager_.addTMEqn2Source(waEqn.ref()); + waEqn.ref().solve(); + + // Adjoint Turbulent kinetic energy equation + tmp kaEqn + ( + fvm::div(-phi, ka()) + + fvm::SuSp(fvc::div(phi), ka()) + - fvm::laplacian(DkEff_, ka()) + + fvm::Sp(betaStar_*omega(), ka()) + - case_2_Pk_()*c1_*betaStar_*omega()()*ka()() + + fvm::SuSp(scalar(2.0/3.0)*divU, ka()) + + kaEqnSourceFromCDkOmega() + + kaEqnSourceFromF1() + + dR_dnut()*dnut_dk_() + - zeroFirstCell_()*gamma_*dGPrime_dk().ref()()*wa() + ); + + kaEqn.ref().relax(); + kaEqn.ref().boundaryManipulate(ka().boundaryFieldRef()); + addWallFunctionTerms(kaEqn.ref(), dR_dnut); + // Add sources from the objective functions + objectiveManager_.addTMEqn1Source(kaEqn.ref()); + + kaEqn.ref().solve(); + + if (adjointVars_.getSolverControl().printMaxMags()) + { + dimensionedScalar maxwa = max(mag(wa())); + dimensionedScalar maxka = max(mag(ka())); + Info<< "Max mag of adjoint dissipation = " << maxwa.value() << endl; + Info<< "Max mag of adjoint kinetic energy = " << maxka.value() << endl; + } +} + + +const boundaryVectorField& adjointkOmegaSST::adjointMomentumBCSource() const +{ + return adjMomentumBCSourcePtr_(); +} + + +const boundaryVectorField& adjointkOmegaSST::wallShapeSensitivities() +{ + boundaryVectorField& wallShapeSens = wallShapeSensitivitiesPtr_(); + volTensorField FITerm(FISensitivityTerm()); + + forAll(mesh_.boundary(), patchi) + { + vectorField nf(mesh_.boundary()[patchi].nf()); + wallShapeSens[patchi] = nf & FITerm.boundaryField()[patchi]; + } + return wallShapeSens; +} + + +const boundaryVectorField& adjointkOmegaSST::wallFloCoSensitivities() +{ + return wallFloCoSensitivitiesPtr_(); +} + + +tmp adjointkOmegaSST::distanceSensitivities() +{ + return tmp::New + ( + IOobject + ( + "adjointEikonalSource" + type(), + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar(dimLength/pow3(dimTime), Zero) + ); +} + + +tmp adjointkOmegaSST::FISensitivityTerm() +{ + const volVectorField& U = primalVars_.U(); + const volScalarField& kInst = + primalVars_.RASModelVariables()->TMVar1Inst(); + const volScalarField& omegaInst = + primalVars_.RASModelVariables()->TMVar2Inst(); + + tmp arg1 = min + ( + min + ( + max + ( + (scalar(1)/betaStar_)*sqrt(k())/(omega()*y_), + scalar(500)*nu()/(sqr(y_)*omega()) + ), + (4*alphaOmega2_)*k()/(CDkOmegaPlus_*sqr(y_)) + ), + scalar(10) + ); + + // Interpolation schemes used by the primal convection terms + auto kScheme(convectionScheme(kInst.name())); + auto omegaScheme(convectionScheme(omegaInst.name())); + const surfaceVectorField& Sf = mesh_.Sf(); + tmp tFISens + ( + tmp::New + ( + IOobject + ( + type() + "FISensTerm", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedTensor(sqr(dimLength)/pow3(dimTime), Zero), + zeroGradientFvPatchTensorField::typeName + ) + ); + volTensorField& FISens = tFISens.ref(); + FISens = + // k convection + - ka()*fvc::div + ( + kScheme().interpolate(k()) + *linear(mesh_).interpolate(U)*Sf + ) + // k diffusion + + ka()*T(fvc::grad(DkEff_*gradK_)) + - DkEff_*(fvc::grad(ka())*gradK_) + // omega convection + - wa()*zeroFirstCell_*fvc::div + ( + omegaScheme().interpolate(omega()) + *linear(mesh_).interpolate(U)*Sf + ) + // omega diffusion + + wa()*zeroFirstCell_*T(fvc::grad(DOmegaEff_*gradOmega_)) + - DOmegaEff_*(fvc::grad(wa()*zeroFirstCell_)*gradOmega_) + // terms including GbyNu0 + + ( + case_1_GPrime_*wa()*gamma_ + + case_1_Pk_*ka()*nutRef() + )*2.*T(gradU_ & dev(twoSymm(gradU_)))*zeroFirstCell_ + // S2 (includes contribution from nut in UEqn as well) + + ( + dR_dnut()*a1_*k()/(b1_*S_*S_*S_*F2_) + + wa()*zeroFirstCell_*gamma_*case_2_GPrime_ + *(c1_/a1_)*betaStar_*omega()*b1_*F2_/S_ + )*T(gradU_ & twoSymm(gradU_))*(1. - case_1_nut_) + // CDkOmega in omegaEqn + + 2.*wa()*(1. - F1_)*alphaOmega2_/omega()*zeroFirstCell_ + *(gradOmega_*gradK_ + gradK_*gradOmega_) + // F1 + - dR_dF1() + *(dF1_dGradK(arg1)*gradK_ + dF1_dGradOmega(arg1)*gradOmega_); + + FISens.correctBoundaryConditions(); + + return tFISens; +} + + +tmp adjointkOmegaSST::topologySensitivities +( + const word& designVarsName +) const +{ + // Missing proper terms - return zero for now + return tmp::New(mesh_.nCells(), Zero); +} + + +void adjointkOmegaSST::nullify() +{ + variablesSet::nullifyField(ka()); + variablesSet::nullifyField(wa()); +} + + +bool adjointkOmegaSST::read() +{ + if (adjointRASModel::read()) + { + kappa_.readIfPresent(coeffDict()); + alphaK1_.readIfPresent(this->coeffDict()); + alphaK2_.readIfPresent(this->coeffDict()); + alphaOmega1_.readIfPresent(this->coeffDict()); + alphaOmega2_.readIfPresent(this->coeffDict()); + gamma1_.readIfPresent(this->coeffDict()); + gamma2_.readIfPresent(this->coeffDict()); + beta1_.readIfPresent(this->coeffDict()); + beta2_.readIfPresent(this->coeffDict()); + betaStar_.readIfPresent(this->coeffDict()); + a1_.readIfPresent(this->coeffDict()); + b1_.readIfPresent(this->coeffDict()); + c1_.readIfPresent(this->coeffDict()); + F3_.readIfPresent("F3", this->coeffDict()); + + return true; + } + else + { + return false; + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace adjointRASModels +} // End namespace incompressible +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.H new file mode 100644 index 0000000000..f670fa7369 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSST.H @@ -0,0 +1,620 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::incompressibleAdjoint::adjointRASModels::adjointkOmegaSST + +Description + Continuous adjoint to the kOmegaSST turbulence model for incompressible + flows. + + Reference: + \verbatim + The code is based on the following reference, with a number of + changes in the numerical implementation + + Kavvadias, I., Papoutsis-Kiachagias, E., Dimitrakopoulos, G., & + Giannakoglou, K. (2014). + The continuous adjoint approach to the k–ω SST turbulence model + with applications in shape optimization + Engineering Optimization, 47(11), 1523-1542. + https://doi.org/10.1080/0305215X.2014.979816 + \endverbatim + +SourceFiles + adjointkOmegaSST.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointkOmegaSST_H +#define adjointkOmegaSST_H + +#include "adjointRASModel.H" +#include "wallDist.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressibleAdjoint +{ +namespace adjointRASModels +{ + +/*---------------------------------------------------------------------------*\ + Class adjointkOmegaSST Declaration +\*---------------------------------------------------------------------------*/ + +class adjointkOmegaSST +: + public adjointRASModel +{ + // Private Member Functions + + //- No copy construct + adjointkOmegaSST(const adjointkOmegaSST&) = delete; + + //- No copy assignment + void operator=(const adjointkOmegaSST&) = delete; + + +protected: + + // Protected data + + // Primal Model coefficients + + dimensionedScalar kappa_; + dimensionedScalar alphaK1_; + dimensionedScalar alphaK2_; + + dimensionedScalar alphaOmega1_; + dimensionedScalar alphaOmega2_; + + dimensionedScalar gamma1_; + dimensionedScalar gamma2_; + + dimensionedScalar beta1_; + dimensionedScalar beta2_; + + dimensionedScalar betaStar_; + + dimensionedScalar a1_; + dimensionedScalar b1_; + dimensionedScalar c1_; + + //- Flag to include the F3 term + Switch F3_; + + + // Fields + + // Primal Fields + + //- Wall distance + // Note: reference to the distance known by the primal model + const volScalarField& y_; + + //- Cached primal gradient fields + volTensorField gradU_; + volVectorField gradOmega_; + volVectorField gradK_; + + //- Primal cached fields involved in the solution of the + // adjoint equations + // Cached to reduce the computational cost + volScalarField S2_; + volScalarField S_; + volScalarField GbyNu0_; + volScalarField CDkOmega_; + volScalarField CDkOmegaPlus_; + volScalarField F1_; + volScalarField F2_; + + // Model Field coefficients + volScalarField alphaK_; + volScalarField alphaOmega_; + volScalarField beta_; + volScalarField gamma_; + + + // Switch fields + + // Switch fields for the differentiation of F1 + volScalarField case_1_F1_; + volScalarField case_2_F1_; + volScalarField case_3_F1_; + volScalarField case_4_F1_; + + //- Switch fields for the production in the k Eqn + volScalarField case_1_Pk_; + volScalarField case_2_Pk_; + volScalarField case_3_Pk_; + + // Switch fields for the differentiation of nut + // Holds also for the differentiation of the second branch of + // GbyNu + volScalarField case_1_nut_; + volScalarField case_2_nut_; + volScalarField case_3_nut_; + + // Switch fields for GPrime + volScalarField case_1_GPrime_; + volScalarField case_2_GPrime_; + + // Zero first cell field and IDs + // Since the omega equation is a two-zonal one, some + // of the terms in the adjoint equations need to ba canceled + // at the cells next to omegaWallFunction BCs + labelList firstCellIDs_; + volScalarField zeroFirstCell_; + + + // Turbulence model multipliers + + //- Nut Jacobian w.r.t. omega + volScalarField dnut_domega_; + + //- Nut Jacobian w.r.t. k + volScalarField dnut_dk_; + + //- Diffusivity of the omega equation + volScalarField DOmegaEff_; + + //- Diffusivity of the k equation + volScalarField DkEff_; + + + // Protected Member Functions + + // Primal functions + + virtual tmp F1() const; + virtual tmp F2() const; + virtual tmp GbyNu + ( + const volScalarField& GbyNu0, + const volScalarField& F2, + const volScalarField& S2 + ) const; + + //- Return G/nu + virtual tmp GbyNu + ( + const volScalarField::Internal& GbyNu0, + const volScalarField::Internal& F2, + const volScalarField::Internal& S2 + ) const; + + tmp blend + ( + const volScalarField& F1, + const dimensionedScalar& psi1, + const dimensionedScalar& psi2 + ) const + { + return F1*(psi1 - psi2) + psi2; + } + + tmp blend + ( + const volScalarField::Internal& F1, + const dimensionedScalar& psi1, + const dimensionedScalar& psi2 + ) const + { + return F1*(psi1 - psi2) + psi2; + } + + tmp alphaK(const volScalarField& F1) const + { + return blend(F1, alphaK1_, alphaK2_); + } + + tmp alphaOmega(const volScalarField& F1) const + { + return blend(F1, alphaOmega1_, alphaOmega2_); + } + + tmp beta + ( + const volScalarField::Internal& F1 + ) const + { + return tmp::New + ( + this->type() + ":beta", + blend(F1, beta1_, beta2_) + ); + } + + tmp beta + ( + const volScalarField& F1 + ) const + { + return tmp::New + ( + this->type() + ":beta", + blend(F1, beta1_, beta2_) + ); + } + + tmp gamma + ( + const volScalarField::Internal& F1 + ) const + { + return tmp::New + ( + this->type() + ":gamma", + blend(F1, gamma1_, gamma2_) + ); + } + + tmp gamma + ( + const volScalarField& F1 + ) const + { + return tmp::New + ( + this->type() + ":gamma", + blend(F1, gamma1_, gamma2_) + ); + } + + tmp zeroFirstCell(); + + + // References to the primal turbulence model variables + + inline const volScalarField& k() const + { + return primalVars_.RASModelVariables()().TMVar1(); + } + + inline volScalarField& k() + { + return primalVars_.RASModelVariables()().TMVar1(); + } + + inline const volScalarField& omega() const + { + return primalVars_.RASModelVariables()().TMVar2(); + } + + inline volScalarField& omega() + { + return primalVars_.RASModelVariables()().TMVar2(); + } + + inline const volScalarField& nutRef() const + { + return primalVars_.RASModelVariables()().nutRef(); + } + + inline volScalarField& nutRef() + { + return primalVars_.RASModelVariables()().nutRef(); + } + + + // Adjoint related protected member functions + + //- Derivative of the primal equations wrt nut + tmp dR_dnut(); + + //- Nut Jacobian wrt omega + tmp dnut_domega() const; + + //- Nut Jacobian wrt k + tmp dnut_dk() const; + + //- F2 Jacobian wrt omega + tmp dF2_domega() const; + + //- F2 Jacobian wrt k + tmp dF2_dk() const; + + //- GbyNu Jacobian wrt omega + tmp dGPrime_domega() const; + + //- GbyNu Jacobian wrt k + tmp dGPrime_dk() const; + + //- Derivative of the primal equations wrt F1 + tmp dR_dF1() const; + + //- F1 Jacobian wrt omega (no contributions from grad(omega)) + tmp dF1_domega(const volScalarField& arg1) const; + + //- F1 Jacobian wrt grad(omega) + tmp dF1_dGradOmega + ( + const volScalarField& arg1 + ) const; + + //- Source to waEqn from the differentiation of F1 + tmp waEqnSourceFromF1() const; + + //- Source to waEqn from the differentiation of CDkOmega + tmp waEqnSourceFromCDkOmega() const; + + //- F1 Jacobian wrt k (no contributions from grad(k)) + tmp dF1_dk(const volScalarField& arg1) const; + + //- F1 Jacobian wrt grad(k) + tmp dF1_dGradK(const volScalarField& arg1) const; + + //- Source to kaEqn from the differentiation of F1 + tmp kaEqnSourceFromF1() const; + + //- Source to kaEqn from the differentiation of CDkOmega + tmp kaEqnSourceFromCDkOmega() const; + + //- Differentiation of the turbulence model diffusion coefficients + tmp coeffsDifferentiation + ( + const volScalarField& primalField, + const volScalarField& adjointField, + const word& schemeName + ) const; + + //- Term multiplying dnut/db, coming from the turbulence model + tmp dNutdbMult + ( + const volScalarField& primalField, + const volScalarField& adjointField, + const volScalarField& coeffField, + const volScalarField& bcField, + const word& schemeName + ) const; + + //- Term multiplying dnut/db, coming from the momentum equations + tmp dNutdbMult + ( + const volVectorField& primalField, + const volVectorField& adjointField, + const volScalarField& bcField, + const word& schemeName + ) const; + + // Functions computing the adjoint mean flow source + + //- Contributions from the turbulence model convection terms + tmp convectionMeanFlowSource + ( + const volScalarField& primalField, + const volScalarField& adjointField + ) const; + + //- Contributions from the G + tmp GMeanFlowSource + ( + tmp& GbyNuMult + ) const; + + //- Contributions from the divU + tmp divUMeanFlowSource + ( + tmp& divUMult + ) const; + + //- Contributions from nut(U), in the diffusion coefficients + //- of the turbulence model + tmp diffusionNutMeanFlowMult + ( + const volScalarField& primalField, + const volScalarField& adjointField, + const volScalarField& coeffField + ) const; + + //- Contributions from nut(U) + tmp nutMeanFlowSource + ( + tmp& mult + ) const; + + + //- Contributions from the differentiation of k existing in + //- nutkWallFunction. + // This could also be implemented in kaqRWallFunction but all + // the fields required for the computation already exist here, + // hence the code complexity is reduced + void addWallFunctionTerms + ( + fvScalarMatrix& kaEqn, + const volScalarField& dR_dnut + ); + + // References to the adjoint turbulence model fields + + inline volScalarField& ka() + { + return adjointTMVariable1Ptr_(); + }; + + inline const volScalarField& ka() const + { + return adjointTMVariable1Ptr_(); + }; + + inline volScalarField& wa() + { + return adjointTMVariable2Ptr_(); + }; + + inline const volScalarField& wa() const + { + return adjointTMVariable2Ptr_(); + }; + + + //- Update of the primal cached fields + void updatePrimalRelatedFields(); + + //- Return the requested interpolation scheme if it exists, + //- otherwise return a reverseLinear scheme + template + tmp> interpolationScheme + ( + const word& schemeName + ) const; + + //- Return the interpolation scheme used by the primal convection + //- term of the equation corresponding to the argument + tmp> convectionScheme + ( + const word& varName + ) const; + + +public: + + //- Runtime type information + TypeName("adjointkOmegaSST"); + + + // Constructors + + //- Construct from components + adjointkOmegaSST + ( + incompressibleVars& primalVars, + incompressibleAdjointMeanFlowVars& adjointVars, + objectiveManager& objManager, + const word& adjointTurbulenceModelName + = adjointTurbulenceModel::typeName, + const word& modelName = typeName + ); + + + //- Destructor + virtual ~adjointkOmegaSST() = default; + + + // Member Functions + + //- Return the effective diffusivity for k + tmp DkEff(const volScalarField& F1) const + { + return tmp + ( + new volScalarField("DkEff", alphaK(F1)*this->nut() + this->nu()) + ); + } + + //- Return the effective diffusivity for omega + tmp DomegaEff(const volScalarField& F1) const + { + return tmp + ( + new volScalarField + ( + "DomegaEff", + alphaOmega(F1)*this->nut() + this->nu() + ) + ); + } + + virtual tmp devReff() const; + + virtual tmp devReff(const volVectorField& U) const; + + //- Return the transpose part of the adjoint momentum stresses + virtual tmp divDevReff(volVectorField& U) const; + + //- Non-conservative part of the terms added to the mean flow equations + virtual tmp nonConservativeMomentumSource() const; + + //- Source term added to the adjoint mean flow due to the + //- differentiation of the turbulence model + virtual tmp adjointMeanFlowSource(); + + //- Jacobian of nut wrt to k + // Needs to be implemented for objectives related to nut, defined in + // the internal field + virtual tmp nutJacobianTMVar1() const; + + //- Jacobian of nut wrt to omega + // Needs to be implemented for objectives related to nut, defined in + // the internal field + virtual tmp nutJacobianTMVar2() const; + + //- Diffusion coeff at the boundary for k + virtual tmp diffusionCoeffVar1(label patchI) const; + + //- Diffusion coeff at the boundary for omega + virtual tmp diffusionCoeffVar2(label patchI) const; + + virtual const boundaryVectorField& adjointMomentumBCSource() const; + + //- Sensitivity derivative contributions when using the (E)SI approach + virtual const boundaryVectorField& wallShapeSensitivities(); + + virtual const boundaryVectorField& wallFloCoSensitivities(); + + //- Contributions to the adjoint eikonal equation (zero for now) + virtual tmp distanceSensitivities(); + + //- Sensitivity derivative contributions when using the FI approach + virtual tmp FISensitivityTerm(); + + virtual tmp topologySensitivities + ( + const word& designVarsName + ) const; + + //- Nullify all adjoint turbulence model fields and their old times + virtual void nullify(); + + //- Solve the adjoint turbulence equations + virtual void correct(); + + //- Read adjointRASProperties dictionary + virtual bool read(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace adjointRASModels +} // End namespace incompressibleAdjoint +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "adjointkOmegaSSTTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C new file mode 100644 index 0000000000..467cd8b01f --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointkOmegaSST/adjointkOmegaSSTTemplates.C @@ -0,0 +1,69 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "reverseLinear.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressibleAdjoint +{ +namespace adjointRASModels +{ + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template +tmp> adjointkOmegaSST::interpolationScheme +( + const word& schemeName +) const +{ + return + tmp> + ( + mesh_.interpolationSchemes().found(schemeName) ? + surfaceInterpolationScheme::New + ( + mesh_, + mesh_.interpolationScheme(schemeName) + ) : + tmp> + (new reverseLinear(mesh_)) + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace adjointRASModels +} // End namespace incompressible +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.C new file mode 100644 index 0000000000..85675f0e35 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.C @@ -0,0 +1,178 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "adjointFarFieldTMVar1FvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" +#include "surfaceFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +adjointFarFieldTMVar1FvPatchScalarField:: +adjointFarFieldTMVar1FvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, word::null) +{} + + +adjointFarFieldTMVar1FvPatchScalarField:: +adjointFarFieldTMVar1FvPatchScalarField +( + const adjointFarFieldTMVar1FvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchScalarField(ptf, p, iF, mapper), + adjointScalarBoundaryCondition(p, iF, ptf.adjointSolverName_) +{} + + +adjointFarFieldTMVar1FvPatchScalarField:: +adjointFarFieldTMVar1FvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, dict.get("solverName")) +{ + fvPatchField::operator= + ( + scalarField("value", dict, p.size()) + ); +} + + +adjointFarFieldTMVar1FvPatchScalarField:: +adjointFarFieldTMVar1FvPatchScalarField +( + const adjointFarFieldTMVar1FvPatchScalarField& tppsf, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(tppsf, iF), + adjointScalarBoundaryCondition(tppsf) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void adjointFarFieldTMVar1FvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + vectorField nf(patch().nf()); + + tmp tnuEff(boundaryContrPtr_->TMVariable1Diffusion()); + const scalarField& nuEff = tnuEff(); + const fvsPatchField& phip = boundaryContrPtr_->phib(); + const scalarField& magSf = patch().magSf(); + + // Patch-adjacent values + tmp intf(patchInternalField()); + + // Patch deltas + const scalarField& delta = patch().deltaCoeffs(); + + operator== + ( + pos(phip) + *( + (nuEff*delta*intf) + /(phip/magSf + nuEff*delta) + ) + ); + + fixedValueFvPatchScalarField::updateCoeffs(); +} + + +tmp> adjointFarFieldTMVar1FvPatchScalarField:: +valueInternalCoeffs +( + const tmp& +) const +{ + const fvsPatchField& phip = boundaryContrPtr_->phib(); + + // For fixedValue patches + return tmp>::New(neg(phip)*pTraits::one); +} + + +tmp> adjointFarFieldTMVar1FvPatchScalarField:: +valueBoundaryCoeffs +( + const tmp& +) const +{ + const fvsPatchField& phip = boundaryContrPtr_->phib(); + + // For zeroGradient patches + return tmp>::New(pos(phip)*(*this)); +} + + +void adjointFarFieldTMVar1FvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + writeEntry("value", os); + os.writeEntry("solverName", adjointSolverName_); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField +( + fvPatchScalarField, + adjointFarFieldTMVar1FvPatchScalarField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.H new file mode 100644 index 0000000000..291ea991b0 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar1/adjointFarFieldTMVar1FvPatchScalarField.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::adjointFarFieldTMVar1FvPatchScalarField + +Description + +SourceFiles + adjointFarFieldTMVar1FvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointFarFieldTMVar1FvPatchScalarField_H +#define adjointFarFieldTMVar1FvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" +#include "adjointBoundaryConditions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class adjointFarFieldTMVar1FvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class adjointFarFieldTMVar1FvPatchScalarField +: + public fixedValueFvPatchScalarField, + public adjointScalarBoundaryCondition +{ +public: + + //- Runtime type information + TypeName("adjointFarFieldTMVar1"); + + + // Constructors + + //- Construct from patch and internal field + adjointFarFieldTMVar1FvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + adjointFarFieldTMVar1FvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given adjointFarFieldTMVar1FvPatchScalarField + //- onto a new patch + adjointFarFieldTMVar1FvPatchScalarField + ( + const adjointFarFieldTMVar1FvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new adjointFarFieldTMVar1FvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + adjointFarFieldTMVar1FvPatchScalarField + ( + const adjointFarFieldTMVar1FvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new adjointFarFieldTMVar1FvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Return the matrix diagonal coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueInternalCoeffs + ( + const tmp& + ) const; + + //- Return the matrix source coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueBoundaryCoeffs + ( + const tmp& + ) const; + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.C new file mode 100644 index 0000000000..1bd84898be --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.C @@ -0,0 +1,178 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "adjointFarFieldTMVar2FvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" +#include "surfaceFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +adjointFarFieldTMVar2FvPatchScalarField:: +adjointFarFieldTMVar2FvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, word::null) +{} + + +adjointFarFieldTMVar2FvPatchScalarField:: +adjointFarFieldTMVar2FvPatchScalarField +( + const adjointFarFieldTMVar2FvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchScalarField(ptf, p, iF, mapper), + adjointScalarBoundaryCondition(p, iF, ptf.adjointSolverName_) +{} + + +adjointFarFieldTMVar2FvPatchScalarField:: +adjointFarFieldTMVar2FvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, dict.get("solverName")) +{ + fvPatchField::operator= + ( + scalarField("value", dict, p.size()) + ); +} + + +adjointFarFieldTMVar2FvPatchScalarField:: +adjointFarFieldTMVar2FvPatchScalarField +( + const adjointFarFieldTMVar2FvPatchScalarField& tppsf, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(tppsf, iF), + adjointScalarBoundaryCondition(tppsf) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void adjointFarFieldTMVar2FvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + vectorField nf(patch().nf()); + + tmp tnuEff(boundaryContrPtr_->TMVariable2Diffusion()); + const scalarField& nuEff = tnuEff(); + const fvsPatchField& phip = boundaryContrPtr_->phib(); + const scalarField& magSf = patch().magSf(); + + // Patch-adjacent values + tmp intf(patchInternalField()); + + // Patch deltas + const scalarField& delta = patch().deltaCoeffs(); + + operator== + ( + pos(phip) + *( + (nuEff*delta*intf) + /(phip/magSf + nuEff*delta) + ) + ); + + fixedValueFvPatchScalarField::updateCoeffs(); +} + + +tmp> adjointFarFieldTMVar2FvPatchScalarField:: +valueInternalCoeffs +( + const tmp& +) const +{ + const fvsPatchField& phip = boundaryContrPtr_->phib(); + + // For fixedValue patches + return tmp>::New(neg(phip)*pTraits::one); +} + + +tmp> adjointFarFieldTMVar2FvPatchScalarField:: +valueBoundaryCoeffs +( + const tmp& +) const +{ + const fvsPatchField& phip = boundaryContrPtr_->phib(); + + // For zeroGradient patches + return tmp>::New(pos(phip)*(*this)); +} + + +void adjointFarFieldTMVar2FvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + writeEntry("value", os); + os.writeEntry("solverName", adjointSolverName_); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField +( + fvPatchScalarField, + adjointFarFieldTMVar2FvPatchScalarField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.H new file mode 100644 index 0000000000..99d5674c26 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointFarFieldTMVar2/adjointFarFieldTMVar2FvPatchScalarField.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::adjointFarFieldTMVar2FvPatchScalarField + +Description + +SourceFiles + adjointFarFieldTMVar2FvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointFarFieldTMVar2FvPatchScalarField_H +#define adjointFarFieldTMVar2FvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" +#include "adjointBoundaryConditions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class adjointFarFieldTMVar2FvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class adjointFarFieldTMVar2FvPatchScalarField +: + public fixedValueFvPatchScalarField, + public adjointScalarBoundaryCondition +{ +public: + + //- Runtime type information + TypeName("adjointFarFieldTMVar2"); + + + // Constructors + + //- Construct from patch and internal field + adjointFarFieldTMVar2FvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + adjointFarFieldTMVar2FvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given adjointFarFieldTMVar2FvPatchScalarField + //- onto a new patch + adjointFarFieldTMVar2FvPatchScalarField + ( + const adjointFarFieldTMVar2FvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new adjointFarFieldTMVar2FvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + adjointFarFieldTMVar2FvPatchScalarField + ( + const adjointFarFieldTMVar2FvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new adjointFarFieldTMVar2FvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Return the matrix diagonal coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueInternalCoeffs + ( + const tmp& + ) const; + + //- Return the matrix source coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueBoundaryCoeffs + ( + const tmp& + ) const; + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.C new file mode 100644 index 0000000000..2ef57ef0c0 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.C @@ -0,0 +1,150 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2022 PCOpt/NTUA + Copyright (C) 2013-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "adjointOutletFluxFvPatchField.H" +#include "addToRunTimeSelectionTable.H" +#include "volFields.H" +#include "fvPatchFieldMapper.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::adjointOutletFluxFvPatchField::adjointOutletFluxFvPatchField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchField(p, iF) +{} + + +template +Foam::adjointOutletFluxFvPatchField::adjointOutletFluxFvPatchField +( + const adjointOutletFluxFvPatchField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchField(ptf, p, iF, mapper) +{} + + +template +Foam::adjointOutletFluxFvPatchField::adjointOutletFluxFvPatchField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchField(p, iF) +{ + fvPatchField::operator= + ( + Field("value", dict, p.size()) + ); +} + + +template +Foam::adjointOutletFluxFvPatchField::adjointOutletFluxFvPatchField +( + const adjointOutletFluxFvPatchField& tppsf, + const DimensionedField& iF +) +: + fixedValueFvPatchField(tppsf, iF) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::adjointOutletFluxFvPatchField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + this->operator==(Field(this->patch().size(), Zero)); + + fixedValueFvPatchField::updateCoeffs(); +} + + +template +Foam::tmp> +Foam::adjointOutletFluxFvPatchField::valueInternalCoeffs +( + const tmp& +) const +{ + return tmp>::New(this->size(), Zero); +} + + +template +Foam::tmp> +Foam::adjointOutletFluxFvPatchField::valueBoundaryCoeffs +( + const tmp& +) const +{ + return tmp>::New(this->size(), Zero); +} + + +template +Foam::tmp> +Foam::adjointOutletFluxFvPatchField::gradientBoundaryCoeffs() const +{ + return tmp>::New(this->size(), Zero); +} + + +template +Foam::tmp> +Foam::adjointOutletFluxFvPatchField::gradientInternalCoeffs() const +{ + return tmp>::New(this->size(), Zero); +} + + +template +void Foam::adjointOutletFluxFvPatchField::write(Ostream& os) const +{ + fvPatchField::write(os); + this->writeEntry("value", os); +} + + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.H new file mode 100644 index 0000000000..888ce70940 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchField.H @@ -0,0 +1,169 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2022 PCOpt/NTUA + Copyright (C) 2013-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::adjointOutletFluxFvPatchField + +Description + +SourceFiles + adjointOutletFluxFvPatchField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointOutletFluxFvPatchField_H +#define adjointOutletFluxFvPatchField_H + +#include "fvPatchFields.H" +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class adjointOutletFluxFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +template +class adjointOutletFluxFvPatchField +: + public fixedValueFvPatchField +{ +public: + + //- Runtime type information + TypeName("adjointOutletFlux"); + + + // Constructors + + //- Construct from patch and internal field + adjointOutletFluxFvPatchField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + adjointOutletFluxFvPatchField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given adjointOutletFluxFvPatchField + //- onto a new patch + adjointOutletFluxFvPatchField + ( + const adjointOutletFluxFvPatchField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp> clone() const + { + return tmp> + ( + new adjointOutletFluxFvPatchField(*this) + ); + } + + //- Construct as copy setting internal field reference + adjointOutletFluxFvPatchField + ( + const adjointOutletFluxFvPatchField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp> clone + ( + const DimensionedField& iF + ) const + { + return tmp> + ( + new adjointOutletFluxFvPatchField(*this, iF) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Return the matrix diagonal coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueInternalCoeffs + ( + const tmp& + ) const; + + //- Return the matrix source coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueBoundaryCoeffs + ( + const tmp& + + ) const; + + //- Return the matrix source coefficients corresponding to the + //- evaluation of the gradient of this patchField + virtual tmp> gradientBoundaryCoeffs() const; + + //- Return the matrix diagonal coefficients corresponding to the + //- evaluation of the gradient of this patchField + virtual tmp> gradientInternalCoeffs() const; + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "adjointOutletFluxFvPatchField.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.C new file mode 100644 index 0000000000..0c4d4242fa --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.C @@ -0,0 +1,46 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2022 PCOpt/NTUA + Copyright (C) 2013-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "adjointOutletFluxFvPatchFields.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +makePatchFields(adjointOutletFlux); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.H new file mode 100644 index 0000000000..40bf9694d0 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFields.H @@ -0,0 +1,53 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2022 PCOpt/NTUA + Copyright (C) 2013-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointOutletFluxFvPatchFields_H +#define adjointOutletFluxFvPatchFields_H + +#include "adjointOutletFluxFvPatchField.H" +#include "fieldTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeFieldTypedefs(adjointOutletFlux); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFieldsFwd.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFieldsFwd.H new file mode 100644 index 0000000000..aa754baa61 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletFlux/adjointOutletFluxFvPatchFieldsFwd.H @@ -0,0 +1,53 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2022 PCOpt/NTUA + Copyright (C) 2013-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointOutletFluxFvPatchFieldsFwd_H +#define adjointOutletFluxFvPatchFieldsFwd_H + +#include "fieldTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template class adjointOutletFluxFvPatchField; + +makePatchTypeFieldTypedefs(adjointOutletFlux); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.C new file mode 100644 index 0000000000..8a0aefde32 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.C @@ -0,0 +1,142 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2020 PCOpt/NTUA + Copyright (C) 2013-2020 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "adjointOutletKaFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +adjointOutletKaFvPatchScalarField::adjointOutletKaFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, word::null) +{} + + +adjointOutletKaFvPatchScalarField::adjointOutletKaFvPatchScalarField +( + const adjointOutletKaFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchScalarField(ptf, p, iF, mapper), + adjointScalarBoundaryCondition(p, iF, ptf.adjointSolverName_) +{} + + +adjointOutletKaFvPatchScalarField::adjointOutletKaFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, dict.get("solverName")) +{ + fvPatchField::operator= + ( + scalarField("value", dict, p.size()) + ); +} + + +adjointOutletKaFvPatchScalarField::adjointOutletKaFvPatchScalarField +( + const adjointOutletKaFvPatchScalarField& tppsf, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(tppsf, iF), + adjointScalarBoundaryCondition(tppsf) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void adjointOutletKaFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + vectorField nf(patch().nf()); + + const fvPatchField& Ub = boundaryContrPtr_->Ub(); + tmp tnuEff(boundaryContrPtr_->TMVariable1Diffusion()); + const scalarField& nuEff = tnuEff(); + + // Patch-adjacent ka + tmp tkaNei(patchInternalField()); + const scalarField& kaNei = tkaNei(); + + const scalarField& delta = patch().deltaCoeffs(); + + // Source from the objective + tmp source = boundaryContrPtr_->adjointTMVariable1Source(); + + operator== + ( + (nuEff*delta*kaNei - source) + /((Ub & nf) + nuEff*delta) + ); + + fixedValueFvPatchScalarField::updateCoeffs(); +} + + +void adjointOutletKaFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + writeEntry("value", os); + os.writeEntry("solverName", adjointSolverName_); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField(fvPatchScalarField, adjointOutletKaFvPatchScalarField); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.H new file mode 100644 index 0000000000..ae73acd7c3 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletKa/adjointOutletKaFvPatchScalarField.H @@ -0,0 +1,138 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2020 PCOpt/NTUA + Copyright (C) 2013-2020 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::adjointOutletKaFvPatchScalarField + +Description + +SourceFiles + adjointOutletKaFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointOutletKaFvPatchScalarField_H +#define adjointOutletKaFvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" +#include "adjointBoundaryConditions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class adjointOutletKaFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class adjointOutletKaFvPatchScalarField +: + public fixedValueFvPatchScalarField, + public adjointScalarBoundaryCondition +{ +public: + + //- Runtime type information + TypeName("adjointOutletKa"); + + + // Constructors + + //- Construct from patch and internal field + adjointOutletKaFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + adjointOutletKaFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given adjointOutletKaFvPatchScalarField + //- onto a new patch + adjointOutletKaFvPatchScalarField + ( + const adjointOutletKaFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new adjointOutletKaFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + adjointOutletKaFvPatchScalarField + ( + const adjointOutletKaFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new adjointOutletKaFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.C new file mode 100644 index 0000000000..e920a5d26b --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.C @@ -0,0 +1,142 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2020 PCOpt/NTUA + Copyright (C) 2013-2020 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "adjointOutletWaFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +adjointOutletWaFvPatchScalarField::adjointOutletWaFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, word::null) +{} + + +adjointOutletWaFvPatchScalarField::adjointOutletWaFvPatchScalarField +( + const adjointOutletWaFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchScalarField(ptf, p, iF, mapper), + adjointScalarBoundaryCondition(p, iF, ptf.adjointSolverName_) +{} + + +adjointOutletWaFvPatchScalarField::adjointOutletWaFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchScalarField(p, iF), + adjointScalarBoundaryCondition(p, iF, dict.get("solverName")) +{ + fvPatchField::operator= + ( + scalarField("value", dict, p.size()) + ); +} + + +adjointOutletWaFvPatchScalarField::adjointOutletWaFvPatchScalarField +( + const adjointOutletWaFvPatchScalarField& tppsf, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(tppsf, iF), + adjointScalarBoundaryCondition(tppsf) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void adjointOutletWaFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + vectorField nf(patch().nf()); + + const fvPatchField& Ub = boundaryContrPtr_->Ub(); + tmp tnuEff(boundaryContrPtr_->TMVariable2Diffusion()); + const scalarField& nuEff = tnuEff(); + + // Patch-adjacent wa + tmp twaNei = patchInternalField(); + const scalarField& waNei = twaNei(); + + const scalarField& delta = patch().deltaCoeffs(); + + // Source from the objective + tmp source = boundaryContrPtr_->adjointTMVariable2Source(); + + operator== + ( + (nuEff*delta*waNei - source) + /((Ub & nf) + nuEff*delta) + ); + + fixedValueFvPatchScalarField::updateCoeffs(); +} + + +void adjointOutletWaFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + writeEntry("value", os); + os.writeEntry("solverName", adjointSolverName_); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField(fvPatchScalarField, adjointOutletWaFvPatchScalarField); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.H new file mode 100644 index 0000000000..dd04e830f3 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/adjointOutletWa/adjointOutletWaFvPatchScalarField.H @@ -0,0 +1,138 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2007-2020 PCOpt/NTUA + Copyright (C) 2013-2020 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::adjointOutletWaFvPatchScalarField + +Description + +SourceFiles + adjointOutletWaFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adjointOutletWaFvPatchScalarField_H +#define adjointOutletWaFvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" +#include "adjointBoundaryConditions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class adjointOutletWaFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class adjointOutletWaFvPatchScalarField +: + public fixedValueFvPatchScalarField, + public adjointScalarBoundaryCondition +{ +public: + + //- Runtime type information + TypeName("adjointOutletWa"); + + + // Constructors + + //- Construct from patch and internal field + adjointOutletWaFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + adjointOutletWaFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given adjointOutletWaFvPatchScalarField + //- onto a new patch + adjointOutletWaFvPatchScalarField + ( + const adjointOutletWaFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new adjointOutletWaFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + adjointOutletWaFvPatchScalarField + ( + const adjointOutletWaFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new adjointOutletWaFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.C new file mode 100644 index 0000000000..0057310ad0 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.C @@ -0,0 +1,176 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "kaqRWallFunctionFvPatchScalarField.H" +#include "fvPatchFieldMapper.H" +#include "nutkWallFunctionFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::kaqRWallFunctionFvPatchScalarField::kaqRWallFunctionFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + kqRWallFunctionFvPatchField(p, iF), + adjointScalarBoundaryCondition(p, iF, word::null) +{} + + +Foam::kaqRWallFunctionFvPatchScalarField::kaqRWallFunctionFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + kqRWallFunctionFvPatchField(p, iF, dict), + adjointScalarBoundaryCondition(p, iF, dict.get("solverName")) +{} + + +Foam::kaqRWallFunctionFvPatchScalarField::kaqRWallFunctionFvPatchScalarField +( + const kaqRWallFunctionFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + kqRWallFunctionFvPatchField(ptf, p, iF, mapper), + adjointScalarBoundaryCondition(p, iF, ptf.adjointSolverName_) +{} + + +Foam::kaqRWallFunctionFvPatchScalarField::kaqRWallFunctionFvPatchScalarField +( + const kaqRWallFunctionFvPatchScalarField& tkqrwfpf +) +: + kqRWallFunctionFvPatchField(tkqrwfpf), + adjointScalarBoundaryCondition(tkqrwfpf) +{} + + +Foam::kaqRWallFunctionFvPatchScalarField::kaqRWallFunctionFvPatchScalarField +( + const kaqRWallFunctionFvPatchScalarField& tkqrwfpf, + const DimensionedField& iF +) +: + kqRWallFunctionFvPatchField(tkqrwfpf, iF), + adjointScalarBoundaryCondition(tkqrwfpf) +{} + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::kaqRWallFunctionFvPatchScalarField::manipulateMatrix +( + fvMatrix& matrix +) +{ + scalarField& source = matrix.source(); + tmp tnutWall(boundaryContrPtr_->turbulentDiffusivity()); + fvPatchScalarField& nutWall = tnutWall.constCast(); + if (isA(nutWall)) + { + const label patchi(patch().index()); + const scalarField& magSf = patch().magSf(); + const turbulenceModel& turbModel = db().lookupObject + ( + IOobject::groupName + ( + turbulenceModel::propertiesName, + internalField().group() + ) + ); + const scalarField& y = turbModel.y()[patchi]; + const tmp tnuw = turbModel.nu(patchi); + const scalarField& nuw = tnuw(); + const nutWallFunctionFvPatchScalarField& nutWF = + refCast(nutWall); + const wallFunctionCoefficients& wallCoeffs = nutWF.wallCoeffs(); + const scalar Cmu = wallCoeffs.Cmu(); + const scalar kappa = wallCoeffs.kappa(); + const scalar E = wallCoeffs.E(); + const scalar yPlusLam = wallCoeffs.yPlusLam(); + + const scalar Cmu25 = pow025(Cmu); + + const labelList& faceCells = patch().faceCells(); + const fvPatchVectorField& Uw = boundaryContrPtr_->Ub(); + const scalarField magGradUw(mag(Uw.snGrad())); + + tmp tdJdnut = boundaryContrPtr_->dJdnut(); + const scalarField& dJdnut = tdJdnut(); + + const tmp tk = turbModel.k(); + const volScalarField& k = tk(); + forAll(dJdnut, facei) + { + const label celli = faceCells[facei]; + + const scalar sqrtkCell(sqrt(k[celli])); + const scalar yPlus = Cmu25*y[facei]*sqrtkCell/nuw[facei]; + if (yPlusLam < yPlus) + { + const scalar logEyPlus = log(E*yPlus); + const scalar dnut_dyPlus = + nuw[facei]*kappa*(logEyPlus - 1)/sqr(logEyPlus); + const scalar dyPlus_dk = + Cmu25*y[facei]/(2*nuw[facei]*sqrtkCell); + const scalar dnut_dk = dnut_dyPlus*dyPlus_dk; + source[celli] -= dJdnut[facei]*dnut_dk*magSf[facei]; + } + } + } +} + + +void Foam::kaqRWallFunctionFvPatchScalarField::write(Ostream& os) const +{ + kqRWallFunctionFvPatchField::write(os); + os.writeEntry("solverName", adjointSolverName_); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchScalarField, + kaqRWallFunctionFvPatchScalarField + ); +} + + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.H new file mode 100644 index 0000000000..7f93602416 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/kaqRWallFunction/kaqRWallFunctionFvPatchScalarField.H @@ -0,0 +1,144 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::kaqRWallFunctionFvPatchScalarField + +Description + +SourceFiles + kaqRWallFunctionFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef kaqRWallFunctionFvPatchScalarField_H +#define kaqRWallFunctionFvPatchScalarField_H + +#include "kqRWallFunctionFvPatchField.H" +#include "adjointBoundaryConditions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class kaqRWallFunctionFvPatchScalarField Declaration +\*---------------------------------------------------------------------------*/ + +class kaqRWallFunctionFvPatchScalarField +: + public kqRWallFunctionFvPatchField, + public adjointScalarBoundaryCondition +{ +public: + + //- Runtime type information + TypeName("kaqRWallFunction"); + + + // Constructors + + //- Construct from patch and internal field + kaqRWallFunctionFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + kaqRWallFunctionFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given kaqRWallFunctionFvPatchScalarField + //- onto a new patch + kaqRWallFunctionFvPatchScalarField + ( + const kaqRWallFunctionFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + kaqRWallFunctionFvPatchScalarField + ( + const kaqRWallFunctionFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new kaqRWallFunctionFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + kaqRWallFunctionFvPatchScalarField + ( + const kaqRWallFunctionFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new kaqRWallFunctionFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + //- Add source terms to the rhs of the first cell centre + virtual void manipulateMatrix(fvMatrix& matrix); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.C new file mode 100644 index 0000000000..a3dc8d51c8 --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.C @@ -0,0 +1,318 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "waWallFunctionFvPatchScalarField.H" +#include "wallFvPatch.H" +#include "omegaWallFunctionFvPatchScalarField.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void waWallFunctionFvPatchScalarField::checkType() +{ + if (!isA(patch())) + { + FatalErrorInFunction + << "Invalid wall function specification" << nl + << " Patch type for patch " << patch().name() + << " must be wall" << nl + << " Current patch type is " << patch().type() << nl << endl + << abort(FatalError); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +waWallFunctionFvPatchScalarField::waWallFunctionFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchField(p, iF), + adjointScalarBoundaryCondition(p, iF, "wa") +{ + checkType(); +} + + +waWallFunctionFvPatchScalarField::waWallFunctionFvPatchScalarField +( + const waWallFunctionFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchField(ptf, p, iF, mapper), + adjointScalarBoundaryCondition(p, iF, ptf.adjointSolverName_) +{ + checkType(); +} + + +waWallFunctionFvPatchScalarField::waWallFunctionFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchField(p, iF, dict), + adjointScalarBoundaryCondition(p, iF, dict.get("solverName")) +{ + checkType(); +} + + +waWallFunctionFvPatchScalarField::waWallFunctionFvPatchScalarField +( + const waWallFunctionFvPatchScalarField& ewfpsf +) +: + fixedValueFvPatchField(ewfpsf), + adjointScalarBoundaryCondition(ewfpsf) +{ + checkType(); +} + + +waWallFunctionFvPatchScalarField::waWallFunctionFvPatchScalarField +( + const waWallFunctionFvPatchScalarField& ewfpsf, + const DimensionedField& iF +) +: + fixedValueFvPatchField(ewfpsf, iF), + adjointScalarBoundaryCondition(ewfpsf) +{ + checkType(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void waWallFunctionFvPatchScalarField::manipulateMatrix +( + fvMatrix& matrix +) +{ + scalarField& Diag = matrix.diag(); + scalarField& lower = matrix.lower(); + scalarField& upper = matrix.upper(); + FieldField& internalCoeffs = matrix.internalCoeffs(); + FieldField& boundaryCoeffs = matrix.boundaryCoeffs(); + const fvMesh& mesh = patch().boundaryMesh().mesh(); + const labelList& faceCells = patch().faceCells(); + + // Add diag term from the omega expression next to the wall + for (const label celli : faceCells) + { + Diag[celli] = 1; + } + + // We want something similar to setValues, but slightly modified. + // The solution of the boundary cell should understand contributions from + // the second cells off the wall but they should see the + // solution of the boundary cell as zero. + // Contributions from neighbouring cells with an omegaWallFunction boundary + // condition should also be zero + + const cellList& cells = mesh.cells(); + const labelUList& own = mesh.owner(); + + /* + const labelUList& nei = mesh.neighbour(); + const turbulenceModel& turbModel = db().lookupObject + ( + IOobject::groupName + ( + turbulenceModel::propertiesName, + internalField().group() + ) + ); + + const tmp tomega = turbModel.omega(); + const volScalarField& omega = tomega(); + typedef omegaWallFunctionFvPatchScalarField omegaWF; + */ + + forAll(faceCells, i) + { + const label celli = faceCells[i]; + const cell& c = cells[celli]; + + forAll(c, j) + { + const label facei = c[j]; + + if (mesh.isInternalFace(facei)) + { + // Neighbouring cells should get no contribution from + // ourselves in all cases + //label cellNei(-1); + if (celli == own[facei]) + { + //cellNei = nei[facei]; + lower[facei] = 0.0; + } + else + { + //cellNei = own[facei]; + upper[facei] = 0.0; + } + // Additionally, if the neighbouring cell is also a boundary + // one with omegaWallFunction in one of its faces, + // contributions between the two cells should be canceled out + // as well. + // Already covered by the above + /* + bool neiHasOmegaWFface(false); + const cell& neiCell = cells[cellNei]; + forAll(neiCell, fNei) + { + const label faceNei = neiCell[fNei]; + + const label patchNei = + mesh.boundaryMesh().whichPatch(faceNei); + if (patchNei != -1) + { + const fvPatchField& omegaNei = + omega.boundaryField()[patchNei]; + if (isA(omegaNei)) + { + neiHasOmegaWFface = true; + break; + } + } + } + if (neiHasOmegaWFface) + { + if (celli == own[facei]) + { + upper[facei] = 0.0; + } + else + { + lower[facei] = 0.0; + } + } + */ + } + // Contributions from boundaries should have already been removed + // using value*Coeffs and boundary*Coeffs + // Just to be safe + else + { + const label patchi = mesh.boundaryMesh().whichPatch(facei); + if (internalCoeffs[patchi].size()) + { + label patchFacei = + mesh.boundaryMesh()[patchi].whichFace(facei); + internalCoeffs[patchi][patchFacei] = Zero; + boundaryCoeffs[patchi][patchFacei] = Zero; + } + } + } + } + + fvPatchField::manipulateMatrix(matrix); +} + + +void waWallFunctionFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + + operator == (scalarField(patch().size(), Zero)); + fixedValueFvPatchField::updateCoeffs(); +} + + +tmp> waWallFunctionFvPatchScalarField::valueInternalCoeffs +( + const tmp& +) const +{ + return tmp>::New(this->size(), Zero); +} + + +tmp> waWallFunctionFvPatchScalarField::valueBoundaryCoeffs +( + const tmp& +) const +{ + return tmp>::New(this->size(), Zero); +} + + +tmp> +waWallFunctionFvPatchScalarField::gradientBoundaryCoeffs() const +{ + return tmp>::New(this->size(), Zero); +} + + +tmp> +waWallFunctionFvPatchScalarField::gradientInternalCoeffs() const +{ + return tmp>::New(this->size(), Zero); +} + + +void waWallFunctionFvPatchScalarField::write(Ostream& os) const +{ + fixedValueFvPatchField::write(os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField +( + fvPatchScalarField, + waWallFunctionFvPatchScalarField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.H new file mode 100644 index 0000000000..c1e95db2aa --- /dev/null +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/derivedFvPatchFields/waWallFunction/waWallFunctionFvPatchScalarField.H @@ -0,0 +1,184 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2014-2022 PCOpt/NTUA + Copyright (C) 2014-2022 FOSS GP +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + + +Class + Foam::waWallFunctionFvPatchScalarField + +Description + +SourceFiles + waWallFunctionFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef waWallFunctionFvPatchScalarField_H +#define waWallFunctionFvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" +#include "adjointBoundaryConditions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class waWallFunctionFvPatchScalarField Declaration +\*---------------------------------------------------------------------------*/ + +class waWallFunctionFvPatchScalarField +: + public fixedValueFvPatchField, + public adjointScalarBoundaryCondition +{ + + // Private member functions + + //- Check the type of the patch + void checkType(); + + +public: + + //- Runtime type information + TypeName("waWallFunction"); + + + // Constructors + + //- Construct from patch and internal field + waWallFunctionFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + waWallFunctionFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given + // waWallFunctionFvPatchScalarField + // onto a new patch + waWallFunctionFvPatchScalarField + ( + const waWallFunctionFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + waWallFunctionFvPatchScalarField + ( + const waWallFunctionFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new waWallFunctionFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + waWallFunctionFvPatchScalarField + ( + const waWallFunctionFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new waWallFunctionFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Evaluation functions + + //- Manipulate waEqn at the cells next to the wall + virtual void manipulateMatrix(fvMatrix& matrix); + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Return the matrix diagonal coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueInternalCoeffs + ( + const tmp& + ) const; + + //- Return the matrix source coefficients corresponding to the + //- evaluation of the value of this patchField with given weights + virtual tmp> valueBoundaryCoeffs + ( + const tmp& + + ) const; + + //- Return the matrix source coefficients corresponding to the + //- evaluation of the gradient of this patchField + virtual tmp> gradientBoundaryCoeffs() const; + + //- Return the matrix diagonal coefficients corresponding to the + //- evaluation of the gradient of this patchField + virtual tmp> gradientInternalCoeffs() const; + + + + // I-O + + //- Write + void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H index 09800b9cb9..9d729bdd59 100644 --- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H @@ -90,8 +90,8 @@ protected: // Protected Member functions - void allocateInitValues(); - void allocateMeanFields(); + virtual void allocateInitValues(); + virtual void allocateMeanFields(); refPtr cloneRefPtr(const refPtr& obj) const; @@ -213,6 +213,13 @@ public: const singlePhaseTransportModel& laminarTransport ) const; + //- Return the turbulence production term + virtual tmp G() + { + NotImplemented; + return nullptr; + } + //- Restore turbulent fields to their initial values void restoreInitValues(); @@ -220,7 +227,7 @@ public: void resetMeanFields(); //- Compute mean fields on the fly - void computeMeanFields(); + virtual void computeMeanFields(); //- Return stress tensor based on the mean flow variables tmp devReff diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C index da8eed8cc1..1d393e3113 100644 --- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C @@ -28,6 +28,7 @@ License \*---------------------------------------------------------------------------*/ #include "kOmegaSST.H" +#include "wallDist.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -44,6 +45,34 @@ namespace RASVariables defineTypeNameAndDebug(kOmegaSST, 0); addToRunTimeSelectionTable(RASModelVariables, kOmegaSST, dictionary); + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void kOmegaSST::allocateMeanFields() +{ + RASModelVariables::allocateMeanFields(); + if (solverControl_.average()) + { + GMean_.reset + ( + new volScalarField::Internal + ( + IOobject + ( + "GMean", + mesh_.time().timeName(), + mesh_, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh_, + dimensionedScalar(dimArea/pow3(dimTime), Zero) + ) + ); + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // kOmegaSST::kOmegaSST @@ -61,13 +90,76 @@ kOmegaSST::kOmegaSST TMVar2Ptr_.ref(mesh_.lookupObjectRef(TMVar2BaseName_)); nutPtr_.ref(mesh_.lookupObjectRef(nutBaseName_)); + distPtr_.ref(const_cast(wallDist::New(mesh_).y())); + allocateInitValues(); allocateMeanFields(); } +tmp kOmegaSST::computeG() +{ + const turbulenceModel& turbModel = mesh_.lookupObject + ( + IOobject::groupName + ( + turbulenceModel::propertiesName, + TMVar2().internalField().group() + ) + ); + // Recompute G and modify values next to the walls + // Ideally, grad(U) should be cached to avoid the overhead + const volVectorField& U = turbModel.U(); + tmp tgradU = fvc::grad(U); + volScalarField::Internal GbyNu0 + ( + this->type() + ":GbyNu", + (tgradU() && dev(twoSymm(tgradU()))) + ); + auto tG = + tmp::New + ( + turbModel.GName(), + nutRefInst()*GbyNu0 + ); + // Use correctBoundaryConditions instead of updateCoeffs to avoid + // messing with updateCoeffs in the next iteration of omegaEqn + TMVar2Inst().correctBoundaryConditions(); + + return tG; +} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +tmp kOmegaSST::G() +{ + if (solverControl_.useAveragedFields()) + { + DebugInfo + << "Using GMean" << endl; + return tmp(GMean_()); + } + DebugInfo + << "Using instantaneous G" << endl; + return computeG(); +} + + +void kOmegaSST::computeMeanFields() +{ + RASModelVariables::computeMeanFields(); + if (solverControl_.doAverageIter()) + { + const label iAverageIter = solverControl_.averageIter(); + scalar avIter(iAverageIter); + scalar oneOverItP1 = 1./(avIter + 1); + scalar mult = avIter*oneOverItP1; + GMean_() = GMean_()*mult + computeG()*oneOverItP1; + } +} + + void kOmegaSST::correctBoundaryConditions ( const incompressible::turbulenceModel& turbulence diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.H index 1672a76bc6..ce62241b44 100644 --- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.H +++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.H @@ -59,6 +59,21 @@ class kOmegaSST public RASModelVariables { +protected: + + // Protected data + + //- Average of the production term + // Avaraged separetely due the bi-zonal treatment next to the wall + autoPtr GMean_; + + + // Protected Member functions + + virtual void allocateMeanFields(); + tmp computeG(); + + public: //- Runtime type information @@ -81,6 +96,12 @@ public: // Member Functions + //- Return the turbulence production term + virtual tmp G(); + + //- Compute mean fields on the fly + virtual void computeMeanFields(); + //- Correct boundary conditions of turbulent fields virtual void correctBoundaryConditions (