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