DAFunctionLocation.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAFunctionLocation.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAFunctionLocation, 0);
16 addToRunTimeSelectionTable(DAFunction, DAFunctionLocation, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAFunctionLocation::DAFunctionLocation(
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 {
32 
33  functionDict_.readEntry<word>("mode", mode_);
34 
35  if (functionDict_.found("coeffKS"))
36  {
37  functionDict_.readEntry<scalar>("coeffKS", coeffKS_);
38  }
39 
40  if (functionDict_.found("axis"))
41  {
42  scalarList axisRead;
43  functionDict_.readEntry<scalarList>("axis", axisRead);
44  axis_ = {axisRead[0], axisRead[1], axisRead[2]};
45  axis_ /= mag(axis_);
46  }
47 
48  if (functionDict_.found("center"))
49  {
50  scalarList centerRead;
51  functionDict_.readEntry<scalarList>("center", centerRead);
52  center_ = {centerRead[0], centerRead[1], centerRead[2]};
53  }
54 
55  snapCenter2Cell_ = functionDict_.lookupOrDefault<label>("snapCenter2Cell", 0);
56  if (snapCenter2Cell_)
57  {
58  point centerPoint = {center_[0], center_[1], center_[2]};
59 
60  // NOTE: we need to call a self-defined findCell func to make it work correctly in ADR
61  snappedCenterCellI_ = DAUtility::myFindCell(mesh_, centerPoint);
62 
63  label foundCellI = 0;
64  if (snappedCenterCellI_ >= 0)
65  {
66  foundCellI = 1;
67  }
68  reduce(foundCellI, sumOp<label>());
69  if (foundCellI != 1)
70  {
71  FatalErrorIn(" ") << "There should be only one cell found globally while "
72  << foundCellI << " was returned."
73  << " Please adjust center such that it is located completely"
74  << " within a cell in the mesh domain. The center should not "
75  << " be outside of the mesh domain or on a mesh face "
76  << abort(FatalError);
77  }
78  vector snappedCenter = vector::zero;
79  this->findGlobalSnappedCenter(snappedCenterCellI_, snappedCenter);
80  Info << "snap to center " << snappedCenter << endl;
81  }
82 
83  if (mode_ == "maxRadius")
84  {
85  // we need to identify the patchI and faceI that has maxR
86  // the proc that does not own the maxR face will have negative patchI and faceI
87  // we assume the patchI and faceI do not change during the optimization
88  // otherwise, we should use maxRadiusKS instead
89  scalar maxR = -100000;
90  label maxRPatchI = -999, maxRFaceI = -999;
91  forAll(faceSources_, idxI)
92  {
93  const label& functionFaceI = faceSources_[idxI];
94  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
95  const label patchI = daIndex_.bFacePatchI[bFaceI];
96  const label faceI = daIndex_.bFaceFaceI[bFaceI];
97 
98  vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center_;
99 
100  tensor faceCTensor(tensor::zero);
101  faceCTensor.xx() = faceC.x();
102  faceCTensor.yy() = faceC.y();
103  faceCTensor.zz() = faceC.z();
104 
105  vector faceCAxial = faceCTensor & axis_;
106  vector faceCRadial = faceC - faceCAxial;
107 
108  scalar radius = mag(faceCRadial);
109 
110  if (radius > maxR)
111  {
112  maxR = radius;
113  maxRPatchI = patchI;
114  maxRFaceI = faceI;
115  }
116  }
117 
118  scalar maxRGlobal = maxR;
119  reduce(maxRGlobal, maxOp<scalar>());
120 
121  if (fabs(maxRGlobal - maxR) < 1e-8)
122  {
123  maxRPatchI_ = maxRPatchI;
124  maxRFaceI_ = maxRFaceI;
125  }
126  }
127 }
128 
129 void DAFunctionLocation::findGlobalSnappedCenter(
130  label snappedCenterCellI,
131  vector& center)
132 {
133  scalar centerX = 0.0;
134  scalar centerY = 0.0;
135  scalar centerZ = 0.0;
136 
137  if (snappedCenterCellI >= 0)
138  {
139  centerX = mesh_.C()[snappedCenterCellI][0];
140  centerY = mesh_.C()[snappedCenterCellI][1];
141  centerZ = mesh_.C()[snappedCenterCellI][2];
142  }
143  reduce(centerX, sumOp<scalar>());
144  reduce(centerY, sumOp<scalar>());
145  reduce(centerZ, sumOp<scalar>());
146 
147  center[0] = centerX;
148  center[1] = centerY;
149  center[2] = centerZ;
150 }
151 
153 scalar DAFunctionLocation::calcFunction()
154 {
155  /*
156  Description:
157  Calculate the Location of a selected patch. The actual computation depends on the mode
158  */
159 
160  // initialize objFunValue
161  scalar functionValue = 0.0;
162 
163  if (mode_ == "maxRadiusKS")
164  {
165  // calculate Location
166  scalar objValTmp = 0.0;
167 
168  vector center = center_;
169  if (snapCenter2Cell_)
170  {
171  this->findGlobalSnappedCenter(snappedCenterCellI_, center);
172  }
173 
174  forAll(faceSources_, idxI)
175  {
176  const label& functionFaceI = faceSources_[idxI];
177  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
178  const label patchI = daIndex_.bFacePatchI[bFaceI];
179  const label faceI = daIndex_.bFaceFaceI[bFaceI];
180 
181  vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
182 
183  tensor faceCTensor(tensor::zero);
184  faceCTensor.xx() = faceC.x();
185  faceCTensor.yy() = faceC.y();
186  faceCTensor.zz() = faceC.z();
187 
188  vector faceCAxial = faceCTensor & axis_;
189  vector faceCRadial = faceC - faceCAxial;
190 
191  scalar radius = mag(faceCRadial);
192 
193  objValTmp += exp(coeffKS_ * radius);
194 
195  if (objValTmp > 1e200)
196  {
197  FatalErrorIn(" ") << "KS function summation term too large! "
198  << "Reduce coeffKS! " << abort(FatalError);
199  }
200  }
201 
202  // need to reduce the sum of force across all processors
203  reduce(objValTmp, sumOp<scalar>());
204 
205  functionValue = log(objValTmp) / coeffKS_;
206  }
207  else if (mode_ == "maxInverseRadiusKS")
208  {
209  // this is essentially minimal radius using KS
210 
211  // calculate Location
212  scalar objValTmp = 0.0;
213 
214  vector center = center_;
215  if (snapCenter2Cell_)
216  {
217  this->findGlobalSnappedCenter(snappedCenterCellI_, center);
218  }
219 
220  forAll(faceSources_, idxI)
221  {
222  const label& functionFaceI = faceSources_[idxI];
223  label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
224  const label patchI = daIndex_.bFacePatchI[bFaceI];
225  const label faceI = daIndex_.bFaceFaceI[bFaceI];
226 
227  vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
228 
229  tensor faceCTensor(tensor::zero);
230  faceCTensor.xx() = faceC.x();
231  faceCTensor.yy() = faceC.y();
232  faceCTensor.zz() = faceC.z();
233 
234  vector faceCAxial = faceCTensor & axis_;
235  vector faceCRadial = faceC - faceCAxial;
236 
237  scalar radius = mag(faceCRadial);
238  scalar iRadius = 1.0 / (radius + 1e-12);
239 
240  objValTmp += exp(coeffKS_ * iRadius);
241 
242  if (objValTmp > 1e200)
243  {
244  FatalErrorIn(" ") << "KS function summation term too large! "
245  << "Reduce coeffKS! " << abort(FatalError);
246  }
247  }
248 
249  // need to reduce the sum of force across all processors
250  reduce(objValTmp, sumOp<scalar>());
251 
252  functionValue = log(objValTmp) / coeffKS_;
253  }
254  else if (mode_ == "maxRadius")
255  {
256  scalar radius = 0.0;
257 
258  vector center = center_;
259  if (snapCenter2Cell_)
260  {
261  this->findGlobalSnappedCenter(snappedCenterCellI_, center);
262  }
263 
264  if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0)
265  {
266 
267  vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center;
268 
269  tensor faceCTensor(tensor::zero);
270  faceCTensor.xx() = faceC.x();
271  faceCTensor.yy() = faceC.y();
272  faceCTensor.zz() = faceC.z();
273 
274  vector faceCAxial = faceCTensor & axis_;
275  vector faceCRadial = faceC - faceCAxial;
276 
277  radius = mag(faceCRadial);
278  }
279 
280  reduce(radius, sumOp<scalar>());
281 
282  functionValue = radius;
283  }
284  else
285  {
286  FatalErrorIn("DAFunctionLocation") << "mode: " << mode_ << " not supported!"
287  << "Options are: maxRadius"
288  << abort(FatalError);
289  }
290 
291  // check if we need to calculate refDiff.
292  this->calcRefVar(functionValue);
293 
294  return functionValue;
295 }
296 
297 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
298 
299 } // End namespace Foam
300 
301 // ************************************************************************* //
Foam::DAFunction::faceSources_
labelList faceSources_
a sorted list of all face sources for the objective function
Definition: DAFunction.H:67
Foam::DAFunctionLocation::snapCenter2Cell_
label snapCenter2Cell_
whether to snap the center to a cell in the mesh if yes the center will move with the mesh
Definition: DAFunctionLocation.H:50
Foam::DAFunctionLocation::snappedCenterCellI_
label snappedCenterCellI_
the cell index for the center if snapCenter2Cell_ = 1
Definition: DAFunctionLocation.H:53
Foam::DAIndex::bFacePatchI
labelList bFacePatchI
Definition: DAIndex.H:110
Foam::DAFunctionLocation::center_
vector center_
center for radius computation
Definition: DAFunctionLocation.H:41
Foam::DAOption
Definition: DAOption.H:29
Foam::DAFunctionLocation::maxRPatchI_
label maxRPatchI_
maxRadius patchI
Definition: DAFunctionLocation.H:44
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
DAFunctionLocation.H
Foam::DAFunctionLocation::axis_
vector axis_
axis for radius computation
Definition: DAFunctionLocation.H:38
Foam::DAFunction::calcRefVar
void calcRefVar(scalar &functionValue)
calculate (var-ref)^2
Definition: DAFunction.C:216
Foam::DAFunction::daIndex_
const DAIndex & daIndex_
DAIndex object.
Definition: DAFunction.H:52
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAFunctionLocation::maxRFaceI_
label maxRFaceI_
maxRadius patchI
Definition: DAFunctionLocation.H:47
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAFunctionLocation::mode_
word mode_
location mode to use
Definition: DAFunctionLocation.H:35
Foam::DAFunctionLocation::calcFunction
virtual scalar calcFunction()
calculate the value of objective function
Definition: DAFunctionLocation.C:153
Foam::DAFunctionLocation::findGlobalSnappedCenter
void findGlobalSnappedCenter(label snappedCenterCellI, vector &center)
Definition: DAFunctionLocation.C:129
Foam::DAModel
Definition: DAModel.H:57
Foam::DAFunction
Definition: DAFunction.H:31
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam
Definition: checkGeometry.C:32
Foam::DAFunctionLocation::DAFunctionLocation
DAFunctionLocation(const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex, const word functionName)
Definition: DAFunctionLocation.C:19
Foam::DAFunctionLocation::coeffKS_
scalar coeffKS_
coefficient for the KS function
Definition: DAFunctionLocation.H:32
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAIndex::nLocalInternalFaces
label nLocalInternalFaces
local internal face size
Definition: DAIndex.H:98
Foam::DAUtility::myFindCell
static label myFindCell(const primitiveMesh &mesh, const point &point)
Definition: DAUtility.C:803
Foam::DAFunction::functionDict_
dictionary functionDict_
dictionary containing the information for the objective function
Definition: DAFunction.H:64
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