DAFvSourceHeatSource.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAFvSourceHeatSource.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAFvSourceHeatSource, 0);
16 addToRunTimeSelectionTable(DAFvSource, DAFvSourceHeatSource, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAFvSourceHeatSource::DAFvSourceHeatSource(
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  printInterval_ = daOption.getOption<label>("printInterval");
30 
31  const dictionary& allOptions = daOption_.getAllOptions();
32 
33  dictionary fvSourceSubDict = allOptions.subDict("fvSource");
34 
35  forAll(fvSourceSubDict.toc(), idxI)
36  {
37  word sourceName = fvSourceSubDict.toc()[idxI];
38  dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
39  word sourceType = sourceSubDict.getWord("source");
40 
41  if (sourceType == "cylinderToCell")
42  {
43  scalarList p1List;
44  sourceSubDict.readEntry<scalarList>("p1", p1List);
45  cylinderP1_.set(sourceName, p1List);
46  scalarList p2List;
47  sourceSubDict.readEntry<scalarList>("p2", p2List);
48  cylinderP2_.set(sourceName, p2List);
49 
50  cylinderRadius_.set(sourceName, sourceSubDict.getScalar("radius"));
51  power_.set(sourceName, sourceSubDict.getScalar("power"));
52 
53  fvSourceCellIndices_.set(sourceName, {});
54  // all available source type are in src/meshTools/sets/cellSources
55  // Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet
56 
57  // create a topoSet
58  autoPtr<topoSet> currentSet(
59  topoSet::New(
60  "cellSet",
61  mesh_,
62  "set0",
63  IOobject::NO_READ));
64  // we need to change the min and max because they need to
65  // be of type point; however, we can't parse point type
66  // in pyDict, we need to change them here.
67 
68  point p1;
69  point p2;
70  p1[0] = cylinderP1_[sourceName][0];
71  p1[1] = cylinderP1_[sourceName][1];
72  p1[2] = cylinderP1_[sourceName][2];
73  p2[0] = cylinderP2_[sourceName][0];
74  p2[1] = cylinderP2_[sourceName][1];
75  p2[2] = cylinderP2_[sourceName][2];
76 
77  dictionary tmpDict;
78  tmpDict.set("p1", p1);
79  tmpDict.set("p2", p2);
80  tmpDict.set("radius", cylinderRadius_[sourceName]);
81 
82  // create the source
83  autoPtr<topoSetSource> sourceSet(
84  topoSetSource::New("cylinderToCell", mesh_, tmpDict));
85 
86  // add the sourceSet to topoSet
87  sourceSet().applyToSet(topoSetSource::NEW, currentSet());
88  // get the face index from currentSet, we need to use
89  // this special for loop
90  for (const label i : currentSet())
91  {
92  fvSourceCellIndices_[sourceName].append(i);
93  }
94 
95  if (daOption_.getOption<label>("debug"))
96  {
97  Info << "fvSourceCellIndices " << fvSourceCellIndices_ << endl;
98  }
99  }
100  else if (sourceType == "cylinderSmooth")
101  {
102  DAGlobalVar& globalVar =
103  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
104  HashTable<List<scalar>>& heatSourcePars = globalVar.heatSourcePars;
105 
106  // eps is a smoothing parameter, it should be the local mesh cell size in meters
107  // near the cylinder region
108  cylinderEps_.set(sourceName, sourceSubDict.getScalar("eps"));
109 
110  snapCenter2Cell_.set(sourceName, sourceSubDict.lookupOrDefault<label>("snapCenter2Cell", 0));
111  if (snapCenter2Cell_[sourceName])
112  {
113  point centerPoint = {heatSourcePars[sourceName][0], heatSourcePars[sourceName][1], heatSourcePars[sourceName][2]};
114 
115  // NOTE: we need to call a self-defined findCell func to make it work correctly in ADR
116  label myCellI = DAUtility::myFindCell(mesh_, centerPoint);
117 
118  snappedCenterCellI_.set(sourceName, myCellI);
119  label foundCellI = 0;
120  if (snappedCenterCellI_[sourceName] >= 0)
121  {
122  foundCellI = 1;
123  //Pout << "snap source " << sourceName << " to center " << mesh_.C()[snappedCenterCellI_[sourceName]] << endl;
124  }
125  reduce(foundCellI, sumOp<label>());
126  if (foundCellI != 1)
127  {
128  FatalErrorIn(" ") << "There should be only one cell found globally while "
129  << foundCellI << " was returned."
130  << " Please adjust center such that it is located completely"
131  << " within a cell in the mesh domain. The center should not "
132  << " be outside of the mesh domain or on a mesh face "
133  << abort(FatalError);
134  }
135 
136  vector snappedCenter = vector::zero;
137  this->findGlobalSnappedCenter(snappedCenterCellI_[sourceName], snappedCenter);
138  Info << "heat source " << sourceName << " snap to center " << snappedCenter << endl;
139  }
140  }
141  else
142  {
143  FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
144  << "Options are: cylinderAnnulusToCell and cylinderSmooth!"
145  << abort(FatalError);
146  }
147  }
148 }
149 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
150 
151 void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
152 {
153  /*
154  Description:
155  Calculate the heat source and assign them to fvSource
156 
157  Example:
158  An example of the fvSource in pyOptions in pyDAFoam can be
159  defOptions =
160  {
161  "fvSource"
162  {
163  "source1"
164  {
165  "type": "heatSource",
166  "source": "cylinderToCell",
167  "p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
168  "p2": [0.5, 0.3, 0.5], # p2-p1 should be the axis of the cylinder
169  "radius": 0.8,
170  "power": 101.0, # here we should prescribe the power in W
171  },
172  "source2"
173  {
174  "type": "heatSource",
175  "source": "cylinder",
176  "p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
177  "p2": [0.5, 0.3, 0.5], # p2-p1 should be the axis of the cylinder
178  "radius": 1.5,
179  "power": 10.5, # here we should prescribe the power in W
180  }
181  }
182  }
183  */
184 
185  forAll(fvSource, idxI)
186  {
187  fvSource[idxI] = 0.0;
188  }
189 
190  const dictionary& allOptions = daOption_.getAllOptions();
191 
192  dictionary fvSourceSubDict = allOptions.subDict("fvSource");
193 
194  forAll(fvSourceSubDict.toc(), idxI)
195  {
196  word sourceName = fvSourceSubDict.toc()[idxI];
197  dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
198  word sourceType = sourceSubDict.getWord("source");
199  // NOTE: here power should be in W. We will evenly divide the power by the
200  // total volume of the source
201  if (sourceType == "cylinderToCell")
202  {
203 
204  scalar power = power_[sourceName];
205 
206  // loop over all cell indices for this source and assign the source term
207  // ----- first loop, calculate the total volume
208  scalar totalV = 0.0;
209  forAll(fvSourceCellIndices_[sourceName], idxJ)
210  {
211  label cellI = fvSourceCellIndices_[sourceName][idxJ];
212  scalar cellV = mesh_.V()[cellI];
213  totalV += cellV;
214  }
215  reduce(totalV, sumOp<scalar>());
216 
217  // ------- second loop, assign power
218  scalar sourceTotal = 0.0;
219  forAll(fvSourceCellIndices_[sourceName], idxJ)
220  {
221  // cell index
222  label cellI = fvSourceCellIndices_[sourceName][idxJ];
223 
224  fvSource[cellI] += power / totalV;
225  sourceTotal += fvSource[cellI] * mesh_.V()[cellI];
226  }
227 
228  reduce(sourceTotal, sumOp<scalar>());
229 
230 #ifndef CODI_ADR
231  if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
232  {
233  Info << "Total volume for " << sourceName << ": " << totalV << " m^3" << endl;
234  Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
235  }
236 #endif
237  }
238  else if (sourceType == "cylinderSmooth")
239  {
240  DAGlobalVar& globalVar =
241  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
242  HashTable<List<scalar>>& heatSourcePars = globalVar.heatSourcePars;
243 
244  vector cylinderCenter =
245  {heatSourcePars[sourceName][0], heatSourcePars[sourceName][1], heatSourcePars[sourceName][2]};
246 
247  if (snapCenter2Cell_[sourceName])
248  {
249  this->findGlobalSnappedCenter(snappedCenterCellI_[sourceName], cylinderCenter);
250  }
251  vector cylinderDir =
252  {heatSourcePars[sourceName][3], heatSourcePars[sourceName][4], heatSourcePars[sourceName][5]};
253  scalar radius = heatSourcePars[sourceName][6];
254  scalar cylinderLen = heatSourcePars[sourceName][7];
255  scalar power = heatSourcePars[sourceName][8];
256  vector cylinderDirNorm = cylinderDir / mag(cylinderDir);
257  scalar eps = cylinderEps_[sourceName];
258 
259  // first loop, calculate the total volume
260  scalar volumeTotal = 0.0;
261  forAll(mesh_.cells(), cellI)
262  {
263  // the cell center coordinates of this cellI
264  vector cellC = mesh_.C()[cellI];
265  // cell center to disk center vector
266  vector cellC2AVec = cellC - cylinderCenter;
267  // tmp tensor for calculating the axial/radial components of cellC2AVec
268  tensor cellC2AVecE(tensor::zero);
269  cellC2AVecE.xx() = cellC2AVec.x();
270  cellC2AVecE.yy() = cellC2AVec.y();
271  cellC2AVecE.zz() = cellC2AVec.z();
272 
273  // now we need to decompose cellC2AVec into axial and radial components
274  // the axial component of cellC2AVec vector
275  vector cellC2AVecA = cellC2AVecE & cylinderDirNorm;
276  // the radial component of cellC2AVec vector
277  vector cellC2AVecR = cellC2AVec - cellC2AVecA;
278 
279  // the magnitude of radial component of cellC2AVecR
280  scalar cellC2AVecRLen = mag(cellC2AVecR);
281  // the magnitude of axial component of cellC2AVecR
282  scalar cellC2AVecALen = mag(cellC2AVecA);
283 
284  scalar axialSmoothC = 1.0;
285  scalar radialSmoothC = 1.0;
286 
287  scalar halfL = cylinderLen / 2.0;
288 
289  if (cellC2AVecALen > halfL)
290  {
291  scalar d2 = (cellC2AVecALen - halfL) * (cellC2AVecALen - halfL);
292  axialSmoothC = exp(-d2 / eps / eps);
293  }
294 
295  if (cellC2AVecRLen > radius)
296  {
297  scalar d2 = (cellC2AVecRLen - radius) * (cellC2AVecRLen - radius);
298  radialSmoothC = exp(-d2 / eps / eps);
299  }
300  volumeTotal += axialSmoothC * radialSmoothC * mesh_.V()[cellI];
301  }
302  reduce(volumeTotal, sumOp<scalar>());
303 
304  // second loop, calculate fvSource
305  scalar sourceTotal = 0.0;
306  forAll(mesh_.cells(), cellI)
307  {
308  // the cell center coordinates of this cellI
309  vector cellC = mesh_.C()[cellI];
310  // cell center to disk center vector
311  vector cellC2AVec = cellC - cylinderCenter;
312  // tmp tensor for calculating the axial/radial components of cellC2AVec
313  tensor cellC2AVecE(tensor::zero);
314  cellC2AVecE.xx() = cellC2AVec.x();
315  cellC2AVecE.yy() = cellC2AVec.y();
316  cellC2AVecE.zz() = cellC2AVec.z();
317 
318  // now we need to decompose cellC2AVec into axial and radial components
319  // the axial component of cellC2AVec vector
320  vector cellC2AVecA = cellC2AVecE & cylinderDirNorm;
321  // the radial component of cellC2AVec vector
322  vector cellC2AVecR = cellC2AVec - cellC2AVecA;
323 
324  // the magnitude of radial component of cellC2AVecR
325  scalar cellC2AVecRLen = mag(cellC2AVecR);
326  // the magnitude of axial component of cellC2AVecR
327  scalar cellC2AVecALen = mag(cellC2AVecA);
328 
329  scalar axialSmoothC = 1.0;
330  scalar radialSmoothC = 1.0;
331 
332  scalar halfL = cylinderLen / 2.0;
333 
334  if (cellC2AVecALen > halfL)
335  {
336  scalar d2 = (cellC2AVecALen - halfL) * (cellC2AVecALen - halfL);
337  axialSmoothC = exp(-d2 / eps / eps);
338  }
339 
340  if (cellC2AVecRLen > radius)
341  {
342  scalar d2 = (cellC2AVecRLen - radius) * (cellC2AVecRLen - radius);
343  radialSmoothC = exp(-d2 / eps / eps);
344  }
345 
346  fvSource[cellI] += power / volumeTotal * axialSmoothC * radialSmoothC;
347  sourceTotal += power / volumeTotal * axialSmoothC * radialSmoothC * mesh_.V()[cellI];
348  }
349 
350  reduce(sourceTotal, sumOp<scalar>());
351 
352 #ifndef CODI_ADR
353  if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
354  {
355  Info << "Total volume for " << sourceName << ": " << volumeTotal << " m^3" << endl;
356  Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
357  }
358 #endif
359  }
360  else
361  {
362  FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
363  << "Options are: cylinderAnnulusToCell and cylinderSmooth!"
364  << abort(FatalError);
365  }
366  }
367 
368  fvSource.correctBoundaryConditions();
369 }
370 
371 void DAFvSourceHeatSource::initFvSourcePars()
372 {
373  /*
374  Description:
375  Initialize the values for all types of fvSource in DAGlobalVar, including
376  actuatorDiskPars, heatSourcePars, etc
377  */
378 
379  // now we need to initialize actuatorDiskPars
380  dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
381 
382  forAll(fvSourceSubDict.toc(), idxI)
383  {
384  word diskName = fvSourceSubDict.toc()[idxI];
385  // sub dictionary with all parameters for this disk
386  dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
387  word type = diskSubDict.getWord("type");
388  word source = diskSubDict.getWord("source");
389 
390  if (type == "heatSource" && source == "cylinderSmooth")
391  {
392  // now read in all parameters for this actuator disk
393  scalarList centerList;
394  diskSubDict.readEntry<scalarList>("center", centerList);
395 
396  scalarList axisList;
397  diskSubDict.readEntry<scalarList>("axis", axisList);
398 
399  // we have 13 design variables for each disk
400  scalarList dvList(9);
401  dvList[0] = centerList[0];
402  dvList[1] = centerList[1];
403  dvList[2] = centerList[2];
404  dvList[3] = axisList[0];
405  dvList[4] = axisList[1];
406  dvList[5] = axisList[2];
407  dvList[6] = diskSubDict.getScalar("radius");
408  dvList[7] = diskSubDict.getScalar("length");
409  dvList[8] = diskSubDict.getScalar("power");
410 
411  // set actuatorDiskPars
412  DAGlobalVar& globalVar =
413  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
414  globalVar.heatSourcePars.set(diskName, dvList);
415  }
416  }
417 }
418 
419 } // End namespace Foam
420 
421 // ************************************************************************* //
DAFvSourceHeatSource.H
Foam::DAFvSource
Definition: DAFvSource.H:34
Foam::DAFvSourceHeatSource::snappedCenterCellI_
HashTable< label > snappedCenterCellI_
the cell index for the center if snapCenter2Cell_ = 1
Definition: DAFvSourceHeatSource.H:48
Foam::DAFvSourceHeatSource::DAFvSourceHeatSource
DAFvSourceHeatSource(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAFvSourceHeatSource.C:19
Foam::DAFvSourceHeatSource::printInterval_
label printInterval_
print interval for primal and adjoint
Definition: DAFvSourceHeatSource.H:42
Foam::DAFvSource::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAFvSource.H:53
Foam::DAOption
Definition: DAOption.H:29
Foam::DAFvSourceHeatSource::initFvSourcePars
virtual void initFvSourcePars()
Initialize the values for all types of fvSource in DAGlobalVar, including actuatorDiskPars,...
Definition: DAFvSourceHeatSource.C:371
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
Foam::DAFvSourceHeatSource::power_
HashTable< scalar > power_
Definition: DAFvSourceHeatSource.H:38
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::DAOption::getAllOptions
const dictionary & getAllOptions() const
return a reference of allOptions_ dictionary
Definition: DAOption.H:56
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAFvSourceHeatSource::fvSourceCellIndices_
HashTable< labelList > fvSourceCellIndices_
HashTable that contains lists of cell indices that are within the actuator disk space.
Definition: DAFvSourceHeatSource.H:33
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAGlobalVar::heatSourcePars
HashTable< List< scalar > > heatSourcePars
the list of parameters for all the actuator disks
Definition: DAGlobalVar.H:74
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::DAUtility::myFindCell
static label myFindCell(const primitiveMesh &mesh, const point &point)
Definition: DAUtility.C:803
Foam::DAFvSourceHeatSource::cylinderP2_
HashTable< scalarList > cylinderP2_
Definition: DAFvSourceHeatSource.H:36
Foam::DAFvSource::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAFvSource.H:50
Foam::DAFvSourceHeatSource::snapCenter2Cell_
HashTable< label > snapCenter2Cell_
whether to snap the center to a cell in the mesh if yes the center will move with the mesh
Definition: DAFvSourceHeatSource.H:45
Foam::DAFvSourceHeatSource::cylinderP1_
HashTable< scalarList > cylinderP1_
Definition: DAFvSourceHeatSource.H:35
Foam::DAFvSourceHeatSource::calcFvSource
virtual void calcFvSource(volScalarField &fvSource)
compute the FvSource term
Definition: DAFvSourceHeatSource.C:151
Foam::DAFvSource::findGlobalSnappedCenter
void findGlobalSnappedCenter(label snappedCenterCellI, vector &center)
Definition: DAFvSource.C:140
Foam::DAFvSourceHeatSource::cylinderEps_
HashTable< scalar > cylinderEps_
Definition: DAFvSourceHeatSource.H:39
Foam::DAGlobalVar
Definition: DAGlobalVar.H:26
Foam::DAFvSourceHeatSource::cylinderRadius_
HashTable< scalar > cylinderRadius_
Definition: DAFvSourceHeatSource.H:37
p2
p2
Definition: p2EqnIrkPimple.H:24

Generated by   doxygen 1.8.17