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