getFIData.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 Description
7  Extract the reference data for field inversion
8 
9 \*---------------------------------------------------------------------------*/
10 
11 #include "fvCFD.H"
12 #include "argList.H"
13 #include "Time.H"
14 #include "fvMesh.H"
15 #include "OFstream.H"
16 
17 using namespace Foam;
18 
19 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
20 
21 int main(int argc, char* argv[])
22 {
23 
24  Info << "Get ref data for field inversion...." << endl;
25 
26  argList::addOption(
27  "refFieldName",
28  "U",
29  "field variables name to extract");
30 
31  argList::addOption(
32  "refFieldType",
33  "vector",
34  "field variables type");
35 
36  argList::addOption(
37  "time",
38  "-1",
39  "prescribe a specific time to extract data");
40 
41 #include "setRootCase.H"
42 #include "createTime.H"
43 #include "createMesh.H"
44 
45  word refFieldName;
46  if (args.optionFound("refFieldName"))
47  {
48  refFieldName = word(args.optionLookup("refFieldName")());
49  }
50  else
51  {
52  Info << "refFieldName not set! Exit." << endl;
53  return 1;
54  }
55 
56  word refFieldType;
57  if (args.optionFound("refFieldType"))
58  {
59  refFieldType = word(args.optionLookup("refFieldType")());
60  }
61  else
62  {
63  Info << "refFieldType not set! Exit." << endl;
64  return 1;
65  }
66 
67  scalar time = -1.0;
68  if (args.optionFound("time"))
69  {
70  time = readScalar(args.optionLookup("time")());
71  if (time == 9999)
72  {
73  Info << "Extract latestTime" << endl;
74  }
75  else
76  {
77  Info << "Extract time = " << time << endl;
78  }
79  }
80  else
81  {
82  Info << "time not set! Extract all time instances." << endl;
83  }
84 
85  scalar endTime = runTime.endTime().value();
86  scalar deltaT = runTime.deltaT().value();
87  label nInstances = -1;
88  if (time == -1.0)
89  {
90  nInstances = round(endTime / deltaT);
91  }
92  else
93  {
94  nInstances = 1;
95  }
96 
97  Info << "Extracting field " << refFieldName << endl;
98 
99  for (label n = 1; n < nInstances + 1; n++)
100  {
101  scalar t = -1.0;
102  if (time == -1.0)
103  {
104  // read all times
105  t = n * deltaT;
106  }
107  else if (time == 9999)
108  {
109  // read from the latestTime (it is not necessarily the endTime)
110  instantList timeDirs = runTime.findTimes(runTime.path(), runTime.constant());
111  t = timeDirs.last().value();
112  }
113  else
114  {
115  // read from the specified time
116  t = time;
117  }
118  runTime.setTime(t, n);
119 
120  if (refFieldName == "wallHeatFlux")
121  {
122 
123  volScalarField hfx(
124  IOobject(
125  "wallHeatFluxData",
126  mesh.time().timeName(),
127  mesh,
128  IOobject::NO_READ,
129  IOobject::NO_WRITE),
130  mesh,
131  dimensionedScalar("wallHeatFluxData", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
132  "fixedValue");
133 
134  volScalarField rho(
135  IOobject(
136  "rho",
137  mesh.time().timeName(),
138  mesh,
139  IOobject::READ_IF_PRESENT,
140  IOobject::NO_WRITE),
141  mesh,
142  dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0, 0, 0), 1.0),
143  "zeroGradient");
144 
145  volScalarField alphat(
146  IOobject(
147  "alphat",
148  mesh.time().timeName(),
149  mesh,
150  IOobject::MUST_READ,
151  IOobject::NO_WRITE),
152  mesh);
153 
154  volScalarField T(
155  IOobject(
156  "T",
157  mesh.time().timeName(),
158  mesh,
159  IOobject::MUST_READ,
160  IOobject::NO_WRITE),
161  mesh);
162 
163  if (rho[0] == 1.0)
164  {
165  // incompressible cases
166 
167  volScalarField alpha(
168  IOobject(
169  "alpha",
170  mesh.time().timeName(),
171  mesh,
172  IOobject::NO_READ,
173  IOobject::NO_WRITE),
174  mesh,
175  dimensionedScalar("alpha", dimensionSet(0, 2, -1, 0, 0, 0, 0), 0.0),
176  "zeroGradient");
177 
178  IOdictionary transportProperties(
179  IOobject(
180  "transportProperties",
181  mesh.time().constant(),
182  mesh,
183  IOobject::MUST_READ,
184  IOobject::NO_WRITE,
185  false));
186 
187  scalar nu = transportProperties.getScalar("nu");
188  scalar Pr = transportProperties.getScalar("Pr");
189  scalar Cp = transportProperties.getScalar("Cp");
190  forAll(alpha, cellI)
191  {
192  alpha[cellI] = nu / Pr;
193  }
194  alpha.correctBoundaryConditions();
195 
196  volScalarField alphaEff("alphaEff", alphat + alpha);
197 
198  forAll(mesh.boundaryMesh(), patchI)
199  {
200  if (mesh.boundaryMesh()[patchI].type() == "wall" && mesh.boundaryMesh()[patchI].size() > 0)
201  {
202  scalarField flux = Cp * alphaEff.boundaryField()[patchI] * T.boundaryField()[patchI].snGrad();
203  forAll(mesh.boundaryMesh()[patchI], faceI)
204  {
205  hfx.boundaryFieldRef()[patchI][faceI] = flux[faceI];
206  }
207  }
208  }
209 
210  hfx.write();
211  }
212  else
213  {
214  // compressible
215 
216  volScalarField alpha(
217  IOobject(
218  "alpha",
219  mesh.time().timeName(),
220  mesh,
221  IOobject::NO_READ,
222  IOobject::NO_WRITE),
223  mesh,
224  dimensionedScalar("alpha", dimensionSet(1, -1, -1, 0, 0, 0, 0), 0.0),
225  "zeroGradient");
226 
227  IOdictionary thermophysicalProperties(
228  IOobject(
229  "thermophysicalProperties",
230  mesh.time().constant(),
231  mesh,
232  IOobject::MUST_READ,
233  IOobject::NO_WRITE,
234  false));
235 
236  dictionary subDict = thermophysicalProperties.subDict("mixture").subDict("transport");
237  scalar mu = subDict.getScalar("mu");
238  scalar Pr = subDict.getScalar("Pr");
239  scalar Cp = thermophysicalProperties.subDict("mixture").subDict("thermodynamics").getScalar("Cp");
240  scalar molWeight = thermophysicalProperties.subDict("mixture").subDict("specie").getScalar("molWeight");
241  forAll(alpha, cellI)
242  {
243  alpha[cellI] = mu * rho[cellI] / Pr;
244  }
245  alpha.correctBoundaryConditions();
246 
247  volScalarField alphaEff("alphaEff", alphat + alpha);
248 
249  volScalarField he(
250  IOobject(
251  "he",
252  mesh.time().timeName(),
253  mesh,
254  IOobject::NO_READ,
255  IOobject::NO_WRITE),
256  mesh,
257  dimensionedScalar("he", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
258  "fixedValue");
259 
260  word energy = thermophysicalProperties.subDict("thermoType").getWord("energy");
261  if (energy == "sensibleInternalEnergy")
262  {
263  // e
264  scalar RR = Foam::constant::thermodynamic::RR;
265  scalar R = RR / molWeight;
266  forAll(he, cellI)
267  {
268  he[cellI] = (Cp - R) * T[cellI];
269  }
270  forAll(he.boundaryField(), patchI)
271  {
272  forAll(he.boundaryField()[patchI], faceI)
273  {
274  he.boundaryFieldRef()[patchI][faceI] = (Cp - R) * T.boundaryField()[patchI][faceI];
275  }
276  }
277  he.correctBoundaryConditions();
278  }
279  else if (energy == "sensibleEnthalpy")
280  {
281  // h
282  forAll(he, cellI)
283  {
284  he[cellI] = Cp * T[cellI];
285  }
286  forAll(he.boundaryField(), patchI)
287  {
288  forAll(he.boundaryField()[patchI], faceI)
289  {
290  he.boundaryFieldRef()[patchI][faceI] = Cp * T.boundaryField()[patchI][faceI];
291  }
292  }
293  he.correctBoundaryConditions();
294  }
295  else
296  {
297  FatalErrorIn("") << "energy type: " << energy << " not supported in thermophysicalProperties" << abort(FatalError);
298  }
299 
300  forAll(mesh.boundaryMesh(), patchI)
301  {
302  if (mesh.boundaryMesh()[patchI].type() == "wall" && mesh.boundaryMesh()[patchI].size() > 0)
303  {
304  scalarField flux = alphaEff.boundaryField()[patchI] * he.boundaryField()[patchI].snGrad();
305  forAll(mesh.boundaryMesh()[patchI], faceI)
306  {
307  hfx.boundaryFieldRef()[patchI][faceI] = flux[faceI];
308  }
309  }
310  }
311 
312  hfx.write();
313  }
314  }
315  else if (refFieldName == "wallShearStress")
316  {
317 
318  // the following codes are adjusted from OpenFOAM/src/functionObjects/field/wallShearStress/wallShearStress.C
319 
320  volVectorField shear(
321  IOobject(
322  "wallShearStressData",
323  mesh.time().timeName(),
324  mesh,
325  IOobject::NO_READ,
326  IOobject::NO_WRITE),
327  mesh,
328  dimensionedVector("wallShearStressData", dimensionSet(0, 0, 0, 0, 0, 0, 0), vector::zero),
329  "fixedValue");
330 
331  volScalarField rho(
332  IOobject(
333  "rho",
334  mesh.time().timeName(),
335  mesh,
336  IOobject::READ_IF_PRESENT,
337  IOobject::NO_WRITE),
338  mesh,
339  dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0, 0, 0), 1.0),
340  "zeroGradient");
341 
342  volScalarField nut(
343  IOobject(
344  "nut",
345  mesh.time().timeName(),
346  mesh,
347  IOobject::MUST_READ,
348  IOobject::NO_WRITE),
349  mesh);
350 
351  volVectorField U(
352  IOobject(
353  "U",
354  mesh.time().timeName(),
355  mesh,
356  IOobject::MUST_READ,
357  IOobject::NO_WRITE),
358  mesh);
359 
360  volScalarField nu(
361  IOobject(
362  "nu",
363  mesh.time().timeName(),
364  mesh,
365  IOobject::NO_READ,
366  IOobject::NO_WRITE),
367  mesh,
368  dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0, 0, 0), 0.0),
369  "zeroGradient");
370 
371  if (rho[0] == 1.0)
372  {
373  // incompressible cases
374  IOdictionary transportProperties(
375  IOobject(
376  "transportProperties",
377  mesh.time().constant(),
378  mesh,
379  IOobject::MUST_READ,
380  IOobject::NO_WRITE,
381  false));
382 
383  scalar nuRead = transportProperties.getScalar("nu");
384  forAll(nu, cellI)
385  {
386  nu[cellI] = nuRead;
387  }
388  nu.correctBoundaryConditions();
389  }
390  else
391  {
392  // compressible
393  IOdictionary thermophysicalProperties(
394  IOobject(
395  "thermophysicalProperties",
396  mesh.time().constant(),
397  mesh,
398  IOobject::MUST_READ,
399  IOobject::NO_WRITE,
400  false));
401 
402  dictionary subDict = thermophysicalProperties.subDict("mixture").subDict("transport");
403  scalar mu = subDict.getScalar("mu");
404  forAll(nu, cellI)
405  {
406  nu[cellI] = mu * rho[cellI];
407  }
408  nu.correctBoundaryConditions();
409  }
410 
411  volScalarField nuEff("nuEff", nut + nu);
412 
413  volSymmTensorField devRhoReff = (-rho * nuEff) * dev(twoSymm(fvc::grad(U)));
414 
415  forAll(shear.boundaryField(), patchI)
416  {
417  vectorField& shearp = shear.boundaryFieldRef()[patchI];
418  const vectorField& Sfp = mesh.Sf().boundaryField()[patchI];
419  const scalarField& magSfp = mesh.magSf().boundaryField()[patchI];
420  const symmTensorField& Reffp = devRhoReff.boundaryField()[patchI];
421 
422  shearp = (-Sfp / magSfp) & Reffp;
423  }
424 
425  shear.write();
426  }
427  else if (refFieldType == "scalar")
428  {
429  volScalarField fieldRead(
430  IOobject(
431  refFieldName,
432  mesh.time().timeName(),
433  mesh,
434  IOobject::MUST_READ,
435  IOobject::NO_WRITE),
436  mesh);
437 
438  fieldRead.rename(refFieldName + "Data");
439  fieldRead.write();
440  }
441  else if (refFieldType == "vector")
442  {
443  volVectorField fieldRead(
444  IOobject(
445  refFieldName,
446  mesh.time().timeName(),
447  mesh,
448  IOobject::MUST_READ,
449  IOobject::NO_WRITE),
450  mesh);
451 
452  fieldRead.rename(refFieldName + "Data");
453  fieldRead.write();
454  }
455  else
456  {
457  Info << "refFieldName and type not supported!" << endl;
458  }
459  }
460 
461  Info << "Finished!" << endl;
462 
463  return 0;
464 }
465 
466 // ************************************************************************* //
mu
volScalarField & mu
Definition: createRefsSolidDisplacement.H:6
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Pr
dimensionedScalar Pr
Definition: TEqnPimpleDyM.H:3
he
volScalarField & he
Definition: EEqnRhoPimple.H:5
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam
Definition: checkGeometry.C:32
T
volScalarField & T
Definition: createRefsHeatTransfer.H:5
rho
volScalarField & rho
Definition: createRefsRhoPimple.H:8
argc
int argc
Definition: setArgs.H:18
alphaEff
volScalarField alphaEff("alphaEff", turbulencePtr_->nu()/Pr+alphat)
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
alphat
volScalarField & alphat
Definition: TEqnPimpleDyM.H:5
main
int main(int argc, char *argv[])
Definition: getFIData.C:21
args
Foam::argList & args
Definition: setRootCasePython.H:42
runTime
Time & runTime
Definition: createRefsHeatTransfer.H:1
argv
char ** argv
Definition: setArgs.H:26

Generated by   doxygen 1.8.17