DAObjFuncFieldInversion.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v3
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAObjFuncFieldInversion.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAObjFuncFieldInversion, 0);
16 addToRunTimeSelectionTable(DAObjFunc, DAObjFuncFieldInversion, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAObjFuncFieldInversion::DAObjFuncFieldInversion(
20  const fvMesh& mesh,
21  const DAOption& daOption,
22  const DAModel& daModel,
23  const DAIndex& daIndex,
24  const DAResidual& daResidual,
25  const word objFuncName,
26  const word objFuncPart,
27  const dictionary& objFuncDict)
28  : DAObjFunc(
29  mesh,
30  daOption,
31  daModel,
32  daIndex,
33  daResidual,
34  objFuncName,
35  objFuncPart,
36  objFuncDict),
37  daTurb_(daModel.getDATurbulenceModel())
38 
39 {
40  /*
41  Description:
42  Calculate the stateErrorNorm
43  f = scale * L2Norm( state-stateRef )
44  Input:
45  objFuncFaceSources: List of face source (index) for this objective
46 
47  objFuncCellSources: List of cell source (index) for this objective
48  Output:
49  objFuncFaceValues: the discrete value of objective for each face source (index).
50  This will be used for computing df/dw in the adjoint.
51 
52  objFuncCellValues: the discrete value of objective on each cell source (index).
53  This will be used for computing df/dw in the adjoint.
54 
55  objFuncValue: the sum of objective, reduced across all processors and scaled by "scale"
56  */
57 
58  objFuncDict_.readEntry<word>("type", objFuncType_);
59  data_ = objFuncDict_.getWord("data");
60  scale_ = objFuncDict_.getScalar("scale");
61  objFuncDict_.readEntry<bool>("weightedSum", weightedSum_);
62  if (weightedSum_ == true)
63  {
64  objFuncDict_.readEntry<scalar>("weight", weight_);
65  }
66 
67  // setup the connectivity, this is needed in Foam::DAJacCondFdW
68  // this objFunc only depends on the state variable at the zero level cell
69  if (DAUtility::isInList<word>(stateName_, daIndex.adjStateNames))
70  {
71  objFuncConInfo_ = {{stateName_}}; // level 0
72  }
73  else
74  {
75  objFuncConInfo_ = {{}}; // level 0
76  }
77 }
78 
80 void DAObjFuncFieldInversion::calcObjFunc(
81  const labelList& objFuncFaceSources,
82  const labelList& objFuncCellSources,
83  scalarList& objFuncFaceValues,
84  scalarList& objFuncCellValues,
85  scalar& objFuncValue)
86 {
87  /*
88  Description:
89  Calculate the stateErrorNorm
90  f = scale * L2Norm( state-stateRef )
91 
92  Input:
93  objFuncFaceSources: List of face source (index) for this objective
94 
95  objFuncCellSources: List of cell source (index) for this objective
96 
97  Output:
98  objFuncFaceValues: the discrete value of objective for each face source (index).
99  This will be used for computing df/dw in the adjoint.
100 
101  objFuncCellValues: the discrete value of objective on each cell source (index).
102  This will be used for computing df/dw in the adjoint.
103 
104  objFuncValue: the sum of objective, reduced across all processors and scaled by "scale"
105  */
106 
107  // initialize cell values to zero
108  forAll(objFuncCellValues, idxI)
109  {
110  objFuncCellValues[idxI] = 0.0;
111  }
112  // initialize objFunValue
113  objFuncValue = 0.0;
114 
115  const objectRegistry& db = mesh_.thisDb();
116 
117  if (data_ == "beta")
118  {
119  stateName_ = "betaFieldInversion";
120  const volScalarField betaFieldInversion_ = db.lookupObject<volScalarField>(stateName_);
121 
122  forAll(objFuncCellSources, idxI)
123  {
124  const label& cellI = objFuncCellSources[idxI];
125  objFuncCellValues[idxI] = scale_ * (sqr(betaFieldInversion_[cellI] - 1.0));
126  objFuncValue += objFuncCellValues[idxI];
127  }
128  // need to reduce the sum of all objectives across all processors
129  reduce(objFuncValue, sumOp<scalar>());
130 
131  if (weightedSum_ == true)
132  {
133  objFuncValue = weight_ * objFuncValue;
134  }
135  }
136  else if (data_ == "UData")
137  {
138  stateName_ = "U";
139  stateRefName_ = data_;
140  const volVectorField state = db.lookupObject<volVectorField>(stateName_);
141  const volVectorField stateRef = db.lookupObject<volVectorField>(stateRefName_);
142  forAll(objFuncCellSources, idxI)
143  {
144  const label& cellI = objFuncCellSources[idxI];
145  scalar stateRefCell(mag(stateRef[cellI]));
146  if (stateRefCell < 1e16)
147  {
148  objFuncCellValues[idxI] = (sqr(mag(scale_ * state[cellI] - stateRef[cellI])));
149  objFuncValue += objFuncCellValues[idxI];
150  }
151  }
152 
153  // need to reduce the sum of all objectives across all processors
154  reduce(objFuncValue, sumOp<scalar>());
155 
156  if (weightedSum_ == true)
157  {
158  objFuncValue = weight_ * objFuncValue;
159  }
160  }
161  else if (data_ == "USingleComponentData")
162  {
163  stateName_ = "U";
164  stateRefName_ = data_;
165  scalarList velocityCompt;
166  objFuncDict_.readEntry<scalarList>("velocityComponent", velocityCompt);
167  vector velocityComponent_;
168  velocityComponent_[0] = velocityCompt[0];
169  velocityComponent_[1] = velocityCompt[1];
170  velocityComponent_[2] = velocityCompt[2];
171 
172  // get the velocity field
173  const volVectorField state = db.lookupObject<volVectorField>(stateName_);
174 
175  // only use data for a specific component
176  const volScalarField stateRef = db.lookupObject<volScalarField>(stateRefName_);
177 
178  forAll(objFuncCellSources, idxI)
179  {
180  const label& cellI = objFuncCellSources[idxI];
181  if (stateRef[cellI] < 1e16)
182  {
183  vector URANS(state[cellI]);
184  scalar UData(stateRef[cellI]); // assume only using one component of velocity
185  objFuncCellValues[idxI] = (sqr(scale_ * (URANS & velocityComponent_) - UData));
186  objFuncValue += objFuncCellValues[idxI];
187  }
188  }
189 
190  // need to reduce the sum of all objectives across all processors
191  reduce(objFuncValue, sumOp<scalar>());
192 
193  if (weightedSum_ == true)
194  {
195  objFuncValue = weight_ * objFuncValue;
196  }
197 
198  }
199  else if (data_ == "pData")
200  {
201  stateName_ = "p";
202  stateRefName_ = data_;
203 
204  // get the pressure field
205  const volScalarField& state = db.lookupObject<volScalarField>(stateName_);
206 
207  // get the reference pressure field
208  const volScalarField& stateRef = db.lookupObject<volScalarField>(stateRefName_);
209 
210  forAll(objFuncCellSources, idxI)
211  {
212  const label& cellI = objFuncCellSources[idxI];
213  if (stateRef[cellI] < 1e16)
214  {
215  objFuncCellValues[idxI] = (sqr(scale_ * state[cellI]- stateRef[cellI]));
216  objFuncValue += objFuncCellValues[idxI];
217  }
218  }
219 
220  // need to reduce the sum of all objectives across all processors
221  reduce(objFuncValue, sumOp<scalar>());
222 
223  if (weightedSum_ == true)
224  {
225  objFuncValue = weight_ * objFuncValue;
226  }
227  }
228  else if (data_ == "surfacePressureData")
229  {
230  stateName_ = "p";
231  stateRefName_ = "pData";
232 
233  const volScalarField surfacePressureRef = db.lookupObject<volScalarField>(stateRefName_);
234  const volScalarField& p = db.lookupObject<volScalarField>(stateName_);
235 
236  wordList patchNames_;
237  objFuncDict_.readEntry<wordList>("patchNames", patchNames_);
238 
239  bool nonZeroPRefFlag_;
240  objFuncDict_.readEntry<bool>("nonZeroPRef", nonZeroPRefFlag_);
241 
242  scalar pRef_(0.0);
243  if (nonZeroPRefFlag_ == true)
244  {
245  scalarList pRefCoords;
246  objFuncDict_.readEntry<scalarList>("pRefCoords", pRefCoords);
247  vector pRefCoords_;
248  pRefCoords_[0] = pRefCoords[0];
249  pRefCoords_[1] = pRefCoords[1];
250  pRefCoords_[2] = pRefCoords[2];
251 
252  pRef_ = 0.0;
253  label cellID = mesh_.findCell(pRefCoords_);
254  // only assign pRef if the required cell is found in processor
255  if (cellID != -1)
256  {
257  pRef_ = p[cellID];
258  }
259  reduce(pRef_, maxOp<scalar>());
260  }
261 
262  forAll(patchNames_, cI)
263  {
264  label patchI = mesh_.boundaryMesh().findPatchID(patchNames_[cI]);
265  const fvPatch& patch = mesh_.boundary()[patchI];
266  forAll(patch, faceI)
267  {
268  scalar bSurfacePressure = scale_ * (p.boundaryField()[patchI][faceI] - pRef_);
269  scalar bSurfacePressureRef = surfacePressureRef.boundaryField()[patchI][faceI];
270  if (bSurfacePressureRef < 1e16)
271  {
272  objFuncValue += sqr(bSurfacePressure - bSurfacePressureRef);
273  }
274  }
275  }
276  // need to reduce the sum of all objectives across all processors
277  reduce(objFuncValue, sumOp<scalar>());
278  if (weightedSum_ == true)
279  {
280  objFuncValue = weight_ * objFuncValue;
281  }
282  }
283  else if(data_ == "surfaceFrictionData")
284  {
285  stateName_ = "surfaceFriction";
286  stateRefName_ = data_;
287  wordList patchNames_;
288  objFuncDict_.readEntry<wordList>("patchNames", patchNames_);
289 
290  scalarList dir;
291  vector wssDir_;
292  objFuncDict_.readEntry<scalarList>("wssDir", dir);
293  wssDir_[0] = dir[0];
294  wssDir_[1] = dir[1];
295  wssDir_[2] = dir[2];
296 
297  volScalarField& surfaceFriction = const_cast<volScalarField&>(db.lookupObject<volScalarField>(stateName_));
298  const volScalarField surfaceFrictionRef = db.lookupObject<volScalarField>(stateRefName_);
299 
300  // ingredients for the computation
301  tmp<volSymmTensorField> Reff = daTurb_.devRhoReff();
302  volSymmTensorField::Boundary bReff = Reff().boundaryField();
303  const surfaceVectorField::Boundary& Sfp = mesh_.Sf().boundaryField();
304  const surfaceScalarField::Boundary& magSfp = mesh_.magSf().boundaryField();
305 
306  forAll(patchNames_, cI)
307  {
308  label patchI = mesh_.boundaryMesh().findPatchID(patchNames_[cI]);
309  const fvPatch& patch = mesh_.boundary()[patchI];
310  forAll(patch, faceI)
311  {
312  scalar bSurfaceFrictionRef = surfaceFrictionRef.boundaryField()[patchI][faceI];
313 
314  // normal vector at wall, use -ve sign to ensure vector pointing into the domain
315  vector normal = -Sfp[patchI][faceI] / magSfp[patchI][faceI];
316 
317  // wall shear stress
318  vector wss = normal & bReff[patchI][faceI];
319 
320  // wallShearStress or surfaceFriction (use surfaceFriction label to match the fields)
321  scalar bSurfaceFriction = scale_ * (wss & wssDir_);
322 
323  surfaceFriction.boundaryFieldRef()[patchI][faceI] = bSurfaceFriction;
324 
325  // The following will allow to only use the Cf data at certain cells.
326  // If you want to exclude cells, then given them a cell value of 1e16.
327  if (bSurfaceFrictionRef < 1e16)
328  {
329  // calculate the objective function
330  objFuncValue += sqr(bSurfaceFriction - bSurfaceFrictionRef);
331  }
332  }
333  }
334 
335  // need to reduce the sum of all objectives across all processors
336  reduce(objFuncValue, sumOp<scalar>());
337 
338  if (weightedSum_ == true)
339  {
340  objFuncValue = weight_ * objFuncValue;
341  }
342  }
343  else if(data_ == "aeroCoeffData")
344  {
345  // get the ingredients for computations
346  scalarList dir; // decides if computing lift or drag
347  objFuncDict_.readEntry<scalarList>("direction", dir);
348  vector forceDir_;
349  forceDir_[0] = dir[0];
350  forceDir_[1] = dir[1];
351  forceDir_[2] = dir[2];
352  scalar aeroCoeffRef_;
353  objFuncDict_.readEntry<scalar>("aeroCoeffRef", aeroCoeffRef_);
354  wordList patchNames_;
355  objFuncDict_.readEntry<wordList>("patchNames", patchNames_);
356 
357  const volScalarField& p = db.lookupObject<volScalarField>("p");
358  const surfaceVectorField::Boundary& Sfb = mesh_.Sf().boundaryField();
359  tmp<volSymmTensorField> tdevRhoReff = daTurb_.devRhoReff();
360  const volSymmTensorField::Boundary& devRhoReffb = tdevRhoReff().boundaryField();
361 
362  vector forces(vector::zero);
363 
364  scalar aeroCoeff(0.0);
365 
366  forAll(patchNames_, cI)
367  {
368  // get the patch id label
369  label patchI = mesh_.boundaryMesh().findPatchID(patchNames_[cI]);
370 
371  // create a shorter handle for the boundary patch
372  const fvPatch& patch = mesh_.boundary()[patchI];
373 
374  // normal force
375  vectorField fN(Sfb[patchI] * p.boundaryField()[patchI]);
376 
377  // tangential force
378  vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
379 
380  forAll(patch, faceI)
381  {
382  forces.x() = fN[faceI].x() + fT[faceI].x();
383  forces.y() = fN[faceI].y() + fT[faceI].y();
384  forces.z() = fN[faceI].z() + fT[faceI].z();
385  aeroCoeff += scale_ * (forces & forceDir_);
386  }
387  }
388  // need to reduce the sum of all forces across all processors
389  reduce(aeroCoeff, sumOp<scalar>());
390 
391  // compute the objective function
392  objFuncValue += sqr(aeroCoeff - aeroCoeffRef_);
393 
394  // scale if performing weighted-sum multi-objective optimisation
395  if (weightedSum_ == true)
396  {
397  objFuncValue = weight_ * objFuncValue;
398  }
399  }
400  else if (data_ == "surfaceFrictionDataModified")
401  {
402  /*
403  In this modified implementation the wallShearStress computation different to the loop starting in line 283.
404  Here the wallShearStress is computed using the equation shown in line 454.
405  (This modified equation is useful because this how the Cf is defined for the popular periodic hill case in
406  literature.)
407  */
408 
409  stateName_ = "surfaceFriction_";
410  stateRefName_ = data_;
411 
412  wordList patchNames_;
413  objFuncDict_.readEntry<wordList>("patchNames", patchNames_);
414 
415  // get surface friction "fields"
416  volScalarField& surfaceFriction = const_cast<volScalarField&>(db.lookupObject<volScalarField>(stateName_));
417  const volScalarField& surfaceFrictionRef = db.lookupObject<volScalarField>(stateRefName_);
418 
419  // ingredients for surface friction computation
420  const volVectorField& U = db.lookupObject<volVectorField>("U");
421  tmp<volTensorField> gradU = fvc::grad(U);
422  const volTensorField::Boundary& bGradU = gradU().boundaryField();
423 
424  const surfaceVectorField::Boundary& Sfp = mesh_.Sf().boundaryField();
425  const surfaceScalarField::Boundary& magSfp = mesh_.magSf().boundaryField();
426 
427  forAll(patchNames_, cI)
428  {
429  label patchI = mesh_.boundaryMesh().findPatchID(patchNames_[cI]);
430  const fvPatch& patch = mesh_.boundary()[patchI];
431  forAll(patch, faceI)
432  {
433  // normal vector at wall, use -ve sign to ensure vector pointing into the domain
434  vector normal = -Sfp[patchI][faceI] / magSfp[patchI][faceI];
435 
436  // tangent vector, computed from normal: tangent = p x normal; where p is a unit vector perpidicular to the plane
437  // NOTE: make this more general
438  vector tangent(normal.y(), -normal.x(), 0.0);
439 
440  // velocity gradient at wall
441  tensor fGradU = bGradU[patchI][faceI];
442 
443  /* need to get transpose of fGradU as it is stored in the following form:
444  fGradU = [ XX XY XZ YX YY YZ ZX ZY ZZ ] (see slide 3 of: https://tinyurl.com/5ops2f5h)
445  but by definition,
446  grad(U) = [ XX YX ZX XY YY ZY XZ YZ ZZ]
447  where XX = partial(u1)/partial(x1), YX = partial(u2)/partial(x1) and so on.
448  => grad(U) = transpose(fGradU)
449  */
450  tensor fGradUT = fGradU.T();
451 
452  /* compute the surface friction assuming incompressible flow and no wall functions! Add run time warning message
453  to reflect this later.
454  Cf = tau_wall/dynPressure,
455  tau_wall = rho * nu * tangent . (grad(U) . normal),
456  dynPressue = 0.5 * rho * mag(U_bulk^2),
457  incorporate rho, mu, and dynamic pressure in scale_ as follows,
458  scale_ = nu / dynPressure,
459  => Cf = scale_ * tangent . (grad(U) . normal) */
460  scalar bSurfaceFriction = scale_ * (tangent & (fGradUT & normal));
461 
462  surfaceFriction.boundaryFieldRef()[patchI][faceI] = bSurfaceFriction;
463 
464  // calculate the objective function
465  // extract the reference surface friction at the boundary
466  scalar bSurfaceFrictionRef = surfaceFrictionRef.boundaryField()[patchI][faceI];
467 
468  if (bSurfaceFrictionRef < 1e16)
469  {
470  // calculate the objective function
471  objFuncValue += sqr(bSurfaceFriction - bSurfaceFrictionRef);
472  }
473  }
474  }
475 
476  // need to reduce the sum of all objectives across all processors
477  reduce(objFuncValue, sumOp<scalar>());
478 
479  if (weightedSum_ == true)
480  {
481  objFuncValue = weight_ * objFuncValue;
482  }
483  }
484  else
485  {
486  FatalErrorIn("") << "dataType: " << data_
487  << " not supported for field inversion! "
488  << "Available options are: UData, pData, surfacePressureData, surfaceFrictionData, aeroCoeffData, and surfaceFrictionDataPeriodicHill."
489  << abort(FatalError);
490  }
491 
492  return;
493 }
494 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
495 
496 } // End namespace Foam
497 
498 // ************************************************************************* //
Foam::DAObjFuncFieldInversion::weightedSum_
bool weightedSum_
weighted objective function (especially useful when performing weighted sum multi-objective optimisat...
Definition: DAObjFuncFieldInversion.H:42
Foam::DAObjFunc::objFuncType_
word objFuncType_
the type of the objective function
Definition: DAObjFunc.H:66
U
U
Definition: pEqnPimpleDyM.H:60
Foam::DAObjFuncFieldInversion::data_
word data_
name of the state variable to compute the error norm
Definition: DAObjFuncFieldInversion.H:33
forAll
forAll(pseudoP.boundaryField(), patchI)
Definition: solvePseudoPEqn.H:10
Foam::DAOption
Definition: DAOption.H:29
Foam::DAObjFuncFieldInversion::scale_
scalar scale_
Definition: DAObjFuncFieldInversion.H:39
Foam::DAObjFuncFieldInversion::DAObjFuncFieldInversion
DAObjFuncFieldInversion(const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex, const DAResidual &daResidual, const word objFuncName, const word objFuncPart, const dictionary &objFuncDict)
Definition: DAObjFuncFieldInversion.C:19
daOption
DAOption daOption(mesh, pyOptions_)
Foam::DAObjFuncFieldInversion::calcObjFunc
virtual void calcObjFunc(const labelList &objFuncFaceSources, const labelList &objFuncCellSources, scalarList &objFuncFaceValues, scalarList &objFuncCellValues, scalar &objFuncValue)
calculate the value of objective function
Definition: DAObjFuncFieldInversion.C:80
p
volScalarField & p
Definition: createRefsPimple.H:6
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAObjFuncFieldInversion::stateName_
word stateName_
Definition: DAObjFuncFieldInversion.H:35
DAObjFuncFieldInversion.H
Foam::DAObjFunc::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAObjFunc.H:45
Foam::DAObjFuncFieldInversion::stateRefName_
word stateRefName_
Definition: DAObjFuncFieldInversion.H:37
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAModel
Definition: DAModel.H:59
Foam
Definition: multiFreqScalarFvPatchField.C:144
Foam::DATurbulenceModel::devRhoReff
tmp< volSymmTensorField > devRhoReff() const
dev terms
Definition: DATurbulenceModel.C:306
Foam::DAObjFuncFieldInversion::weight_
scalar weight_
Definition: DAObjFuncFieldInversion.H:43
Foam::DAResidual
Definition: DAResidual.H:35
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFvSource, 0)
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFvSource, DAFvSourceActuatorDisk, dictionary)
daModel
DAModel daModel(mesh, daOption)
Foam::DAObjFunc::objFuncConInfo_
List< List< word > > objFuncConInfo_
the connectivity information for the objective function
Definition: DAObjFunc.H:93
daIndex
DAIndex daIndex(mesh, daOption, daModel)
Foam::DAObjFunc::objFuncDict_
const dictionary & objFuncDict_
dictionary containing the information for the objective function
Definition: DAObjFunc.H:69
Foam::DAObjFunc
Definition: DAObjFunc.H:33
Foam::DAObjFuncFieldInversion::daTurb_
const DATurbulenceModel & daTurb_
DATurbulenceModel object.
Definition: DAObjFuncFieldInversion.H:46

Generated by   doxygen 1.8.17