DAResidualSimpleFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAResidualSimpleFoam.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAResidualSimpleFoam, 0);
16 addToRunTimeSelectionTable(DAResidual, DAResidualSimpleFoam, 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(0, 1, -2, 0, 0, 0, 0)),
28  setResidualClassMemberScalar(p, dimensionSet(0, 0, -1, 0, 0, 0, 0)),
30  TResPtr_(nullptr),
31  alphaPorosity_(const_cast<volScalarField&>(
32  mesh_.thisDb().lookupObject<volScalarField>("alphaPorosity"))),
33  fvSource_(const_cast<volVectorField&>(
34  mesh_.thisDb().lookupObject<volVectorField>("fvSource"))),
35  fvOptions_(fv::options::New(mesh)),
36  daTurb_(const_cast<DATurbulenceModel&>(daModel.getDATurbulenceModel())),
37  // create simpleControl
38  simple_(const_cast<fvMesh&>(mesh)),
39  MRF_(const_cast<IOMRFZoneListDF&>(
40  mesh_.thisDb().lookupObject<IOMRFZoneListDF>("MRFProperties")))
41 {
42  // initialize fvSource
43  const dictionary& allOptions = daOption.getAllOptions();
44  if (allOptions.subDict("fvSource").toc().size() != 0)
45  {
46  hasFvSource_ = 1;
47  }
48 
49  // check whether to include the temperature field
50  hasTField_ = DAUtility::isFieldReadable(mesh, "T", "volScalarField");
51  if (hasTField_)
52  {
53 
54  TResPtr_.reset(new volScalarField(
55  IOobject(
56  "TRes",
57  mesh.time().timeName(),
58  mesh,
59  IOobject::NO_READ,
60  IOobject::NO_WRITE),
61  mesh,
62  dimensionedScalar("TRes", dimensionSet(0, 0, -1, 1, 0, 0, 0), 0.0),
63  zeroGradientFvPatchField<scalar>::typeName));
64 
65  // initialize the Prandtl number from transportProperties
66  IOdictionary transportProperties(
67  IOobject(
68  "transportProperties",
69  mesh.time().constant(),
70  mesh,
71  IOobject::MUST_READ,
72  IOobject::NO_WRITE,
73  false));
74  Pr_ = readScalar(transportProperties.lookup("Pr"));
75  Prt_ = readScalar(transportProperties.lookup("Prt"));
76  }
77 
78  // this is just a dummy call because we need to run the constrain once
79  // to initialize fvOptions, before we can use it. Otherwise, we may get
80  // a seg fault when we call fvOptions_.correct(U_) in updateIntermediateVars
81  fvVectorMatrix UEqn(
82  fvm::div(phi_, U_)
83  - fvOptions_(U_));
84  fvOptions_.constrain(UEqn);
85 }
86 
87 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
88 
90 {
91  /*
92  Description:
93  Clear all members to avoid memory leak because we will initalize
94  multiple objects of DAResidual. Here we need to delete all members
95  in the parent and child classes
96  */
97  URes_.clear();
98  pRes_.clear();
99  phiRes_.clear();
100  if (hasTField_)
101  {
102  TResPtr_->clear();
103  }
104 }
105 
106 void DAResidualSimpleFoam::calcResiduals(const dictionary& options)
107 {
108  /*
109  Description:
110  This is the function to compute residuals.
111 
112  Input:
113  options.isPC: 1 means computing residuals for preconditioner matrix.
114  This essentially use the first order scheme for div(phi,U)
115 
116  p_, U_, phi_, etc: State variables in OpenFOAM
117 
118  Output:
119  URes_, pRes_, phiRes_: residual field variables
120  */
121 
122  // ******** U Residuals **********
123  // copied and modified from UEqn.H
124 
125  word divUScheme = "div(phi,U)";
126 
127  label isPC = options.getLabel("isPC");
128 
129  if (isPC)
130  {
131  divUScheme = "div(pc)";
132  }
133 
134  if (hasFvSource_)
135  {
136  DAFvSource& daFvSource(const_cast<DAFvSource&>(
137  mesh_.thisDb().lookupObject<DAFvSource>("DAFvSource")));
138  daFvSource.calcFvSource(fvSource_);
139  }
140 
141  tmp<fvVectorMatrix> tUEqn(
142  fvm::div(phi_, U_, divUScheme)
143  + fvm::Sp(alphaPorosity_, U_)
144  + MRF_.DDt(U_)
146  - fvSource_
147  - fvOptions_(U_));
148  fvVectorMatrix& UEqn = tUEqn.ref();
149 
150  UEqn.relax();
151 
152  fvOptions_.constrain(UEqn);
153 
154  URes_ = (UEqn & U_) + fvc::grad(p_);
155  normalizeResiduals(URes);
156 
157  // ******** p Residuals **********
158  // copied and modified from pEqn.H
159  // NOTE manually set pRefCell and pRefValue
160  label pRefCell = 0;
161  scalar pRefValue = 0.0;
162 
163  volScalarField rAU(1.0 / UEqn.A());
164  //***************** NOTE *******************
165  // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
166  // because constraining variables will create discontinuity. Here we have a option to use the old
167  // implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
168  autoPtr<volVectorField> HbyAPtr = nullptr;
169  label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
170  if (useConstrainHbyA)
171  {
172  HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U_, p_)));
173  }
174  else
175  {
176  HbyAPtr.reset(new volVectorField("HbyA", U_));
177  HbyAPtr() = rAU * UEqn.H();
178  }
179  volVectorField& HbyA = HbyAPtr();
180 
181  surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
182 
184 
185  adjustPhi(phiHbyA, U_, p_);
186 
187  tmp<volScalarField> rAtU(rAU);
188 
189  if (simple_.consistent())
190  {
191  rAtU = 1.0 / (1.0 / rAU - UEqn.H1());
192  phiHbyA += fvc::interpolate(rAtU() - rAU) * fvc::snGrad(p_) * mesh_.magSf();
193  HbyA -= (rAU - rAtU()) * fvc::grad(p_);
194  }
195 
196  tUEqn.clear();
197 
198  // Update the pressure BCs to ensure flux consistency
200 
201  fvScalarMatrix pEqn(
202  fvm::laplacian(rAtU(), p_)
203  == fvc::div(phiHbyA));
204  pEqn.setReference(pRefCell, pRefValue);
205 
206  pRes_ = pEqn & p_;
207  normalizeResiduals(pRes);
208 
209  // ******** phi Residuals **********
210  // copied and modified from pEqn.H
211  phiRes_ = phiHbyA - pEqn.flux() - phi_;
212  // need to normalize phiRes
213  normalizePhiResiduals(phiRes);
214 
215  if (hasTField_)
216  {
217  volScalarField& alphat = const_cast<volScalarField&>(
218  mesh_.thisDb().lookupObject<volScalarField>("alphat"));
219 
220  volScalarField& T = const_cast<volScalarField&>(
221  mesh_.thisDb().lookupObject<volScalarField>("T"));
222 
223  volScalarField& TRes_ = TResPtr_();
224 
225  // ******** T Residuals **************
226  volScalarField alphaEff("alphaEff", daTurb_.nu() / Pr_ + alphat);
227 
228  fvScalarMatrix TEqn(
229  fvm::div(phi_, T)
230  - fvm::laplacian(alphaEff, T));
231 
232  TEqn.relax();
233 
234  TRes_ = TEqn & T;
235  normalizeResiduals(TRes);
236  }
237 }
238 
240 {
241  /*
242  Description:
243  Update the intermediate variables that depend on the state variables
244  */
245 
247  fvOptions_.correct(U_);
248 }
249 
251 {
252  /*
253  Description:
254  Update the boundary condition for all the states in the selected solver
255  */
256 
257  U_.correctBoundaryConditions();
258  p_.correctBoundaryConditions();
259  if (hasTField_)
260  {
261  volScalarField& T = const_cast<volScalarField&>(
262  mesh_.thisDb().lookupObject<volScalarField>("T"));
263  T.correctBoundaryConditions();
264  }
265 }
266 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
267 
268 } // End namespace Foam
269 
270 // ************************************************************************* //
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
pRefValue
scalar & pRefValue
Definition: createRefsPimpleDyM.H:12
Foam::DAResidualSimpleFoam::calcResiduals
virtual void calcResiduals(const dictionary &options)
compute residual
Definition: DAResidualSimpleFoam.C:106
Foam::DAResidualSimpleFoam::alphaPorosity_
volScalarField & alphaPorosity_
alpha porosity term
Definition: DAResidualSimpleFoam.H:48
DAResidualSimpleFoam.H
Foam::DAResidualSimpleFoam::updateIntermediateVariables
virtual void updateIntermediateVariables()
update any intermediate variables that are dependent on state variables and are used in calcResiduals
Definition: DAResidualSimpleFoam.C:239
Foam::DAResidualSimpleFoam::p_
volScalarField & p_
Definition: DAResidualSimpleFoam.H:38
useConstrainHbyA
label useConstrainHbyA
Definition: pEqnPimpleDyM.H:7
Foam::DAResidual::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAResidual.H:49
Foam::DATurbulenceModel::divDevReff
tmp< fvVectorMatrix > divDevReff(volVectorField &U)
divDev terms
Definition: DATurbulenceModel.C:411
Foam::DAOption
Definition: DAOption.H:29
Foam::DAResidualSimpleFoam::correctBoundaryConditions
virtual void correctBoundaryConditions()
update the boundary condition for all the states in the selected solver
Definition: DAResidualSimpleFoam.C:250
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
UEqn
fvVectorMatrix & UEqn
Definition: UEqnPimpleDyM.H:15
Foam::DAResidualSimpleFoam::hasFvSource_
label hasFvSource_
whether to have fvSource term
Definition: DAResidualSimpleFoam.H:66
Foam::DAFvSource::calcFvSource
virtual void calcFvSource(volVectorField &fvSource)
compute the FvSource term
Definition: DAFvSource.C:90
Foam::DAResidualSimpleFoam::URes_
volVectorField URes_
Definition: DAResidualSimpleFoam.H:36
Foam::DAResidualSimpleFoam::fvSource_
volVectorField & fvSource_
fvSource term
Definition: DAResidualSimpleFoam.H:51
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::DAResidualSimpleFoam::simple_
simpleControl simple_
simpleControl object which will be initialized in this class
Definition: DAResidualSimpleFoam.H:60
Foam::DAResidualSimpleFoam::hasTField_
label hasTField_
whether to include the temperature field
Definition: DAResidualSimpleFoam.H:69
pRefCell
label & pRefCell
Definition: createRefsPimpleDyM.H:11
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
phiHbyA
phiHbyA
Definition: pEqnRhoSimpleC.H:36
Foam::DAIndex
Definition: DAIndex.H:32
p
volScalarField & p
Definition: createRefsPimpleDyM.H:6
setResidualClassMemberVector
#define setResidualClassMemberVector(stateName, stateUnit)
Definition: DAMacroFunctions.H:68
Foam::DAResidualSimpleFoam::DAResidualSimpleFoam
DAResidualSimpleFoam(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAResidualSimpleFoam.C:19
Foam::DAResidualSimpleFoam::MRF_
IOMRFZoneListDF & MRF_
Multiple Reference Frame.
Definition: DAResidualSimpleFoam.H:63
Foam::DAResidualSimpleFoam::Prt_
scalar Prt_
Definition: DAResidualSimpleFoam.H:73
Foam::DAModel
Definition: DAModel.H:57
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
phi
surfaceScalarField & phi
Definition: createRefsPimpleDyM.H:8
Foam::DAResidualSimpleFoam::TResPtr_
autoPtr< volScalarField > TResPtr_
Definition: DAResidualSimpleFoam.H:44
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)
rAtU
tmp< volScalarField > rAtU(rAU)
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
Foam::DATurbulenceModel::nu
tmp< volScalarField > nu() const
get the nu field
Definition: DATurbulenceModel.C:267
adjustPhi
adjustPhi(phiHbyA, U, p)
alphaEff
volScalarField alphaEff("alphaEff", turbulencePtr_->nu()/Pr+alphat)
Foam::DAResidual
Definition: DAResidual.H:36
Foam::DAResidualSimpleFoam::Pr_
scalar Pr_
Definition: DAResidualSimpleFoam.H:71
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:52
alphat
volScalarField & alphat
Definition: TEqnPimpleDyM.H:5
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
Foam::DAUtility::isFieldReadable
static label isFieldReadable(const fvMesh &mesh, const word fieldName, const word fieldType)
Definition: DAUtility.C:817
HbyAPtr
autoPtr< volVectorField > HbyAPtr
Definition: pEqnPimpleDyM.H:6
Foam::DAResidualSimpleFoam::clear
virtual void clear()
clear the members
Definition: DAResidualSimpleFoam.C:89
Foam::DAResidualSimpleFoam::daTurb_
DATurbulenceModel & daTurb_
DATurbulenceModel object.
Definition: DAResidualSimpleFoam.H:57
Foam::DAResidualSimpleFoam::phiRes_
surfaceScalarField phiRes_
Definition: DAResidualSimpleFoam.H:42
setResidualClassMemberScalar
#define setResidualClassMemberScalar(stateName, stateUnit)
Definition: DAMacroFunctions.H:53
Foam::DAResidualSimpleFoam::U_
volVectorField & U_
Definition: DAResidualSimpleFoam.H:35
TEqn
fvScalarMatrix TEqn(fvm::ddt(T)+fvm::div(phi, T) - fvm::laplacian(alphaEff, T))
Foam::DAResidualSimpleFoam::fvOptions_
fv::options & fvOptions_
fvOptions term
Definition: DAResidualSimpleFoam.H:54
normalizePhiResiduals
#define normalizePhiResiduals(resName)
Definition: DAMacroFunctions.H:37
Foam::DAResidualSimpleFoam::phi_
surfaceScalarField & phi_
Definition: DAResidualSimpleFoam.H:41
Foam::DAResidualSimpleFoam::pRes_
volScalarField pRes_
Definition: DAResidualSimpleFoam.H:39