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 (