DAFunctionVariance.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAFunctionVariance.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAFunctionVariance, 0);
16 addToRunTimeSelectionTable(DAFunction, DAFunctionVariance, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAFunctionVariance::DAFunctionVariance(
20  const fvMesh& mesh,
21  const DAOption& daOption,
22  const DAModel& daModel,
23  const DAIndex& daIndex,
24  const word functionName)
25  : DAFunction(
26  mesh,
27  daOption,
28  daModel,
29  daIndex,
30  functionName),
31  daTurb_(const_cast<DATurbulenceModel&>(daModel.getDATurbulenceModel()))
32 {
33 
34  // read the parameters
35 
36  functionDict_.readEntry<word>("mode", mode_);
37 
38  if (mode_ == "surface")
39  {
40  // check if the faceSources is set
41  label size = faceSources_.size();
42  reduce(size, sumOp<label>());
43  if (size == 0)
44  {
45  FatalErrorIn("") << "surface mode is used but patchToFace is not set!"
46  << abort(FatalError);
47  }
48  }
49 
50  if (mode_ == "probePoint")
51  {
52  functionDict_.readEntry<List<List<scalar>>>("probePointCoords", probePointCoords_);
53  }
54 
55  functionDict_.readEntry<word>("varName", varName_);
56 
57  if (varName_ == "wallHeatFlux")
58  {
59  if (daTurb_.getTurbModelType() == "incompressible")
60  {
61  // initialize the Prandtl number from transportProperties
62  IOdictionary transportProperties(
63  IOobject(
64  "transportProperties",
65  mesh.time().constant(),
66  mesh,
67  IOobject::MUST_READ,
68  IOobject::NO_WRITE,
69  false));
70  // for incompressible flow, we need to read Cp from transportProperties
71  Cp_ = readScalar(transportProperties.lookup("Cp"));
72  }
73  }
74 
75  functionDict_.readEntry<word>("varType", varType_);
76 
77  if (varType_ == "vector")
78  {
79  functionDict_.readEntry<labelList>("indices", indices_);
80  }
81 
82  timeDependentRefData_ = functionDict_.getLabel("timeDependentRefData");
83 
84  // get the time information
85  scalar endTime = mesh_.time().endTime().value();
86  scalar deltaT = mesh_.time().deltaT().value();
87  label nTimeSteps = round(endTime / deltaT);
88 
89  word checkRefDataFolder;
90  label nRefValueInstances;
91  if (timeDependentRefData_)
92  {
93  // check if we can find the ref data in the endTime folder
94  checkRefDataFolder = Foam::name(endTime);
95  nRefValueInstances = nTimeSteps;
96  }
97  else
98  {
99  // check if we can find the ref data in the 0 folder
100  checkRefDataFolder = Foam::name(0);
101  nRefValueInstances = 1;
102  }
103 
104  // check if the reference data files exist
105  isRefData_ = 1;
106  if (varType_ == "scalar")
107  {
108  volScalarField varData(
109  IOobject(
110  varName_ + "Data",
111  checkRefDataFolder,
112  mesh_,
113  IOobject::READ_IF_PRESENT,
114  IOobject::NO_WRITE),
115  mesh_,
116  dimensionedScalar("dummy", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1e16),
117  "zeroGradient");
118 
119  if (gSumMag(varData) > 1e16)
120  {
121  isRefData_ = 0;
122  }
123  }
124  else if (varType_ == "vector")
125  {
126  volVectorField varData(
127  IOobject(
128  varName_ + "Data",
129  checkRefDataFolder,
130  mesh_,
131  IOobject::READ_IF_PRESENT,
132  IOobject::NO_WRITE),
133  mesh_,
134  dimensionedVector("dummy", dimensionSet(0, 0, 0, 0, 0, 0, 0), {1e16, 1e16, 1e16}),
135  "zeroGradient");
136 
137  if (mag(gSumCmptMag(varData)) > 1e16)
138  {
139  isRefData_ = 0;
140  }
141  }
142  else
143  {
144  FatalErrorIn("") << "varType " << varType_ << " not supported!"
145  << "Options are: scalar or vector"
146  << abort(FatalError);
147  }
148 
149  // if no varData file found, we print out a warning
150  if (isRefData_ == 0)
151  {
152  Info << endl;
153  Info << "**************************************************************************** " << endl;
154  Info << "* WARNING! Can't find data files or can't find valid * " << endl;
155  Info << "* values in data files for the variance function " << varName_ << " * " << endl;
156  Info << "**************************************************************************** " << endl;
157  Info << endl;
158  }
159  else
160  {
161  // varData file found, we need to read in the ref values for all time instances
162 
163  refValue_.setSize(nRefValueInstances);
164 
165  // set refValue
166  if (varType_ == "scalar")
167  {
168  for (label n = 0; n < nRefValueInstances; n++)
169  {
170  word timeName;
171  if (timeDependentRefData_)
172  {
173  scalar t = (n + 1) * deltaT;
174  timeName = Foam::name(t);
175  }
176  else
177  {
178  timeName = Foam::name(0);
179  }
180 
181  volScalarField varData(
182  IOobject(
183  varName_ + "Data",
184  timeName,
185  mesh_,
186  IOobject::MUST_READ,
187  IOobject::NO_WRITE),
188  mesh_);
189 
190  nRefPoints_ = 0;
191 
192  if (mode_ == "probePoint")
193  {
194  probeCellIndex_.setSize(0);
195 
196  forAll(probePointCoords_, idxI)
197  {
198  point pointCoord = {probePointCoords_[idxI][0], probePointCoords_[idxI][1], probePointCoords_[idxI][2]};
199  label cellI = DAUtility::myFindCell(mesh_, pointCoord);
200  if (cellI >= 0)
201  {
202  probeCellIndex_.append(cellI);
203  refValue_[n].append(varData[cellI]);
204  nRefPoints_++;
205  }
206  }
207  }
208  else if (mode_ == "surface")
209  {
210  forAll(faceSources_, idxI)
211  {
212  const label& functionFaceI = faceSources_[idxI];
213  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
214  const label patchI = daIndex_.bFacePatchI[bFaceI];
215  const label faceI = daIndex_.bFaceFaceI[bFaceI];
216  refValue_[n].append(varData.boundaryField()[patchI][faceI]);
217  nRefPoints_++;
218  }
219  }
220  else if (mode_ == "field")
221  {
222  forAll(cellSources_, idxI)
223  {
224  label cellI = cellSources_[idxI];
225  refValue_[n].append(varData[cellI]);
226  nRefPoints_++;
227  }
228  }
229  else
230  {
231  FatalErrorIn("") << "mode " << mode_ << " not supported!"
232  << "Options are: probePoint, field, or surface"
233  << abort(FatalError);
234  }
235  }
236  }
237  else if (varType_ == "vector")
238  {
239  for (label n = 0; n < nRefValueInstances; n++)
240  {
241  word timeName;
242  if (timeDependentRefData_)
243  {
244  scalar t = (n + 1) * deltaT;
245  timeName = Foam::name(t);
246  }
247  else
248  {
249  timeName = Foam::name(0);
250  }
251 
252  volVectorField varData(
253  IOobject(
254  varName_ + "Data",
255  timeName,
256  mesh_,
257  IOobject::MUST_READ,
258  IOobject::NO_WRITE),
259  mesh_);
260 
261  nRefPoints_ = 0;
262 
263  if (mode_ == "probePoint")
264  {
265  probeCellIndex_.setSize(0);
266 
267  forAll(probePointCoords_, idxI)
268  {
269  point pointCoord = {probePointCoords_[idxI][0], probePointCoords_[idxI][1], probePointCoords_[idxI][2]};
270  label cellI = DAUtility::myFindCell(mesh_, pointCoord);
271  if (cellI >= 0)
272  {
273  probeCellIndex_.append(cellI);
274  forAll(indices_, idxJ)
275  {
276  label compI = indices_[idxJ];
277  refValue_[n].append(varData[cellI][compI]);
278  nRefPoints_++;
279  }
280  }
281  }
282  }
283  else if (mode_ == "surface")
284  {
285  forAll(faceSources_, idxI)
286  {
287  const label& functionFaceI = faceSources_[idxI];
288  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
289  const label patchI = daIndex_.bFacePatchI[bFaceI];
290  const label faceI = daIndex_.bFaceFaceI[bFaceI];
291 
292  forAll(indices_, idxJ)
293  {
294  label compI = indices_[idxJ];
295  refValue_[n].append(varData.boundaryField()[patchI][faceI][compI]);
296  nRefPoints_++;
297  }
298  }
299  }
300  else if (mode_ == "field")
301  {
302  forAll(cellSources_, idxI)
303  {
304  label cellI = cellSources_[idxI];
305  forAll(indices_, idxJ)
306  {
307  label compI = indices_[idxJ];
308  refValue_[n].append(varData[cellI][compI]);
309  nRefPoints_++;
310  }
311  }
312  }
313  else
314  {
315  FatalErrorIn("") << "mode " << mode_ << " not supported!"
316  << "Options are: probePoint, field, or surface"
317  << abort(FatalError);
318  }
319  }
320  }
321 
322  reduce(nRefPoints_, sumOp<label>());
323 
324  Info << "Find " << nRefPoints_ << " reference points for variance of " << varName_ << endl;
325  if (nRefPoints_ == 0)
326  {
327  FatalErrorIn("") << "varData field exists but one can not find any valid data!"
328  << abort(FatalError);
329  }
330  }
331 }
332 
334 scalar DAFunctionVariance::calcFunction()
335 {
336  /*
337  Description:
338  Calculate the variance
339  */
340 
341  // initialize objFunValue
342  scalar functionValue = 0.0;
343 
344  if (isRefData_)
345  {
346 
347  const objectRegistry& db = mesh_.thisDb();
348 
349  label timeIndex;
350  if (timeDependentRefData_)
351  {
352  timeIndex = mesh_.time().timeIndex();
353  }
354  else
355  {
356  timeIndex = 1;
357  }
358 
359  if (varName_ == "wallShearStress")
360  {
361  volSymmTensorField devRhoReff = daTurb_.devRhoReff();
362 
363  label pointI = 0;
364  forAll(faceSources_, idxI)
365  {
366  const label& functionFaceI = faceSources_[idxI];
367  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
368  const label patchI = daIndex_.bFacePatchI[bFaceI];
369  const label faceI = daIndex_.bFaceFaceI[bFaceI];
370 
371  const vectorField& SfB = mesh_.Sf().boundaryField()[patchI];
372  const scalarField& magSfB = mesh_.magSf().boundaryField()[patchI];
373  const symmTensorField& ReffB = devRhoReff.boundaryField()[patchI];
374 
375  vectorField shearB = (-SfB / magSfB) & ReffB;
376 
377  forAll(indices_, idxJ)
378  {
379  label compI = indices_[idxJ];
380  scalar varDif = (shearB[faceI][compI] - refValue_[timeIndex - 1][pointI]);
381  functionValue += scale_ * varDif * varDif;
382  pointI++;
383  }
384  }
385  }
386  else if (varName_ == "wallHeatFlux")
387  {
388  if (daTurb_.getTurbModelType() == "incompressible")
389  {
390  // incompressible flow does not have he, so we do H = Cp * alphaEff * dT/dz
391  const volScalarField& T = mesh_.thisDb().lookupObject<volScalarField>("T");
392  volScalarField alphaEff = daTurb_.alphaEff();
393  const volScalarField::Boundary& TBf = T.boundaryField();
394  const volScalarField::Boundary& alphaEffBf = alphaEff.boundaryField();
395 
396  label pointI = 0;
397  forAll(faceSources_, idxI)
398  {
399  const label& functionFaceI = faceSources_[idxI];
400  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
401  const label patchI = daIndex_.bFacePatchI[bFaceI];
402  const label faceI = daIndex_.bFaceFaceI[bFaceI];
403 
404  scalarField hfx = Cp_ * alphaEffBf[patchI] * TBf[patchI].snGrad();
405 
406  scalar varDif = (hfx[faceI] - refValue_[timeIndex - 1][pointI]);
407  functionValue += scale_ * varDif * varDif;
408  pointI++;
409  }
410  }
411 
412  if (daTurb_.getTurbModelType() == "compressible")
413  {
414  // compressible flow, H = alphaEff * dHE/dz
415  fluidThermo& thermo_(const_cast<fluidThermo&>(
416  mesh_.thisDb().lookupObject<fluidThermo>("thermophysicalProperties")));
417  volScalarField& he = thermo_.he();
418  const volScalarField::Boundary& heBf = he.boundaryField();
419  volScalarField alphaEff = daTurb_.alphaEff();
420  const volScalarField::Boundary& alphaEffBf = alphaEff.boundaryField();
421 
422  label pointI = 0;
423  forAll(faceSources_, idxI)
424  {
425  const label& functionFaceI = faceSources_[idxI];
426  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
427  const label patchI = daIndex_.bFacePatchI[bFaceI];
428  const label faceI = daIndex_.bFaceFaceI[bFaceI];
429 
430  scalarField hfx = alphaEffBf[patchI] * heBf[patchI].snGrad();
431 
432  scalar varDif = (hfx[faceI] - refValue_[timeIndex - 1][pointI]);
433  functionValue += scale_ * varDif * varDif;
434  pointI++;
435  }
436  }
437  }
438  else
439  {
440  if (varType_ == "scalar")
441  {
442  const volScalarField& var = db.lookupObject<volScalarField>(varName_);
443 
444  if (mode_ == "probePoint")
445  {
446  forAll(probeCellIndex_, idxI)
447  {
448  label cellI = probeCellIndex_[idxI];
449  scalar varDif = (var[cellI] - refValue_[timeIndex - 1][idxI]);
450  functionValue += scale_ * varDif * varDif;
451  }
452  }
453  else if (mode_ == "surface")
454  {
455  label pointI = 0;
456  forAll(faceSources_, idxI)
457  {
458  const label& functionFaceI = faceSources_[idxI];
459  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
460  const label patchI = daIndex_.bFacePatchI[bFaceI];
461  const label faceI = daIndex_.bFaceFaceI[bFaceI];
462  scalar varDif = (var.boundaryField()[patchI][faceI] - refValue_[timeIndex - 1][pointI]);
463  functionValue += scale_ * varDif * varDif;
464  pointI++;
465  }
466  }
467  else if (mode_ == "field")
468  {
469  forAll(cellSources_, idxI)
470  {
471  label cellI = cellSources_[idxI];
472  scalar varDif = (var[cellI] - refValue_[timeIndex - 1][cellI]);
473  functionValue += scale_ * varDif * varDif;
474  }
475  }
476  }
477  else if (varType_ == "vector")
478  {
479  const volVectorField& var = db.lookupObject<volVectorField>(varName_);
480 
481  if (mode_ == "probePoint")
482  {
483  label pointI = 0;
484  forAll(probeCellIndex_, idxI)
485  {
486  label cellI = probeCellIndex_[idxI];
487  forAll(indices_, idxJ)
488  {
489  label compI = indices_[idxJ];
490  scalar varDif = (var[cellI][compI] - refValue_[timeIndex - 1][pointI]);
491  functionValue += scale_ * varDif * varDif;
492  pointI++;
493  }
494  }
495  }
496  else if (mode_ == "surface")
497  {
498  label pointI = 0;
499  forAll(faceSources_, idxI)
500  {
501  const label& functionFaceI = faceSources_[idxI];
502  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
503  const label patchI = daIndex_.bFacePatchI[bFaceI];
504  const label faceI = daIndex_.bFaceFaceI[bFaceI];
505  forAll(indices_, idxJ)
506  {
507  label compI = indices_[idxJ];
508  scalar varDif = (var.boundaryField()[patchI][faceI][compI] - refValue_[timeIndex - 1][pointI]);
509  functionValue += scale_ * varDif * varDif;
510  pointI++;
511  }
512  }
513  }
514  else if (mode_ == "field")
515  {
516  label pointI = 0;
517  forAll(cellSources_, idxI)
518  {
519  label cellI = cellSources_[idxI];
520  forAll(indices_, idxJ)
521  {
522  label compI = indices_[idxJ];
523  scalar varDif = (var[cellI][compI] - refValue_[timeIndex - 1][pointI]);
524  functionValue += scale_ * varDif * varDif;
525  pointI++;
526  }
527  }
528  }
529  }
530  else
531  {
532  FatalErrorIn("") << "varType " << varType_ << " not supported!"
533  << "Options are: scalar or vector"
534  << abort(FatalError);
535  }
536  }
537  // need to reduce the sum of force across all processors
538  reduce(functionValue, sumOp<scalar>());
539 
540  if (nRefPoints_ != 0)
541  {
542  functionValue /= nRefPoints_;
543  }
544  }
545 
546  return functionValue;
547 }
548 
549 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
550 
551 } // End namespace Foam
552 
553 // ************************************************************************* //
Foam::DAFunction::faceSources_
labelList faceSources_
a sorted list of all face sources for the objective function
Definition: DAFunction.H:67
Foam::DAFunctionVariance::indices_
List< label > indices_
components/indices of the vector variable
Definition: DAFunctionVariance.H:41
Foam::DAFunctionVariance::daTurb_
DATurbulenceModel & daTurb_
DATurbulenceModel object.
Definition: DAFunctionVariance.H:71
Foam::DAFunctionVariance::varName_
word varName_
name of the variable
Definition: DAFunctionVariance.H:32
Foam::DAIndex::bFacePatchI
labelList bFacePatchI
Definition: DAIndex.H:110
Foam::DAFunctionVariance::nRefPoints_
label nRefPoints_
total number of reference points
Definition: DAFunctionVariance.H:59
Foam::DAOption
Definition: DAOption.H:29
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
he
volScalarField & he
Definition: EEqnRhoPimple.H:5
Foam::DAFunctionVariance::timeDependentRefData_
label timeDependentRefData_
whether the ref data is time dependent if yes we need data in all time folders otherwise get it from ...
Definition: DAFunctionVariance.H:68
Foam::DAFunctionVariance::probeCellIndex_
labelList probeCellIndex_
cell index of the probe points
Definition: DAFunctionVariance.H:44
Foam::DAFunctionVariance::DAFunctionVariance
DAFunctionVariance(const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex, const word functionName)
Definition: DAFunctionVariance.C:19
Foam::DAFunctionVariance::refValue_
List< List< scalar > > refValue_
probe value
Definition: DAFunctionVariance.H:56
Foam::DAFunctionVariance::Cp_
scalar Cp_
Cp used in incompressible heatFlux calculation.
Definition: DAFunctionVariance.H:74
Foam::DAFunction::daIndex_
const DAIndex & daIndex_
DAIndex object.
Definition: DAFunction.H:52
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAFunctionVariance::calcFunction
virtual scalar calcFunction()
calculate the value of objective function
Definition: DAFunctionVariance.C:334
Foam::DAFunctionVariance::mode_
word mode_
data mode options are probePoint surface or field
Definition: DAFunctionVariance.H:35
Foam::DAModel
Definition: DAModel.H:57
Foam::DAFunction
Definition: DAFunction.H:31
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam
Definition: checkGeometry.C:32
T
volScalarField & T
Definition: createRefsHeatTransfer.H:5
Foam::DATurbulenceModel::devRhoReff
tmp< volSymmTensorField > devRhoReff() const
dev terms
Definition: DATurbulenceModel.C:371
Foam::DAFunctionVariance::isRefData_
label isRefData_
whether we find the reference data
Definition: DAFunctionVariance.H:65
alphaEff
volScalarField alphaEff("alphaEff", turbulencePtr_->nu()/Pr+alphat)
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DATurbulenceModel::getTurbModelType
word getTurbModelType() const
Definition: DATurbulenceModel.H:260
Foam::DAIndex::nLocalInternalFaces
label nLocalInternalFaces
local internal face size
Definition: DAIndex.H:98
Foam::DAFunctionVariance::probePointCoords_
List< List< scalar > > probePointCoords_
coordinates of the probe points
Definition: DAFunctionVariance.H:47
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:52
Foam::DAUtility::myFindCell
static label myFindCell(const primitiveMesh &mesh, const point &point)
Definition: DAUtility.C:803
Foam::DATurbulenceModel::alphaEff
tmp< volScalarField > alphaEff()
return effective thermal diffusivity
Definition: DATurbulenceModel.C:233
DAFunctionVariance.H
Foam::DAFunction::functionDict_
dictionary functionDict_
dictionary containing the information for the objective function
Definition: DAFunction.H:64
Foam::DAFunctionVariance::varType_
word varType_
type of the variable either vector or scalar
Definition: DAFunctionVariance.H:38
Foam::DAFunction::cellSources_
labelList cellSources_
a sorted list of all cell sources for the objective function
Definition: DAFunction.H:70
Foam::DAFunction::scale_
scalar scale_
scale of the objective function
Definition: DAFunction.H:73
Foam::DAFunction::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAFunction.H:43
Foam::DAIndex::bFaceFaceI
labelList bFaceFaceI
Definition: DAIndex.H:116

Generated by   doxygen 1.8.17