DAFvSourceActuatorDisk.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAFvSourceActuatorDisk.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAFvSourceActuatorDisk, 0);
16 addToRunTimeSelectionTable(DAFvSource, DAFvSourceActuatorDisk, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAFvSourceActuatorDisk::DAFvSourceActuatorDisk(
20  const word modelType,
21  const fvMesh& mesh,
22  const DAOption& daOption,
23  const DAModel& daModel,
24  const DAIndex& daIndex)
25  : DAFvSource(modelType, mesh, daOption, daModel, daIndex)
26 {
27  this->initFvSourcePars();
28 
29  this->calcFvSourceCellIndices(fvSourceCellIndices_);
30 
31  printInterval_ = daOption.getOption<label>("printInterval");
32 }
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
36 {
37  /*
38  Description:
39  Compute the actuator disk source term.
40  We follow: Hoekstra, A RANS-based analysis tool for ducted propeller systems
41  in open water condition, International Shipbuilding Progress
42  source = rStar * sqrt(1.0 - rStar) * scale
43  where rStar is the normalized radial location and scale is used to
44  make sure the integral force equals the desired total thrust
45 
46  NOTE: rotDir = right means propeller rotates clockwise viewed from
47  the tail of the aircraft looking forward
48 
49  There are two options to assign the source term:
50  1. cylinderAnnulusToCell. Users prescribe a cylinderAnnulus and the fvSource will be
51  added to all the cells inside this cylinderAnnulus
52  2. cylinderAnnulusSmooth. Users prescribe the cylinderAnnulus and the Gaussian function
53  will be used to smoothly assign fvSource term. This allows us to use the actuatorDisk
54  parameters as design variables and move the actuator during optimization
55 
56  Example:
57  An example of the fvSource in pyOptions in pyDAFoam can be
58  defOptions =
59  {
60  "fvSource"
61  {
62  "disk1"
63  {
64  "type": "actuatorDisk",
65  "source": "cylinderAnnulusToCell",
66  "p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
67  "p2": [0.5, 0.3, 0.5], # p2-p1 should be streamwise
68  "innerRadius": 0.1,
69  "outerRadius": 0.8,
70  "rotDir": "left",
71  "scale": 1.0,
72  "POD": 0.7
73  },
74  "disk2"
75  {
76  "type": "actuatorDisk",
77  "source": "cylinderAnnulusSmooth",
78  "center": [0.0, 0.0, 0.0],
79  "direction": [1.0, 0.0, 0.0],
80  "innerRadius": 0.1,
81  "outerRadius": 0.8,
82  "rotDir": "right",
83  "scale": 0.1,
84  "POD": 1.0,
85  "eps": 0.05 # eps should be of cell size
86  "expM": 1.0,
87  "expN": 0.5,
88  }
89  }
90  }
91  */
92 
93  forAll(fvSource, idxI)
94  {
95  fvSource[idxI] = vector::zero;
96  }
97 
98  const dictionary& allOptions = daOption_.getAllOptions();
99 
100  dictionary fvSourceSubDict = allOptions.subDict("fvSource");
101 
102  forAll(fvSourceSubDict.toc(), idxI)
103  {
104  word diskName = fvSourceSubDict.toc()[idxI];
105  dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
106  word source = diskSubDict.getWord("source");
107 
108  if (source == "cylinderAnnulusToCell")
109  {
110  // loop over all the cell indices for all actuator disks
111 
112  // now read in all parameters for this actuator disk
113  scalarList point1;
114  scalarList point2;
115  diskSubDict.readEntry<scalarList>("p1", point1);
116  diskSubDict.readEntry<scalarList>("p2", point2);
117  vector p1 = {point1[0], point1[1], point1[2]};
118  vector p2 = {point2[0], point2[1], point2[2]};
119  vector diskCenter = (p1 + p2) / 2.0;
120  vector diskDir = p2 - p1; // NOTE: p2 - p1 should be streamwise
121  vector diskDirNorm = diskDir / mag(diskDir);
122  scalar outerRadius = diskSubDict.getScalar("outerRadius");
123  scalar innerRadius = diskSubDict.getScalar("innerRadius");
124  word rotDir = diskSubDict.getWord("rotDir");
125  scalar scale = diskSubDict.getScalar("scale");
126  scalar POD = diskSubDict.getScalar("POD");
127 
128  // loop over all cell indices for this disk and computer the source term
129  scalar thrustSourceSum = 0.0;
130  scalar torqueSourceSum = 0.0;
131  forAll(fvSourceCellIndices_[diskName], idxJ)
132  {
133  // cell index
134  label cellI = fvSourceCellIndices_[diskName][idxJ];
135 
136  // the cell center coordinates of this cellI
137  vector cellC = mesh_.C()[cellI];
138  // cell center to disk center vector
139  vector cellC2AVec = cellC - diskCenter;
140  // tmp tensor for calculating the axial/radial components of cellC2AVec
141  tensor cellC2AVecE(tensor::zero);
142  cellC2AVecE.xx() = cellC2AVec.x();
143  cellC2AVecE.yy() = cellC2AVec.y();
144  cellC2AVecE.zz() = cellC2AVec.z();
145 
146  // now we need to decompose cellC2AVec into axial and radial components
147  // the axial component of cellC2AVec vector
148  vector cellC2AVecA = cellC2AVecE & diskDirNorm;
149  // the radial component of cellC2AVec vector
150  vector cellC2AVecR = cellC2AVec - cellC2AVecA;
151 
152  // now we can use the cross product to compute the tangential
153  // (circ) direction of cellI
154  vector cellC2AVecC(vector::zero);
155  if (rotDir == "left")
156  {
157  // propeller rotates counter-clockwise viewed from the tail of the aircraft looking forward
158  cellC2AVecC = cellC2AVecR ^ diskDirNorm; // circ
159  }
160  else if (rotDir == "right")
161  {
162  // propeller rotates clockwise viewed from the tail of the aircraft looking forward
163  cellC2AVecC = diskDirNorm ^ cellC2AVecR; // circ
164  }
165  else
166  {
167  FatalErrorIn(" ") << "rotDir not valid" << abort(FatalError);
168  }
169 
170  // the magnitude of radial component of cellC2AVecR
171  scalar cellC2AVecRLen = mag(cellC2AVecR);
172  // the magnitude of tangential component of cellC2AVecR
173  scalar cellC2AVecCLen = mag(cellC2AVecC);
174  // the normalized cellC2AVecC (tangential) vector
175  vector cellC2AVecCNorm = cellC2AVecC / cellC2AVecCLen;
176 
177  // now we can use Hoekstra's formulation to compute source
178  scalar rPrime = cellC2AVecRLen / outerRadius;
179  scalar rPrimeHub = innerRadius / outerRadius;
180  // rStar is normalized radial location
181  scalar rStar = (rPrime - rPrimeHub) / (1.0 - rPrimeHub);
182  // axial force, NOTE: user need to prescribe "scale" such that the integrated
183  // axial force matches the desired thrust
184  scalar fAxial = rStar * sqrt(1.0 - rStar) * scale;
185  // we use Hoekstra's method to calculate the fCirc based on fAxial
186  scalar fCirc = fAxial * POD / constant::mathematical::pi / rPrime;
187  vector sourceVec = (fAxial * diskDirNorm + fCirc * cellC2AVecCNorm);
188  // the source is the force normalized by the cell volume
189  fvSource[cellI] += sourceVec;
190  thrustSourceSum += fAxial * mesh_.V()[cellI];
191  torqueSourceSum += fCirc * mesh_.V()[cellI];
192  }
193 
194  reduce(thrustSourceSum, sumOp<scalar>());
195  reduce(torqueSourceSum, sumOp<scalar>());
196 
197 #ifndef CODI_ADR
198  if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
199  {
200  Info << "ThrustCoeff Source Term for " << diskName << ": " << thrustSourceSum << endl;
201  Info << "TorqueCoeff Source Term for " << diskName << ": " << torqueSourceSum << endl;
202  }
203 #endif
204  }
205  else if (source == "cylinderAnnulusSmooth")
206  {
207  DAGlobalVar& globalVar =
208  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
209  HashTable<List<scalar>>& actuatorDiskPars = globalVar.actuatorDiskPars;
210 
211  vector center = {actuatorDiskPars[diskName][0], actuatorDiskPars[diskName][1], actuatorDiskPars[diskName][2]};
212  vector dirNorm = {actuatorDiskPars[diskName][3], actuatorDiskPars[diskName][4], actuatorDiskPars[diskName][5]};
213  dirNorm = dirNorm / mag(dirNorm);
214  scalar innerRadius = actuatorDiskPars[diskName][6];
215  scalar outerRadius = actuatorDiskPars[diskName][7];
216  word rotDir = diskSubDict.getWord("rotDir");
217  // we will calculate or read scale later
218  scalar scale;
219  scalar POD = actuatorDiskPars[diskName][9];
220  scalar eps = diskSubDict.getScalar("eps");
221  scalar expM = actuatorDiskPars[diskName][10];
222  scalar expN = actuatorDiskPars[diskName][11];
223  // Now we need to compute normalized eps in the radial direction, i.e. epsRStar this is because
224  // we need to smooth the radial distribution of the thrust, here the radial location is
225  // normalized as rStar = (r - rInner) / (rOuter - rInner), so to make epsRStar consistent with this
226  // we need to normalize eps with the demoninator of rStar, i.e. Outer - rInner
227  scalar epsRStar = eps / (outerRadius - innerRadius);
228  scalar rStarMin = epsRStar;
229  scalar rStarMax = 1.0 - epsRStar;
230  scalar fRMin = pow(rStarMin, expM) * pow(1.0 - rStarMin, expN);
231  scalar fRMax = pow(rStarMax, expM) * pow(1.0 - rStarMax, expN);
232 
233  label adjustThrust = diskSubDict.getLabel("adjustThrust");
234  // if adjustThrust = False, we just read "scale" from daOption
235  // if we want to adjust thrust, we calculate scale, instead of reading from daOption
236  // to calculate the scale, we just compute the fAxial with scale = 1, then we find
237  // the correct scale = targetThrust / thrust_with_scale_1
238  if (adjustThrust)
239  {
240  scale = 1.0;
241  scalar tmpThrustSumAll = 0.0;
242  forAll(mesh_.cells(), cellI)
243  {
244  // the cell center coordinates of this cellI
245  vector cellC = mesh_.C()[cellI];
246  // cell center to disk center vector
247  vector cellC2AVec = cellC - center;
248  // tmp tensor for calculating the axial/radial components of cellC2AVec
249  tensor cellC2AVecE(tensor::zero);
250  cellC2AVecE.xx() = cellC2AVec.x();
251  cellC2AVecE.yy() = cellC2AVec.y();
252  cellC2AVecE.zz() = cellC2AVec.z();
253 
254  // now we need to decompose cellC2AVec into axial and radial components
255  // the axial component of cellC2AVec vector
256  vector cellC2AVecA = cellC2AVecE & dirNorm;
257  // the radial component of cellC2AVec vector
258  vector cellC2AVecR = cellC2AVec - cellC2AVecA;
259 
260  // the magnitude of radial component of cellC2AVecR
261  scalar cellC2AVecRLen = mag(cellC2AVecR);
262  // the magnitude of axial component of cellC2AVecR
263  scalar cellC2AVecALen = mag(cellC2AVecA);
264 
265  // now we can use the smoothed formulation to compute source
266  scalar rPrime = cellC2AVecRLen / outerRadius;
267  scalar rPrimeHub = innerRadius / outerRadius;
268  // rStar is normalized radial location
269  scalar rStar = (rPrime - rPrimeHub) / (1.0 - rPrimeHub);
270 
271  scalar fAxial = 0.0;
272 
273  scalar dA2 = cellC2AVecALen * cellC2AVecALen;
274 
275  if (rStar < rStarMin)
276  {
277  scalar dR2 = (rStar - rStarMin) * (rStar - rStarMin);
278  scalar fR = fRMin * exp(-dR2 / epsRStar / epsRStar) * scale;
279  fAxial = fR * exp(-dA2 / eps / eps);
280  }
281  else if (rStar >= rStarMin && rStar <= rStarMax)
282  {
283  scalar fR = pow(rStar, expM) * pow(1.0 - rStar, expN) * scale;
284  fAxial = fR * exp(-dA2 / eps / eps);
285  }
286  else
287  {
288  scalar dR2 = (rStar - rStarMax) * (rStar - rStarMax);
289  scalar fR = fRMax * exp(-dR2 / epsRStar / epsRStar) * scale;
290  fAxial = fR * exp(-dA2 / eps / eps);
291  }
292 
293  tmpThrustSumAll += fAxial * mesh_.V()[cellI];
294  }
295  reduce(tmpThrustSumAll, sumOp<scalar>());
296  scalar targetThrust = actuatorDiskPars[diskName][12];
297  scale = targetThrust / tmpThrustSumAll;
298  }
299  else
300  {
301  scale = actuatorDiskPars[diskName][8];
302  }
303 
304  // now we have the correct scale, repeat the loop to assign fvSource
305  scalar thrustSourceSum = 0.0;
306  scalar torqueSourceSum = 0.0;
307  forAll(mesh_.cells(), cellI)
308  {
309  // the cell center coordinates of this cellI
310  vector cellC = mesh_.C()[cellI];
311  // cell center to disk center vector
312  vector cellC2AVec = cellC - center;
313  // tmp tensor for calculating the axial/radial components of cellC2AVec
314  tensor cellC2AVecE(tensor::zero);
315  cellC2AVecE.xx() = cellC2AVec.x();
316  cellC2AVecE.yy() = cellC2AVec.y();
317  cellC2AVecE.zz() = cellC2AVec.z();
318 
319  // now we need to decompose cellC2AVec into axial and radial components
320  // the axial component of cellC2AVec vector
321  vector cellC2AVecA = cellC2AVecE & dirNorm;
322  // the radial component of cellC2AVec vector
323  vector cellC2AVecR = cellC2AVec - cellC2AVecA;
324 
325  // now we can use the cross product to compute the tangential
326  // (circ) direction of cellI
327  vector cellC2AVecC(vector::zero);
328  if (rotDir == "left")
329  {
330  // propeller rotates counter-clockwise viewed from the tail of the aircraft looking forward
331  cellC2AVecC = cellC2AVecR ^ dirNorm; // circ
332  }
333  else if (rotDir == "right")
334  {
335  // propeller rotates clockwise viewed from the tail of the aircraft looking forward
336  cellC2AVecC = dirNorm ^ cellC2AVecR; // circ
337  }
338  else
339  {
340  FatalErrorIn(" ") << "rotDir not valid" << abort(FatalError);
341  }
342 
343  // the magnitude of radial component of cellC2AVecR
344  scalar cellC2AVecRLen = mag(cellC2AVecR);
345  // the magnitude of tangential component of cellC2AVecR
346  scalar cellC2AVecCLen = mag(cellC2AVecC);
347  // the magnitude of axial component of cellC2AVecR
348  scalar cellC2AVecALen = mag(cellC2AVecA);
349  // the normalized cellC2AVecC (tangential) vector
350  vector cellC2AVecCNorm = cellC2AVecC / cellC2AVecCLen;
351 
352  // now we can use the smoothed formulation to compute source
353  scalar rPrime = cellC2AVecRLen / outerRadius;
354  scalar rPrimeHub = innerRadius / outerRadius;
355  // rStar is normalized radial location
356  scalar rStar = (rPrime - rPrimeHub) / (1.0 - rPrimeHub);
357 
358  scalar fAxial = 0.0;
359  scalar dA2 = cellC2AVecALen * cellC2AVecALen;
360 
361  if (rStar < rStarMin)
362  {
363  scalar dR2 = (rStar - rStarMin) * (rStar - rStarMin);
364  scalar fR = fRMin * exp(-dR2 / epsRStar / epsRStar) * scale;
365  fAxial = fR * exp(-dA2 / eps / eps);
366  }
367  else if (rStar >= rStarMin && rStar <= rStarMax)
368  {
369  scalar fR = pow(rStar, expM) * pow(1.0 - rStar, expN) * scale;
370  fAxial = fR * exp(-dA2 / eps / eps);
371  }
372  else
373  {
374  scalar dR2 = (rStar - rStarMax) * (rStar - rStarMax);
375  scalar fR = fRMax * exp(-dR2 / epsRStar / epsRStar) * scale;
376  fAxial = fR * exp(-dA2 / eps / eps);
377  }
378  // we use Hoekstra's method to calculate the fCirc based on fAxial
379  // here we add 0.01*eps/outerRadius to avoid diving a zero rPrime
380  // this might happen if a cell center is very close to actuator center
381  scalar fCirc = fAxial * POD / constant::mathematical::pi / (rPrime + 0.01 * eps / outerRadius);
382 
383  vector sourceVec = (fAxial * dirNorm + fCirc * cellC2AVecCNorm);
384  // the source is the force normalized by the cell volume
385  fvSource[cellI] += sourceVec;
386  thrustSourceSum += fAxial * mesh_.V()[cellI];
387  torqueSourceSum += fCirc * mesh_.V()[cellI];
388  }
389 
390  reduce(thrustSourceSum, sumOp<scalar>());
391  reduce(torqueSourceSum, sumOp<scalar>());
392 
393 #ifndef CODI_ADR
394  if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
395  {
396  Info << "ThrustCoeff Source Term for " << diskName << ": " << thrustSourceSum << endl;
397  Info << "TorqueCoeff Source Term for " << diskName << ": " << torqueSourceSum << endl;
398  if (adjustThrust)
399  {
400  Info << "Dynamically adjusted scale for " << diskName << ": " << scale << endl;
401  }
402  if (daOption_.getOption<label>("debug"))
403  {
404  Info << "adjustThrust for " << diskName << ": " << adjustThrust << endl;
405  Info << "center for " << diskName << ": " << center << endl;
406  Info << "innerRadius for " << diskName << ": " << innerRadius << endl;
407  Info << "outerRadius for " << diskName << ": " << outerRadius << endl;
408  Info << "scale for " << diskName << ": " << scale << endl;
409  Info << "POD for " << diskName << ": " << POD << endl;
410  Info << "eps for " << diskName << ": " << eps << endl;
411  Info << "expM for " << diskName << ": " << expM << endl;
412  Info << "expN for " << diskName << ": " << expN << endl;
413  Info << "epsRStar for " << diskName << ": " << epsRStar << endl;
414  Info << "rStarMin for " << diskName << ": " << rStarMin << endl;
415  Info << "rStarMax for " << diskName << ": " << rStarMax << endl;
416  }
417  }
418 #endif
419  }
420  else
421  {
422  FatalErrorIn("calcFvSourceCells") << "source: " << source << " not supported!"
423  << "Options are: cylinderAnnulusToCell and cylinderAnnulusSmooth!"
424  << abort(FatalError);
425  }
426  }
427 
428  fvSource.correctBoundaryConditions();
429 }
430 
431 void DAFvSourceActuatorDisk::initFvSourcePars()
432 {
433  /*
434  Description:
435  Initialize the values for all types of fvSource in DAGlobalVar, including
436  actuatorDiskPars, heatSourcePars, etc
437  */
438 
439  // now we need to initialize actuatorDiskPars
440  dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
441 
442  forAll(fvSourceSubDict.toc(), idxI)
443  {
444  word diskName = fvSourceSubDict.toc()[idxI];
445  // sub dictionary with all parameters for this disk
446  dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
447  word type = diskSubDict.getWord("type");
448  word source = diskSubDict.getWord("source");
449 
450  if (type == "actuatorDisk" && source == "cylinderAnnulusSmooth")
451  {
452 
453  // now read in all parameters for this actuator disk
454  scalarList centerList;
455  diskSubDict.readEntry<scalarList>("center", centerList);
456 
457  scalarList dirList;
458  diskSubDict.readEntry<scalarList>("direction", dirList);
459 
460  // we have 13 design variables for each disk
461  scalarList dvList(13);
462  dvList[0] = centerList[0];
463  dvList[1] = centerList[1];
464  dvList[2] = centerList[2];
465  dvList[3] = dirList[0];
466  dvList[4] = dirList[1];
467  dvList[5] = dirList[2];
468  dvList[6] = diskSubDict.getScalar("innerRadius");
469  dvList[7] = diskSubDict.getScalar("outerRadius");
470  dvList[8] = diskSubDict.getScalar("scale");
471  dvList[9] = diskSubDict.getScalar("POD");
472  dvList[10] = diskSubDict.getScalar("expM");
473  dvList[11] = diskSubDict.getScalar("expN");
474  dvList[12] = diskSubDict.getScalar("targetThrust");
475 
476  // set actuatorDiskPars
477  DAGlobalVar& globalVar =
478  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
479  globalVar.actuatorDiskPars.set(diskName, dvList);
480  }
481  }
482 }
483 
484 void DAFvSourceActuatorDisk::calcFvSourceCellIndices(HashTable<labelList>& fvSourceCellIndices)
485 {
486  /*
487  Description:
488  Calculate the lists of cell indices that are within the actuator disk space
489  NOTE: we support multiple actuator disks.
490 
491  Output:
492  fvSourceCellIndices: Hash table that contains the lists of cell indices that
493  are within the actuator disk space. The hash key is the name of the disk. We support
494  multiple disks. An example of fvSourceCellIndices can be:
495 
496  fvSourceCellIndices=
497  {
498  "disk1": {1,2,3,4,5},
499  "disk2": {6,7,8,9,10,11,12}
500  }
501 
502  Example:
503  An example of the fvSource in pyOptions in pyDAFoam can be
504  defOptions =
505  {
506  "fvSource"
507  {
508  "disk1"
509  {
510  "type": "actuatorDisk",
511  "source": cylinderAnnulusToCell,
512  "p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
513  "p2": [0.5, 0.3, 0.5],
514  "innerRadius": 0.1,
515  "outerRadius": 0.8,
516  "rotDir": "left",
517  "scale": 1.0,
518  "POD": 0.7
519  },
520  "disk2"
521  {
522  "type": "actuatorDisk",
523  "source": cylinderAnnulusToCell,
524  "p1": [0.0, 0.0, 0.1],
525  "p2": [0.0, 0.0, 0.5],
526  "innerRadius": 0.1,
527  "outerRadius": 0.8,
528  "rotDir": "right",
529  "scale": 0.1,
530  "POD": 1.0
531  }
532  }
533  }
534  */
535 
536  const dictionary& allOptions = daOption_.getAllOptions();
537 
538  dictionary fvSourceSubDict = allOptions.subDict("fvSource");
539 
540  if (fvSourceSubDict.toc().size() == 0)
541  {
542  FatalErrorIn("calcSourceCells") << "actuatorDisk is selected as fvSource "
543  << " but the options are empty!"
544  << abort(FatalError);
545  }
546 
547  forAll(fvSourceSubDict.toc(), idxI)
548  {
549  word diskName = fvSourceSubDict.toc()[idxI];
550 
551  fvSourceCellIndices.set(diskName, {});
552 
553  dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
554  word sourceType = diskSubDict.getWord("source");
555  // all avaiable source type are in src/meshTools/sets/cellSources
556  // Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet
557  if (sourceType == "cylinderAnnulusToCell")
558  {
559  // create a topoSet
560  autoPtr<topoSet> currentSet(
561  topoSet::New(
562  "cellSet",
563  mesh_,
564  "set0",
565  IOobject::NO_READ));
566  // we need to change the min and max because they need to
567  // be of type point; however, we can't parse point type
568  // in pyDict, we need to change them here.
569 
570  scalarList point1;
571  scalarList point2;
572  diskSubDict.readEntry<scalarList>("p1", point1);
573  diskSubDict.readEntry<scalarList>("p2", point2);
574 
575  point p1;
576  point p2;
577  p1[0] = point1[0];
578  p1[1] = point1[1];
579  p1[2] = point1[2];
580  p2[0] = point2[0];
581  p2[1] = point2[1];
582  p2[2] = point2[2];
583 
584  scalar outerRadius = diskSubDict.getScalar("outerRadius");
585  scalar innerRadius = diskSubDict.getScalar("innerRadius");
586 
587  dictionary tmpDict;
588  tmpDict.set("p1", p1);
589  tmpDict.set("p2", p2);
590  tmpDict.set("innerRadius", innerRadius);
591  tmpDict.set("outerRadius", outerRadius);
592 
593  // create the source
594  autoPtr<topoSetSource> sourceSet(
595  topoSetSource::New(sourceType, mesh_, tmpDict));
596 
597  // add the sourceSet to topoSet
598  sourceSet().applyToSet(topoSetSource::NEW, currentSet());
599  // get the face index from currentSet, we need to use
600  // this special for loop
601  for (const label i : currentSet())
602  {
603  fvSourceCellIndices[diskName].append(i);
604  }
605  }
606  else if (sourceType == "cylinderAnnulusSmooth")
607  {
608  // do nothing, no need to compute the cell indices since
609  // we are using Gaussian function to compute a smooth
610  // distribution of fvSource term
611  }
612  else
613  {
614  FatalErrorIn("calcFvSourceCells") << "source: " << sourceType << " not supported!"
615  << "Options are: cylinderAnnulusToCell and cylinderAnnulusSmooth!"
616  << abort(FatalError);
617  }
618  }
619 
620  if (daOption_.getOption<label>("debug"))
621  {
622  Info << "fvSourceCellIndices " << fvSourceCellIndices << endl;
623  }
624 }
625 
626 } // End namespace Foam
627 
628 // ************************************************************************* //
Foam::DAFvSource
Definition: DAFvSource.H:34
Foam::DAFvSourceActuatorDisk::printInterval_
label printInterval_
print interval for primal and adjoint
Definition: DAFvSourceActuatorDisk.H:39
Foam::DAGlobalVar::actuatorDiskPars
HashTable< List< scalar > > actuatorDiskPars
the list of parameters for all the actuator disks
Definition: DAGlobalVar.H:65
DAFvSourceActuatorDisk.H
Foam::DAFvSource::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAFvSource.H:53
Foam::DAOption
Definition: DAOption.H:29
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
Foam::DAOption::getOption
classType getOption(const word key) const
get an option from subDict and key
Definition: DAOption.H:92
fvSource
volScalarField & fvSource
Definition: createRefsHeatTransfer.H:7
Foam::DAFvSourceActuatorDisk::fvSourceCellIndices_
HashTable< labelList > fvSourceCellIndices_
HashTable that contains lists of cell indices that are within the actuator disk space.
Definition: DAFvSourceActuatorDisk.H:33
Foam::DAFvSourceActuatorDisk::initFvSourcePars
virtual void initFvSourcePars()
Initialize the values for all types of fvSource in DAGlobalVar, including actuatorDiskPars,...
Definition: DAFvSourceActuatorDisk.C:431
Foam::DAOption::getAllOptions
const dictionary & getAllOptions() const
return a reference of allOptions_ dictionary
Definition: DAOption.H:56
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAFvSourceActuatorDisk::DAFvSourceActuatorDisk
DAFvSourceActuatorDisk(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAFvSourceActuatorDisk.C:19
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAModel
Definition: DAModel.H:57
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam
Definition: checkGeometry.C:32
allOptions
const dictionary & allOptions
Definition: createRefsPimpleDyM.H:14
p1
p1
Definition: p1EqnIrkPimple.H:24
source
pseudoPEqn source()
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAFvSourceActuatorDisk::calcFvSource
virtual void calcFvSource(volVectorField &fvSource)
compute the FvSource term
Definition: DAFvSourceActuatorDisk.C:35
Foam::DAFvSource::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAFvSource.H:50
Foam::DAFvSourceActuatorDisk::calcFvSourceCellIndices
void calcFvSourceCellIndices(HashTable< labelList > &fvSourceCellIndices)
calculate DAFvSourceActuatorDisk::fvSourceCellIndices_
Definition: DAFvSourceActuatorDisk.C:484
Foam::DAGlobalVar
Definition: DAGlobalVar.H:26
p2
p2
Definition: p2EqnIrkPimple.H:24

Generated by   doxygen 1.8.17