DAResidualTurboFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAResidualTurboFoam.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAResidualTurboFoam, 0);
16 addToRunTimeSelectionTable(DAResidual, DAResidualTurboFoam, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
20  const word modelType,
21  const fvMesh& mesh,
22  const DAOption& daOption,
23  const DAModel& daModel,
24  const DAIndex& daIndex)
25  : DAResidual(modelType, mesh, daOption, daModel, daIndex),
26  // initialize and register state variables and their residuals, we use macros defined in macroFunctions.H
27  setResidualClassMemberVector(U, dimensionSet(1, -2, -2, 0, 0, 0, 0)),
28  setResidualClassMemberScalar(p, dimensionSet(1, -3, -1, 0, 0, 0, 0)),
29  setResidualClassMemberScalar(T, dimensionSet(1, -1, -3, 0, 0, 0, 0)),
31  thermo_(const_cast<fluidThermo&>(
32  mesh_.thisDb().lookupObject<fluidThermo>("thermophysicalProperties"))),
33  URel_(const_cast<volVectorField&>(
34  mesh_.thisDb().lookupObject<volVectorField>("URel"))),
35  he_(thermo_.he()),
36  rho_(const_cast<volScalarField&>(
37  mesh_.thisDb().lookupObject<volScalarField>("rho"))),
38  alphat_(const_cast<volScalarField&>(
39  mesh_.thisDb().lookupObject<volScalarField>("alphat"))),
40  psi_(const_cast<volScalarField&>(
41  mesh_.thisDb().lookupObject<volScalarField>("thermo:psi"))),
42  daTurb_(const_cast<DATurbulenceModel&>(daModel.getDATurbulenceModel())),
43  // create simpleControl
44  simple_(const_cast<fvMesh&>(mesh)),
45  pressureControl_(p_, rho_, simple_.dict()),
46  MRF_(const_cast<IOMRFZoneListDF&>(
47  mesh_.thisDb().lookupObject<IOMRFZoneListDF>("MRFProperties")))
48 {
49 
50  // get molWeight and Cp from thermophysicalProperties
51  const IOdictionary& thermoDict = mesh.thisDb().lookupObject<IOdictionary>("thermophysicalProperties");
52  dictionary mixSubDict = thermoDict.subDict("mixture");
53  dictionary specieSubDict = mixSubDict.subDict("specie");
54  molWeight_ = specieSubDict.getScalar("molWeight");
55  dictionary thermodynamicsSubDict = mixSubDict.subDict("thermodynamics");
56  Cp_ = thermodynamicsSubDict.getScalar("Cp");
57 
58  if (daOption_.getOption<label>("debug"))
59  {
60  Info << "molWeight " << molWeight_ << endl;
61  Info << "Cp " << Cp_ << endl;
62  }
63 }
64 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
65 
67 {
68  /*
69  Description:
70  Clear all members to avoid memory leak because we will initalize
71  multiple objects of DAResidual. Here we need to delete all members
72  in the parent and child classes
73  */
74  URes_.clear();
75  pRes_.clear();
76  TRes_.clear();
77  phiRes_.clear();
78 }
79 
80 void DAResidualTurboFoam::calcResiduals(const dictionary& options)
81 {
82  /*
83  Description:
84  This is the function to compute residuals.
85 
86  Input:
87  options.isPC: 1 means computing residuals for preconditioner matrix.
88  This essentially use the first order scheme for div(phi,U), div(phi,e)
89 
90  p_, T_, U_, phi_, etc: State variables in OpenFOAM
91 
92  Output:
93  URes_, pRes_, TRes_, phiRes_: residual field variables
94  */
95 
96  // We dont support MRF and fvOptions so all the related lines are commented
97  // out for now
98 
99  label isPC = options.getLabel("isPC");
100 
101  word divUScheme = "div(phi,U)";
102  word divHEScheme = "div(phi,e)";
103  word divPhidPScheme = "div(phid,p)";
104 
105  if (he_.name() == "h")
106  {
107  divHEScheme = "div(phi,h)";
108  }
109 
110  if (isPC)
111  {
112  divUScheme = "div(pc)";
113  divHEScheme = "div(pc)";
114  divPhidPScheme = "div(pc)";
115  }
116 
117  // ******** U Residuals **********
118  // copied and modified from UEqn.H
119 
120  tmp<fvVectorMatrix> tUEqn(
121  fvm::div(phi_, U_, divUScheme)
122  + MRF_.DDt(rho_, U_)
124  fvVectorMatrix& UEqn = tUEqn.ref();
125 
126  UEqn.relax();
127 
128  URes_ = (UEqn & U_) + fvc::grad(p_);
129  normalizeResiduals(URes);
130 
131  // ******** e Residuals **********
132  // copied and modified from EEqn.H
133  volSymmTensorField Teff = -daTurb_.devRhoReff();
134  volScalarField alphaEff("alphaEff", thermo_.alphaEff(alphat_));
135 
136  URel_ == U_;
138 
139  fvScalarMatrix EEqn(
140  fvm::div(phi_, he_, divHEScheme)
141  + (he_.name() == "e"
142  ? fvc::div(phi_, volScalarField("Ekp", 0.5 * magSqr(U_) + p_ / rho_))
143  : fvc::div(phi_, volScalarField("K", 0.5 * magSqr(U_))) - fvc::div(Teff.T() & U_) + fvc::div(p_ * (U_ - URel_)))
144  - fvm::Sp(fvc::div(phi_), he_)
145  - fvm::laplacian(alphaEff, he_));
146 
147  EEqn.relax();
148 
149  TRes_ = EEqn & he_;
150  normalizeResiduals(TRes);
151 
152  // ******** p and phi Residuals **********
153  // copied and modified from pEqn.H
154  volScalarField AU(UEqn.A());
155  volScalarField AtU(AU - UEqn.H1());
156  volVectorField HbyA("HbyA", U_);
157  HbyA = UEqn.H() / AU;
158 
159  volScalarField rAU(1.0 / UEqn.A());
160  tUEqn.clear();
161 
162  if (simple_.transonic())
163  {
164  surfaceScalarField phid(
165  "phid",
166  fvc::interpolate(psi_) * (fvc::interpolate(HbyA) & mesh_.Sf()));
167 
168  MRF_.makeRelative(fvc::interpolate(psi_), phid);
169 
170  fvScalarMatrix pEqn(
171  fvm::div(phid, p_, divPhidPScheme)
172  - fvm::laplacian(rho_ * rAU, p_));
173 
174  // for PC we do not include the div(phid, p) term, this improves the convergence
175  if (isPC && daOption_.getOption<label>("transonicPCOption") == 1)
176  {
177  pEqn -= fvm::div(phid, p_, divPhidPScheme);
178  }
179 
180  // Relax the pressure equation to maintain diagonal dominance
181  pEqn.relax();
182 
183  pEqn.setReference(pressureControl_.refCell(), pressureControl_.refValue());
184 
185  pRes_ = pEqn & p_;
186  normalizeResiduals(pRes);
187 
188  // ******** phi Residuals **********
189  // copied and modified from pEqn.H
190  if (isPC && daOption_.getOption<label>("transonicPCOption") == 2)
191  {
192  // transonic PC option 2, we ignore all the off-diagonal
193  // terms for the phiRes
194  phiRes_ = phi_;
195  }
196  else
197  {
198  phiRes_ == pEqn.flux() - phi_;
199  }
200 
201  // need to normalize phiRes
202  normalizePhiResiduals(phiRes);
203  }
204  else
205  {
206 
207  surfaceScalarField phiHbyA(
208  "phiHbyA",
209  fvc::interpolate(rho_ * HbyA) & mesh_.Sf());
210 
211  MRF_.makeRelative(fvc::interpolate(rho_), phiHbyA);
212 
213  adjustPhi(phiHbyA, U_, p_);
214  phiHbyA += fvc::interpolate(rho_ / AtU - rho_ / AU) * fvc::snGrad(p_) * mesh_.magSf();
215 
216  fvScalarMatrix pEqn(
217  fvc::div(phiHbyA)
218  - fvm::laplacian(rho_ / AtU, p_));
219 
220  pEqn.setReference(pressureControl_.refCell(), pressureControl_.refValue());
221 
222  pRes_ = pEqn & p_;
223  normalizeResiduals(pRes);
224 
225  // ******** phi Residuals **********
226  // copied and modified from pEqn.H
227  // TODO: the phiRes is not zero, need to fix
228  phiRes_ == phiHbyA + pEqn.flux() - phi_;
229  normalizePhiResiduals(phiRes);
230  }
231 }
232 
234 {
235  /*
236  Description:
237  Update the intermediate variables that depend on the state variables
238 
239  ********************** NOTE *****************
240  we assume hePsiThermo, pureMixture, perfectGas, hConst, and const transport
241  TODO: need to do this using built-in openfoam functions.
242 
243  we need to:
244  1, update psi based on T, psi=1/(R*T)
245  2, update rho based on p and psi, rho=psi*p
246  3, update E based on T, p and rho, E=Cp*T-p/rho
247  4, update velocity boundary based on MRF
248  */
249 
250  // 8314.4700665 gas constant in OpenFOAM
251  // src/OpenFOAM/global/constants/thermodynamic/thermodynamicConstants.H
252  scalar RR = Foam::constant::thermodynamic::RR;
253 
254  // R = RR/molWeight
255  // Foam::specie::R() function in src/thermophysicalModels/specie/specie/specieI.H
256  dimensionedScalar R(
257  "R1",
258  dimensionSet(0, 2, -2, -1, 0, 0, 0),
259  RR / molWeight_);
260 
261  // psi = 1/T/R
262  // see src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGasI.H
263  psi_ = 1.0 / T_ / R;
264 
265  // rho = psi*p
266  // see src/thermophysicalModels/basic/psiThermo/psiThermo.C
267  rho_ = psi_ * p_;
268 
269  // **************** NOTE ****************
270  // need to relax rho to be consistent with the primal solver
271  // However, the rho.relax() will mess up perturbation
272  // That being said, we comment out the rho.relax() call to
273  // get the correct perturbed rho; however, the E residual will
274  // be a bit off compared with the ERes at the converged state
275  // from the primal solver. TODO. Need to figure out how to improve this
276  // **************** NOTE ****************
277  // rho_.relax();
278 
279  dimensionedScalar Cp(
280  "Cp1",
281  dimensionSet(0, 2, -2, -1, 0, 0, 0),
282  Cp_);
283 
284  // Hs = Cp*T
285  // see Hs() in src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
286  // here the H departure EquationOfState::H(p, T) will be zero for perfectGas
287  // Es = Hs - p/rho = Hs - T * R;
288  // see Es() in src/thermophysicalModels/specie/thermo/thermo/thermoI.H
289  // **************** NOTE ****************
290  // See the comment from the rho.relax() call, if we write he_=Cp*T-p/rho, the
291  // accuracy of he_ may be impact by the inaccurate rho. So here we want to
292  // rewrite he_ as he_ = Cp * T_ - T_ * R instead, such that we dont include rho
293  // **************** NOTE ****************
294  if (he_.name() == "e")
295  {
296  he_ = Cp * T_ - T_ * R;
297  }
298  else
299  {
300  he_ = Cp * T_;
301  }
302  he_.correctBoundaryConditions();
303 
304  // NOTE: alphat is updated in the correctNut function in DATurbulenceModel child classes
305 
307 }
308 
310 {
311  /*
312  Description:
313  Update the boundary condition for all the states in the selected solver
314  */
315 
316  U_.correctBoundaryConditions();
317  p_.correctBoundaryConditions();
318  T_.correctBoundaryConditions();
319 }
320 
321 } // End namespace Foam
322 
323 // ************************************************************************* //
phid
surfaceScalarField phid("phid",(fvc::interpolate(psi)/fvc::interpolate(rho)) *phiHbyA)
Foam::MRFZoneListDF::DDt
tmp< volVectorField > DDt(const volVectorField &U) const
Definition: MRFZoneListDF.C:120
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Foam::DAResidualTurboFoam::pRes_
volScalarField pRes_
Definition: DAResidualTurboFoam.H:41
Foam::DAResidualTurboFoam::URes_
volVectorField URes_
Definition: DAResidualTurboFoam.H:38
Foam::DATurbulenceModel::divDevRhoReff
tmp< fvVectorMatrix > divDevRhoReff(volVectorField &U)
divDev terms
Definition: DATurbulenceModel.C:389
Foam::DAResidual::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAResidual.H:49
Foam::DAResidualTurboFoam::clear
virtual void clear()
clear the members
Definition: DAResidualTurboFoam.C:66
Foam::DAResidualTurboFoam::TRes_
volScalarField TRes_
Definition: DAResidualTurboFoam.H:44
Foam::DAResidualTurboFoam::p_
volScalarField & p_
Definition: DAResidualTurboFoam.H:40
Foam::DAOption
Definition: DAOption.H:29
Foam::DAResidualTurboFoam::calcResiduals
virtual void calcResiduals(const dictionary &options)
compute residual
Definition: DAResidualTurboFoam.C:80
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
UEqn
fvVectorMatrix & UEqn
Definition: UEqnPimpleDyM.H:15
Foam::DAResidualTurboFoam::DAResidualTurboFoam
DAResidualTurboFoam(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAResidualTurboFoam.C:19
he
volScalarField & he
Definition: EEqnRhoPimple.H:5
normalizeResiduals
#define normalizeResiduals(resName)
Definition: DAMacroFunctions.H:28
Foam::DAResidualTurboFoam::simple_
simpleControl simple_
simpleControl object which will be initialized in this class
Definition: DAResidualTurboFoam.H:68
Foam::DAResidualTurboFoam::phi_
surfaceScalarField & phi_
Definition: DAResidualTurboFoam.H:46
Foam::DAOption::getOption
classType getOption(const word key) const
get an option from subDict and key
Definition: DAOption.H:92
Foam::DAResidualTurboFoam::daTurb_
DATurbulenceModel & daTurb_
DATurbulenceModel object.
Definition: DAResidualTurboFoam.H:65
Foam::DAResidual::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAResidual.H:52
HbyA
volVectorField & HbyA
Definition: pEqnPimpleDyM.H:17
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
AU
volScalarField AU(UEqn.A())
setResidualClassMemberPhi
#define setResidualClassMemberPhi(stateName)
Definition: DAMacroFunctions.H:83
phiHbyA
phiHbyA
Definition: pEqnRhoSimpleC.H:36
Foam::DAIndex
Definition: DAIndex.H:32
p
volScalarField & p
Definition: createRefsPimpleDyM.H:6
DAResidualTurboFoam.H
setResidualClassMemberVector
#define setResidualClassMemberVector(stateName, stateUnit)
Definition: DAMacroFunctions.H:68
Foam::DAResidualTurboFoam::MRF_
IOMRFZoneListDF & MRF_
Multiple Reference Frame.
Definition: DAResidualTurboFoam.H:74
Foam::DAModel
Definition: DAModel.H:57
Foam::DAResidualTurboFoam::T_
volScalarField & T_
Definition: DAResidualTurboFoam.H:43
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
phi
surfaceScalarField & phi
Definition: createRefsPimpleDyM.H:8
Foam
Definition: checkGeometry.C:32
T
volScalarField & T
Definition: createRefsHeatTransfer.H:5
tUEqn
tmp< fvVectorMatrix > tUEqn(fvm::ddt(U)+fvm::div(phi, U)+turbulencePtr_->divDevReff(U) - fvSource)
Foam::DATurbulenceModel::devRhoReff
tmp< volSymmTensorField > devRhoReff() const
dev terms
Definition: DATurbulenceModel.C:371
AtU
volScalarField AtU(AU - UEqn.H1())
Foam::IOMRFZoneListDF
Definition: IOMRFZoneListDF.H:49
Foam::DAResidualTurboFoam::alphat_
volScalarField & alphat_
Definition: DAResidualTurboFoam.H:60
Foam::MRFZoneListDF::makeRelative
void makeRelative(volVectorField &U) const
Definition: MRFZoneListDF.C:148
adjustPhi
adjustPhi(phiHbyA, U, p)
alphaEff
volScalarField alphaEff("alphaEff", turbulencePtr_->nu()/Pr+alphat)
Foam::DAResidual
Definition: DAResidual.H:36
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:52
Foam::DAResidualTurboFoam::molWeight_
scalar molWeight_
Definition: DAResidualTurboFoam.H:78
Foam::DAResidualTurboFoam::correctBoundaryConditions
virtual void correctBoundaryConditions()
update the boundary condition for all the states in the selected solver
Definition: DAResidualTurboFoam.C:309
Foam::DAResidualTurboFoam::pressureControl_
pressureControl pressureControl_
pressureControl object to set ref pressure
Definition: DAResidualTurboFoam.H:71
rAU
volScalarField rAU(1.0/UEqn.A())
Foam::DAResidualTurboFoam::updateIntermediateVariables
virtual void updateIntermediateVariables()
update any intermediate variables that are dependent on state variables and are used in calcResiduals
Definition: DAResidualTurboFoam.C:233
Foam::DAResidualTurboFoam::phiRes_
surfaceScalarField phiRes_
Definition: DAResidualTurboFoam.H:47
Foam::MRFZoneListDF::correctBoundaryVelocity
void correctBoundaryVelocity(volVectorField &U) const
Definition: MRFZoneListDF.C:207
Foam::DAResidualTurboFoam::he_
volScalarField & he_
Definition: DAResidualTurboFoam.H:58
Foam::DAResidualTurboFoam::rho_
volScalarField & rho_
Definition: DAResidualTurboFoam.H:59
Foam::DAResidualTurboFoam::thermo_
fluidThermo & thermo_
thermophysical property
Definition: DAResidualTurboFoam.H:51
setResidualClassMemberScalar
#define setResidualClassMemberScalar(stateName, stateUnit)
Definition: DAMacroFunctions.H:53
Foam::DAResidualTurboFoam::U_
volVectorField & U_
Definition: DAResidualTurboFoam.H:37
Foam::DAResidualTurboFoam::Cp_
scalar Cp_
Definition: DAResidualTurboFoam.H:79
normalizePhiResiduals
#define normalizePhiResiduals(resName)
Definition: DAMacroFunctions.H:37
Foam::DAResidualTurboFoam::URel_
volVectorField & URel_
relative velocity
Definition: DAResidualTurboFoam.H:54
Teff
volSymmTensorField Teff
Definition: EEqnTurbo.H:3
Foam::DAResidualTurboFoam::psi_
volScalarField & psi_
Definition: DAResidualTurboFoam.H:61
EEqn
fvScalarMatrix EEqn(fvm::ddt(rho, he)+fvm::div(phi, he)+fvc::ddt(rho, K)+fvc::div(phi, K)+(he.name()=="e" ? fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)") :-dpdt) - fvm::laplacian(turbulencePtr_->alphaEff(), he) - fvSourceEnergy)