DAField.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAField.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
16 
17 DAField::DAField(
18  const fvMesh& mesh,
19  const DAOption& daOption,
20  const DAModel& daModel,
21  const DAIndex& daIndex)
22  : mesh_(mesh),
23  daOption_(daOption),
24  daModel_(daModel),
25  daIndex_(daIndex)
26 {
27  // initialize stateInfo_
28  word solverName = daOption.getOption<word>("solverName");
29  autoPtr<DAStateInfo> daStateInfo(DAStateInfo::New(solverName, mesh, daOption, daModel));
30  stateInfo_ = daStateInfo->getStateInfo();
31 
32  // check if we have special boundary conditions that need special treatment
33  this->checkSpecialBCs();
34 }
35 
36 void DAField::ofField2StateVec(Vec stateVec) const
37 {
38  /*
39  Description:
40  Assign values for the state variable vector based on the
41  latest OpenFOAM field values
42 
43  Input:
44  OpenFOAM field variables
45 
46  Output:
47  stateVec: state variable vector
48 
49  Example:
50  Image we have two state variables (p,T) and five cells, running on two CPU
51  processors, the proc0 owns two cells and the proc1 owns three cells,
52  then calling this function gives the state vector (state-by-state ordering):
53 
54  stateVec = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
55  0 1 2 3 | 4 5 6 7 8 9 <- global state vec index
56  ---- proc0 -----|--------- proc1 -------
57  */
58 
59  const objectRegistry& db = mesh_.thisDb();
60  PetscScalar* stateVecArray;
61  VecGetArray(stateVec, &stateVecArray);
62 
63  forAll(stateInfo_["volVectorStates"], idxI)
64  {
65  // lookup state from meshDb
66  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
67 
68  forAll(mesh_.cells(), cellI)
69  {
70  for (label comp = 0; comp < 3; comp++)
71  {
72  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
73  assignValueCheckAD(stateVecArray[localIdx], state[cellI][comp]);
74  }
75  }
76  }
77 
78  forAll(stateInfo_["volScalarStates"], idxI)
79  {
80  // lookup state from meshDb
81  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
82 
83  forAll(mesh_.cells(), cellI)
84  {
85  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
86  assignValueCheckAD(stateVecArray[localIdx], state[cellI]);
87  }
88  }
89 
90  forAll(stateInfo_["modelStates"], idxI)
91  {
92  // lookup state from meshDb
93  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
94 
95  forAll(mesh_.cells(), cellI)
96  {
97  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
98  assignValueCheckAD(stateVecArray[localIdx], state[cellI]);
99  }
100  }
101 
102  forAll(stateInfo_["surfaceScalarStates"], idxI)
103  {
104  // lookup state from meshDb
105  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
106 
107  forAll(mesh_.faces(), faceI)
108  {
109  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
110  if (faceI < daIndex_.nLocalInternalFaces)
111  {
112  assignValueCheckAD(stateVecArray[localIdx], state[faceI]);
113  }
114  else
115  {
116  label relIdx = faceI - daIndex_.nLocalInternalFaces;
117  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
118  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
119  assignValueCheckAD(stateVecArray[localIdx], state.boundaryField()[patchIdx][faceIdx]);
120  }
121  }
122  }
123  VecRestoreArray(stateVec, &stateVecArray);
124 }
125 
126 void DAField::state2OFField(const scalar* states) const
127 {
128  /*
129  Description:
130  Assign values OpenFOAM field values based on the state variable array
131 
132  Input:
133  states: state variable array
134 
135  Output:
136  OpenFoam field variables
137 
138  Example:
139  Image we have two state variables (p,T) and five cells, running on two CPU
140  processors, the proc0 owns two cells and the proc1 owns three cells,
141  then calling this function will assign the p, and T based on the the state
142  vector (state-by-state ordering):
143 
144  stateVec = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
145  0 1 2 3 | 4 5 6 7 8 9 <- global state vec index
146  ---- proc0 -----|--------- proc1 -------
147  */
148 
149  const objectRegistry& db = mesh_.thisDb();
150 
151  forAll(stateInfo_["volVectorStates"], idxI)
152  {
153  // lookup state from meshDb
154  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
155 
156  forAll(mesh_.cells(), cellI)
157  {
158  for (label comp = 0; comp < 3; comp++)
159  {
160  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
161  state[cellI][comp] = states[localIdx];
162  }
163  }
164  }
165 
166  forAll(stateInfo_["volScalarStates"], idxI)
167  {
168  // lookup state from meshDb
169  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
170 
171  forAll(mesh_.cells(), cellI)
172  {
173  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
174  state[cellI] = states[localIdx];
175  }
176  }
177 
178  forAll(stateInfo_["modelStates"], idxI)
179  {
180  // lookup state from meshDb
181  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
182 
183  forAll(mesh_.cells(), cellI)
184  {
185  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
186  state[cellI] = states[localIdx];
187  }
188  }
189 
190  forAll(stateInfo_["surfaceScalarStates"], idxI)
191  {
192  // lookup state from meshDb
193  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
194 
195  forAll(mesh_.faces(), faceI)
196  {
197  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
198  if (faceI < daIndex_.nLocalInternalFaces)
199  {
200  state[faceI] = states[localIdx];
201  }
202  else
203  {
204  label relIdx = faceI - daIndex_.nLocalInternalFaces;
205  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
206  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
207  state.boundaryFieldRef()[patchIdx][faceIdx] = states[localIdx];
208  }
209  }
210  }
211 }
212 
213 void DAField::ofField2State(scalar* states) const
214 {
215  /*
216  Description:
217  Assign values OpenFOAM field values to the state variable array
218 
219  Input:
220  OpenFoam field variables
221 
222  Output:
223  states: state variable array
224 
225  Example:
226  Image we have two state variables (p,T) and five cells, running on two CPU
227  processors, the proc0 owns two cells and the proc1 owns three cells,
228  then calling this function will assign the p, and T based on the the state
229  vector (state-by-state ordering):
230 
231  stateVec = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
232  0 1 2 3 | 4 5 6 7 8 9 <- global state vec index
233  ---- proc0 -----|--------- proc1 -------
234  */
235 
236  const objectRegistry& db = mesh_.thisDb();
237 
238  forAll(stateInfo_["volVectorStates"], idxI)
239  {
240  // lookup state from meshDb
241  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
242 
243  forAll(mesh_.cells(), cellI)
244  {
245  for (label comp = 0; comp < 3; comp++)
246  {
247  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
248  states[localIdx] = state[cellI][comp];
249  }
250  }
251  }
252 
253  forAll(stateInfo_["volScalarStates"], idxI)
254  {
255  // lookup state from meshDb
256  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
257 
258  forAll(mesh_.cells(), cellI)
259  {
260  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
261  states[localIdx] = state[cellI];
262  }
263  }
264 
265  forAll(stateInfo_["modelStates"], idxI)
266  {
267  // lookup state from meshDb
268  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
269 
270  forAll(mesh_.cells(), cellI)
271  {
272  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
273  states[localIdx] = state[cellI];
274  }
275  }
276 
277  forAll(stateInfo_["surfaceScalarStates"], idxI)
278  {
279  // lookup state from meshDb
280  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
281 
282  forAll(mesh_.faces(), faceI)
283  {
284  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
285  if (faceI < daIndex_.nLocalInternalFaces)
286  {
287  states[localIdx] = state[faceI];
288  }
289  else
290  {
291  label relIdx = faceI - daIndex_.nLocalInternalFaces;
292  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
293  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
294  states[localIdx] = state.boundaryField()[patchIdx][faceIdx];
295  }
296  }
297  }
298 }
299 
300 void DAField::stateVec2OFField(const Vec stateVec) const
301 {
302  /*
303  Description:
304  Assign values OpenFOAM field values based on the state variable vector
305 
306  Input:
307  stateVec: state variable vector
308 
309  Output:
310  OpenFoam field variables
311 
312  Example:
313  Image we have two state variables (p,T) and five cells, running on two CPU
314  processors, the proc0 owns two cells and the proc1 owns three cells,
315  then calling this function will assign the p, and T based on the the state
316  vector (state-by-state ordering):
317 
318  stateVec = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
319  0 1 2 3 | 4 5 6 7 8 9 <- global state vec index
320  ---- proc0 -----|--------- proc1 -------
321  */
322 
323  const objectRegistry& db = mesh_.thisDb();
324  const PetscScalar* stateVecArray;
325  VecGetArrayRead(stateVec, &stateVecArray);
326 
327  forAll(stateInfo_["volVectorStates"], idxI)
328  {
329  // lookup state from meshDb
330  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
331 
332  forAll(mesh_.cells(), cellI)
333  {
334  for (label comp = 0; comp < 3; comp++)
335  {
336  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
337  state[cellI][comp] = stateVecArray[localIdx];
338  }
339  }
340  }
341 
342  forAll(stateInfo_["volScalarStates"], idxI)
343  {
344  // lookup state from meshDb
345  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
346 
347  forAll(mesh_.cells(), cellI)
348  {
349  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
350  state[cellI] = stateVecArray[localIdx];
351  }
352  }
353 
354  forAll(stateInfo_["modelStates"], idxI)
355  {
356  // lookup state from meshDb
357  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
358 
359  forAll(mesh_.cells(), cellI)
360  {
361  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
362  state[cellI] = stateVecArray[localIdx];
363  }
364  }
365 
366  forAll(stateInfo_["surfaceScalarStates"], idxI)
367  {
368  // lookup state from meshDb
369  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
370 
371  forAll(mesh_.faces(), faceI)
372  {
373  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
374  if (faceI < daIndex_.nLocalInternalFaces)
375  {
376  state[faceI] = stateVecArray[localIdx];
377  }
378  else
379  {
380  label relIdx = faceI - daIndex_.nLocalInternalFaces;
381  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
382  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
383  state.boundaryFieldRef()[patchIdx][faceIdx] = stateVecArray[localIdx];
384  }
385  }
386  }
387  VecRestoreArrayRead(stateVec, &stateVecArray);
388 }
389 
390 void DAField::point2OFMesh(const scalar* volCoords) const
391 {
392  /*
393  Description:
394  Assign the points in fvMesh of OpenFOAM based on the volCoords array
395 
396  Input:
397  volCoords: a vector that stores the x, y, and z coordinates for all
398  points in the fvMesh mesh
399 
400  Output:
401  New mesh metrics in fvMesh, effectively by calling mesh.movePoints
402 
403  Example:
404  Image we have three points in fvMesh, running on two CPU
405  processors, the proc0 owns one point and the proc1 owns two points,
406  then calling this function will assign xvVec based on the the points
407  coordinates in fvMesh
408 
409  xvVec = [x0, y0, z0 | x0, y0, z0, x1, y1, z1] <- x0 means x coordinate for the 0th point on local processor
410  0 1 2 | 3 4 5 6 7 8 <- global point vec index
411  --- proc0 --|--------- proc1 -------
412  */
413 
414  pointField meshPoints(mesh_.points());
415 
416  forAll(mesh_.points(), pointI)
417  {
418  for (label comp = 0; comp < 3; comp++)
419  {
420  label localIdx = daIndex_.getLocalXvIndex(pointI, comp);
421  meshPoints[pointI][comp] = volCoords[localIdx];
422  }
423  }
424 
425  // movePoints update the mesh metrics such as volume, surface area and cell centers
426  fvMesh& mesh = const_cast<fvMesh&>(mesh_);
427  mesh.movePoints(meshPoints);
428  if (daOption_.getAllOptions().subDict("dynamicMesh").getLabel("active"))
429  {
430  mesh.moving(true);
431  }
432  else
433  {
434  mesh.moving(false);
435  }
436 }
437 
438 void DAField::ofMesh2PointVec(Vec xvVec) const
439 {
440  /*
441  Description:
442  Assign the point vector based on the points in fvMesh of OpenFOAM
443 
444  Input:
445  Mesh coordinates in fvMesh
446 
447  Output:
448  xvVec: a vector that stores the x, y, and z coordinates for all
449  points in the fvMesh mesh
450 
451  Example:
452  Image we have three points in fvMesh, running on two CPU
453  processors, the proc0 owns one point and the proc1 owns two points,
454  then calling this function will assign xvVec based on the the points
455  coordinates in fvMesh
456 
457  xvVec = [x0, y0, z0 | x0, y0, z0, x1, y1, z1] <- x0 means x coordinate for the 0th point on local processor
458  0 1 2 | 3 4 5 6 7 8 <- global point vec index
459  --- proc0 --|--------- proc1 -------
460  */
461 
462  PetscScalar* xvVecArray;
463  VecGetArray(xvVec, &xvVecArray);
464 
465  forAll(mesh_.points(), pointI)
466  {
467  for (label comp = 0; comp < 3; comp++)
468  {
469  label localIdx = daIndex_.getLocalXvIndex(pointI, comp);
470  assignValueCheckAD(xvVecArray[localIdx], mesh_.points()[pointI][comp]);
471  }
472  }
473 
474  VecRestoreArray(xvVec, &xvVecArray);
475 }
476 
477 void DAField::ofResField2ResVec(Vec resVec) const
478 {
479  /*
480  Description:
481  Assign values for the residual vector based on the
482  latest OpenFOAM residual field values
483 
484  Input:
485  OpenFOAM residual field variables
486 
487  Output:
488  resVec: state residual vector
489 
490  Example:
491  Image we have two state residuals (pRes,TRes) and five cells, running on two CPU
492  processors, the proc0 owns two cells and the proc1 owns three cells,
493  then calling this function gives the residual vector (state-by-state ordering):
494 
495  resVec = [pRes0, pRes1, TRes0, TRes1 | pRes0, pRes1, pRes2, TRes0, TRes1, TRes2]
496  0 1 2 3 | 4 5 6 7 8 9 <- global residual vec index
497  ---------- proc0 ---------|------------- proc1 ----------------------
498  NOTE: pRes0 means p residual for the 0th cell on local processor
499  */
500 
501  const objectRegistry& db = mesh_.thisDb();
502  PetscScalar* stateResVecArray;
503  VecGetArray(resVec, &stateResVecArray);
504 
505  forAll(stateInfo_["volVectorStates"], idxI)
506  {
507  // lookup state from meshDb
508  makeStateRes(stateInfo_["volVectorStates"][idxI], volVectorField, db);
509 
510  forAll(mesh_.cells(), cellI)
511  {
512  for (label comp = 0; comp < 3; comp++)
513  {
514  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
515  assignValueCheckAD(stateResVecArray[localIdx], stateRes[cellI][comp]);
516  }
517  }
518  }
519 
520  forAll(stateInfo_["volScalarStates"], idxI)
521  {
522  // lookup state from meshDb
523  makeStateRes(stateInfo_["volScalarStates"][idxI], volScalarField, db);
524 
525  forAll(mesh_.cells(), cellI)
526  {
527  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
528  assignValueCheckAD(stateResVecArray[localIdx], stateRes[cellI]);
529  }
530  }
531 
532  forAll(stateInfo_["modelStates"], idxI)
533  {
534  // lookup state from meshDb
535  makeStateRes(stateInfo_["modelStates"][idxI], volScalarField, db);
536 
537  forAll(mesh_.cells(), cellI)
538  {
539  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
540  assignValueCheckAD(stateResVecArray[localIdx], stateRes[cellI]);
541  }
542  }
543 
544  forAll(stateInfo_["surfaceScalarStates"], idxI)
545  {
546  // lookup state from meshDb
547  makeStateRes(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
548 
549  forAll(mesh_.faces(), faceI)
550  {
551  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
552  if (faceI < daIndex_.nLocalInternalFaces)
553  {
554  assignValueCheckAD(stateResVecArray[localIdx], stateRes[faceI]);
555  }
556  else
557  {
558  label relIdx = faceI - daIndex_.nLocalInternalFaces;
559  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
560  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
561  assignValueCheckAD(stateResVecArray[localIdx], stateRes.boundaryField()[patchIdx][faceIdx]);
562  }
563  }
564  }
565  VecRestoreArray(resVec, &stateResVecArray);
566 }
567 
568 void DAField::checkSpecialBCs()
569 {
570  /*
571  Description:
572  Check if we need to do special treatment for boundary conditions
573  If any special BC is detected, append their names to specialBCs list
574  */
575 
576  // *******************************************************************
577  // pressureInletVelocity
578  // *******************************************************************
579  // Note we need to read the U field, instead of getting it from db
580  // this is because coloringSolver does not read U
581  // Also we need to set read_if_present for solid solvers in which
582  // there is no U field
583  volVectorField U(
584  IOobject(
585  "U",
586  mesh_.time().timeName(),
587  mesh_,
588  IOobject::READ_IF_PRESENT,
589  IOobject::NO_WRITE,
590  false),
591  mesh_,
592  dimensionedVector("U", dimensionSet(0, 1, -1, 0, 0, 0, 0), vector::zero),
593  zeroGradientFvPatchField<vector>::typeName);
594  forAll(U.boundaryField(), patchI)
595  {
596  if (U.boundaryFieldRef()[patchI].type() == "pressureInletVelocity")
597  {
598  specialBCs.append("pressureInletVelocity");
599  break;
600  }
601  }
602 
603  // *******************************************************************
604  // append more special BCs
605  // *******************************************************************
606 }
607 
608 void DAField::specialBCTreatment()
609 {
610  /*
611  Description:
612  Apply special treatment for boundary conditions
613  */
614 
615  // *******************************************************************
616  // pressureInletVelocity
617  // *******************************************************************
618  // for pressureInletVelocity, the inlet U depends on
619  // rho and phi, so we need to call U.correctBoundaryConditions again
620  if (specialBCs.found("pressureInletVelocity"))
621  {
622 
623  volVectorField& U(const_cast<volVectorField&>(
624  mesh_.thisDb().lookupObject<volVectorField>("U")));
625  U.correctBoundaryConditions();
626  }
627 
628  // *******************************************************************
629  // append more special BCs
630  // *******************************************************************
631 }
632 
633 void DAField::setPrimalBoundaryConditions(const label printInfo)
634 {
635  /*
636  Description:
637  A general function to read the inlet/outlet values from DAOption, and set
638  the corresponding values to the boundary field. It also setup turbulence
639  wall boundary condition
640  Note: this function should be called before running the primal solver
641  If nothing is set, the BC will remain unchanged
642  Example
643  "primalBC":
644  {
645  "bc0":
646  {
647  "patches": ["inlet", "side"],
648  "variable": "U",
649  "value": [10.0 0.0, 0.0],
650  },
651  "bc1":
652  {
653  "patches": ["outlet"],
654  "variable": "p",
655  "value": [101325.0],
656  },
657  "useWallFunction": True,
658  "MRF": 1000.0,
659  "fvSource": {"value": 0.01, "comp": 0}
660  }
661  */
662 
663  const objectRegistry& db = mesh_.thisDb();
664 
665  label setTurbWallBCs = 0;
666  label useWallFunction = 0;
667 
668  const dictionary& allOptions = daOption_.getAllOptions();
669  dictionary bcDict = allOptions.subDict("primalBC");
670 
671  forAll(bcDict.toc(), idxI)
672  {
673  word bcKey = bcDict.toc()[idxI];
674 
675  if (bcKey == "useWallFunction")
676  {
677  setTurbWallBCs = 1;
678  useWallFunction = bcDict.getLabel("useWallFunction");
679  continue;
680  }
681  else if (bcKey == "MRF")
682  {
683  // change the rotation speed in MRF
684  scalar omegaNew = bcDict.getScalar("MRF");
685  const IOMRFZoneListDF& MRF = db.lookupObject<IOMRFZoneListDF>("MRFProperties");
686  scalar& omega = const_cast<scalar&>(MRF.getOmegaRef());
687  omega = omegaNew;
688 
689  if (printInfo)
690  {
691  Info << "Setting MRF omega to " << omegaNew << endl;
692  }
693 
694  continue;
695  }
696  else if (bcKey == "transport:nu")
697  {
698  // change the nu field
699  scalar nu = bcDict.getScalar("transport:nu");
700  volScalarField& nuField = const_cast<volScalarField&>(
701  db.lookupObject<volScalarField>("nu"));
702  forAll(nuField, cellI)
703  {
704  nuField[cellI] = nu;
705  }
706  forAll(nuField.boundaryField(), patchI)
707  {
708  forAll(nuField.boundaryField()[patchI], faceI)
709  {
710  nuField.boundaryFieldRef()[patchI][faceI] = nu;
711  }
712  }
713  nuField.correctBoundaryConditions();
714 
715  if (printInfo)
716  {
717  Info << "Setting transportProperties nu to " << nu << endl;
718  }
719 
720  continue;
721  }
722  else if (bcKey == "thermo:mu")
723  {
724  // change the nu field
725  scalar mu = bcDict.getScalar("thermo:mu");
726  volScalarField& muField = const_cast<volScalarField&>(
727  db.lookupObject<volScalarField>("thermo:mu"));
728  forAll(muField, cellI)
729  {
730  muField[cellI] = mu;
731  }
732  forAll(muField.boundaryField(), patchI)
733  {
734  forAll(muField.boundaryField()[patchI], faceI)
735  {
736  muField.boundaryFieldRef()[patchI][faceI] = mu;
737  }
738  }
739  muField.correctBoundaryConditions();
740 
741  if (printInfo)
742  {
743  Info << "Setting thermalphysicalProperties mu to " << mu << endl;
744  }
745 
746  continue;
747  }
748 
749  dictionary bcSubDict = bcDict.subDict(bcKey);
750 
751  wordList patches;
752  bcSubDict.readEntry<wordList>("patches", patches);
753  word variable = bcSubDict.getWord("variable");
754  scalarList value;
755  bcSubDict.readEntry<scalarList>("value", value);
756 
757  // loop over all patches and set values
758  forAll(patches, idxI)
759  {
760  word patch = patches[idxI];
761 
762  // it should be a scalar
763  if (value.size() == 1)
764  {
765  if (!db.foundObject<volScalarField>(variable))
766  {
767  if (printInfo)
768  {
769  Info << variable << " not found, skip it." << endl;
770  }
771  continue;
772  }
773  // it is a scalar
774  volScalarField& state(const_cast<volScalarField&>(
775  db.lookupObject<volScalarField>(variable)));
776 
777  if (printInfo)
778  {
779  Info << "Setting primal boundary conditions..." << endl;
780  Info << "Setting " << variable << " = " << value[0] << " at " << patch << endl;
781  }
782 
783  label patchI = mesh_.boundaryMesh().findPatchID(patch);
784 
785  // for decomposed domain, don't set BC if the patch is empty
786  if (mesh_.boundaryMesh()[patchI].size() > 0)
787  {
788  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
789  {
790  forAll(state.boundaryFieldRef()[patchI], faceI)
791  {
792  state.boundaryFieldRef()[patchI][faceI] = value[0];
793  }
794  }
795  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
796  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
797  {
798  // set value
799  forAll(state.boundaryFieldRef()[patchI], faceI)
800  {
801  state.boundaryFieldRef()[patchI][faceI] = value[0];
802  }
803  // set inletValue
804  mixedFvPatchField<scalar>& inletOutletPatch =
805  refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
806 
807  inletOutletPatch.refValue() = value[0];
808  }
809  else if (state.boundaryFieldRef()[patchI].type() == "fixedGradient")
810  {
811  fixedGradientFvPatchField<scalar>& patchBC =
812  refCast<fixedGradientFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
813  scalarField& grad = const_cast<scalarField&>(patchBC.gradient());
814  forAll(grad, idxI)
815  {
816  grad[idxI] = value[0];
817  }
818  }
819  else
820  {
821  FatalErrorIn("") << "only support fixedValues, inletOutlet, "
822  << "outletInlet, fixedGradient!" << abort(FatalError);
823  }
824  }
825  }
826  else if (value.size() == 3)
827  {
828  if (!db.foundObject<volVectorField>(variable))
829  {
830  if (printInfo)
831  {
832  Info << variable << " not found, skip it." << endl;
833  }
834  continue;
835  }
836  // it is a vector
837  volVectorField& state(const_cast<volVectorField&>(
838  db.lookupObject<volVectorField>(variable)));
839 
840  vector valVec = {value[0], value[1], value[2]};
841  if (printInfo)
842  {
843  Info << "Setting primal boundary conditions..." << endl;
844  Info << "Setting " << variable << " = (" << value[0] << " "
845  << value[1] << " " << value[2] << ") at " << patch << endl;
846  }
847 
848  label patchI = mesh_.boundaryMesh().findPatchID(patch);
849 
850  // for decomposed domain, don't set BC if the patch is empty
851  if (mesh_.boundaryMesh()[patchI].size() > 0)
852  {
853  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
854  {
855  forAll(state.boundaryFieldRef()[patchI], faceI)
856  {
857  state.boundaryFieldRef()[patchI][faceI] = valVec;
858  }
859  }
860  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
861  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
862  {
863  // set value
864  forAll(state.boundaryFieldRef()[patchI], faceI)
865  {
866  state.boundaryFieldRef()[patchI][faceI] = valVec;
867  }
868  // set inletValue
869  mixedFvPatchField<vector>& inletOutletPatch =
870  refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
871  inletOutletPatch.refValue() = valVec;
872  }
873  else if (state.boundaryFieldRef()[patchI].type() == "fixedGradient")
874  {
875  fixedGradientFvPatchField<vector>& patchBC =
876  refCast<fixedGradientFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
877  vectorField& grad = const_cast<vectorField&>(patchBC.gradient());
878  forAll(grad, idxI)
879  {
880  grad[idxI][0] = value[0];
881  grad[idxI][1] = value[1];
882  grad[idxI][2] = value[2];
883  }
884  }
885  else
886  {
887  FatalErrorIn("") << "only support fixedValues, inletOutlet, "
888  << "fixedGradient, and tractionDisplacement!"
889  << abort(FatalError);
890  }
891  }
892  }
893  else
894  {
895  FatalErrorIn("") << "value should be a list of either 1 (scalar) "
896  << "or 3 (vector) elements" << abort(FatalError);
897  }
898  }
899  }
900 
901  // we also set wall boundary conditions for turbulence variables
902  if (setTurbWallBCs)
903  {
904 
905  IOdictionary turbDict(
906  IOobject(
907  "turbulenceProperties",
908  mesh_.time().constant(),
909  mesh_,
910  IOobject::MUST_READ,
911  IOobject::NO_WRITE,
912  false));
913 
914  dictionary coeffDict(turbDict.subDict("RAS"));
915  word turbModelType = word(coeffDict["RASModel"]);
916  //Info<<"turbModelType: "<<turbModelType<<endl;
917 
918  // create word regular expression for SA model
919  wordReList SAModelWordReList {
920  {"SpalartAllmaras.*", wordRe::REGEX}};
921  wordRes SAModelWordRes(SAModelWordReList);
922 
923  // ------ nut ----------
924  if (db.foundObject<volScalarField>("nut"))
925  {
926 
927  volScalarField& nut(const_cast<volScalarField&>(
928  db.lookupObject<volScalarField>("nut")));
929 
930  forAll(nut.boundaryField(), patchI)
931  {
932  if (mesh_.boundaryMesh()[patchI].type() == "wall")
933  {
934  if (printInfo)
935  {
936  Info << "Setting nut wall BC for "
937  << mesh_.boundaryMesh()[patchI].name() << ". ";
938  }
939 
940  if (useWallFunction)
941  {
942  // wall function for SA
943  if (SAModelWordRes(turbModelType))
944  {
945  nut.boundaryFieldRef().set(
946  patchI,
947  fvPatchField<scalar>::New(
948  "nutUSpaldingWallFunction",
949  mesh_.boundary()[patchI],
950  nut));
951 
952  if (printInfo)
953  {
954  Info << "BCType=nutUSpaldingWallFunction" << endl;
955  }
956  }
957  else // wall function for kOmega and kEpsilon
958  {
959  nut.boundaryFieldRef().set(
960  patchI,
961  fvPatchField<scalar>::New(
962  "nutkWallFunction",
963  mesh_.boundary()[patchI],
964  nut));
965 
966  if (printInfo)
967  {
968  Info << "BCType=nutkWallFunction" << endl;
969  }
970  }
971 
972  // set boundary values
973  // for decomposed domain, don't set BC if the patch is empty
974  if (mesh_.boundaryMesh()[patchI].size() > 0)
975  {
976  scalar wallVal = nut[0];
977  forAll(nut.boundaryFieldRef()[patchI], faceI)
978  {
979  // assign uniform field
980  nut.boundaryFieldRef()[patchI][faceI] = wallVal;
981  }
982  }
983  }
984  else
985  {
986  nut.boundaryFieldRef().set(
987  patchI,
988  fvPatchField<scalar>::New(
989  "nutLowReWallFunction",
990  mesh_.boundary()[patchI],
991  nut));
992 
993  if (printInfo)
994  {
995  Info << "BCType=nutLowReWallFunction" << endl;
996  }
997 
998  // set boundary values
999  // for decomposed domain, don't set BC if the patch is empty
1000  if (mesh_.boundaryMesh()[patchI].size() > 0)
1001  {
1002  forAll(nut.boundaryFieldRef()[patchI], faceI)
1003  {
1004  // assign uniform field
1005  nut.boundaryFieldRef()[patchI][faceI] = 1e-14;
1006  }
1007  }
1008  }
1009  }
1010  }
1011  }
1012 
1013  // ------ k ----------
1014  if (db.foundObject<volScalarField>("k"))
1015  {
1016 
1017  volScalarField& k(const_cast<volScalarField&>(
1018  db.lookupObject<volScalarField>("k")));
1019 
1020  forAll(k.boundaryField(), patchI)
1021  {
1022  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1023  {
1024  if (printInfo)
1025  {
1026  Info << "Setting k wall BC for "
1027  << mesh_.boundaryMesh()[patchI].name() << ". ";
1028  }
1029 
1030  if (useWallFunction)
1031  {
1032  // wall function for SA
1033  k.boundaryFieldRef().set(
1034  patchI,
1035  fvPatchField<scalar>::New("kqRWallFunction", mesh_.boundary()[patchI], k));
1036 
1037  if (printInfo)
1038  {
1039  Info << "BCType=kqRWallFunction" << endl;
1040  }
1041 
1042  // set boundary values
1043  // for decomposed domain, don't set BC if the patch is empty
1044  if (mesh_.boundaryMesh()[patchI].size() > 0)
1045  {
1046  scalar wallVal = k[0];
1047  forAll(k.boundaryFieldRef()[patchI], faceI)
1048  {
1049  k.boundaryFieldRef()[patchI][faceI] = wallVal; // assign uniform field
1050  }
1051  }
1052  }
1053  else
1054  {
1055  k.boundaryFieldRef().set(
1056  patchI,
1057  fvPatchField<scalar>::New("fixedValue", mesh_.boundary()[patchI], k));
1058 
1059  if (printInfo)
1060  {
1061  Info << "BCType=fixedValue" << endl;
1062  }
1063 
1064  // set boundary values
1065  // for decomposed domain, don't set BC if the patch is empty
1066  if (mesh_.boundaryMesh()[patchI].size() > 0)
1067  {
1068  forAll(k.boundaryFieldRef()[patchI], faceI)
1069  {
1070  k.boundaryFieldRef()[patchI][faceI] = 1e-14; // assign uniform field
1071  }
1072  }
1073  }
1074  }
1075  }
1076  }
1077 
1078  // ------ omega ----------
1079  if (db.foundObject<volScalarField>("omega"))
1080  {
1081 
1082  volScalarField& omega(const_cast<volScalarField&>(
1083  db.lookupObject<volScalarField>("omega")));
1084 
1085  forAll(omega.boundaryField(), patchI)
1086  {
1087  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1088  {
1089  if (printInfo)
1090  {
1091  Info << "Setting omega wall BC for "
1092  << mesh_.boundaryMesh()[patchI].name() << ". ";
1093  }
1094 
1095  // always use omegaWallFunction
1096  omega.boundaryFieldRef().set(
1097  patchI,
1098  fvPatchField<scalar>::New("omegaWallFunction", mesh_.boundary()[patchI], omega));
1099 
1100  if (printInfo)
1101  {
1102  Info << "BCType=omegaWallFunction" << endl;
1103  }
1104 
1105  // set boundary values
1106  // for decomposed domain, don't set BC if the patch is empty
1107  if (mesh_.boundaryMesh()[patchI].size() > 0)
1108  {
1109  scalar wallVal = omega[0];
1110  forAll(omega.boundaryFieldRef()[patchI], faceI)
1111  {
1112  omega.boundaryFieldRef()[patchI][faceI] = wallVal; // assign uniform field
1113  }
1114  }
1115  }
1116  }
1117  }
1118 
1119  // ------ epsilon ----------
1120  if (db.foundObject<volScalarField>("epsilon"))
1121  {
1122 
1123  volScalarField& epsilon(const_cast<volScalarField&>(
1124  db.lookupObject<volScalarField>("epsilon")));
1125 
1126  forAll(epsilon.boundaryField(), patchI)
1127  {
1128  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1129  {
1130 
1131  if (printInfo)
1132  {
1133  Info << "Setting epsilon wall BC for "
1134  << mesh_.boundaryMesh()[patchI].name() << ". ";
1135  }
1136 
1137  if (useWallFunction)
1138  {
1139  epsilon.boundaryFieldRef().set(
1140  patchI,
1141  fvPatchField<scalar>::New("epsilonWallFunction", mesh_.boundary()[patchI], epsilon));
1142 
1143  if (printInfo)
1144  {
1145  Info << "BCType=epsilonWallFunction" << endl;
1146  }
1147 
1148  // set boundary values
1149  // for decomposed domain, don't set BC if the patch is empty
1150  if (mesh_.boundaryMesh()[patchI].size() > 0)
1151  {
1152  scalar wallVal = epsilon[0];
1153  forAll(epsilon.boundaryFieldRef()[patchI], faceI)
1154  {
1155  epsilon.boundaryFieldRef()[patchI][faceI] = wallVal; // assign uniform field
1156  }
1157  }
1158  }
1159  else
1160  {
1161  epsilon.boundaryFieldRef().set(
1162  patchI,
1163  fvPatchField<scalar>::New("fixedValue", mesh_.boundary()[patchI], epsilon));
1164 
1165  if (printInfo)
1166  {
1167  Info << "BCType=fixedValue" << endl;
1168  }
1169 
1170  // set boundary values
1171  // for decomposed domain, don't set BC if the patch is empty
1172  if (mesh_.boundaryMesh()[patchI].size() > 0)
1173  {
1174  forAll(epsilon.boundaryFieldRef()[patchI], faceI)
1175  {
1176  epsilon.boundaryFieldRef()[patchI][faceI] = 1e-14; // assign uniform field
1177  }
1178  }
1179  }
1180  }
1181  }
1182  }
1183 
1184  // ------ alphat ----------
1185  // TODO: need to figure out a way to pass the Prt parameter
1186  // to the wall function. The default is 0.85 and is hard coded
1187  if (db.foundObject<volScalarField>("alphat"))
1188  {
1189 
1190  volScalarField& alphat(const_cast<volScalarField&>(
1191  db.lookupObject<volScalarField>("alphat")));
1192 
1193  forAll(alphat.boundaryField(), patchI)
1194  {
1195  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1196  {
1197 
1198  if (printInfo)
1199  {
1200  Info << "Setting alphat wall BC for "
1201  << mesh_.boundaryMesh()[patchI].name();
1202  }
1203 
1204  if (useWallFunction)
1205  {
1206 
1207  if (mesh_.thisDb().foundObject<incompressible::turbulenceModel>(incompressible::turbulenceModel::propertiesName))
1208  {
1209  // incompressible case
1210  alphat.boundaryFieldRef().set(
1211  patchI,
1212  fvPatchField<scalar>::New("incompressible::alphatWallFunction", mesh_.boundary()[patchI], alphat));
1213 
1214  if (printInfo)
1215  {
1216  Info << "BCType=incompressible::alphatWallFunction. Default Prt=0.85" << endl;
1217  }
1218  }
1219  else if (mesh_.thisDb().foundObject<compressible::turbulenceModel>(compressible::turbulenceModel::propertiesName))
1220  {
1221  // compressible case
1222  // incompressible case
1223  alphat.boundaryFieldRef().set(
1224  patchI,
1225  fvPatchField<scalar>::New("compressible::alphatWallFunction", mesh_.boundary()[patchI], alphat));
1226 
1227  if (printInfo)
1228  {
1229  Info << "BCType=compressible::alphatWallFunction. Default Prt=0.85" << endl;
1230  }
1231  }
1232  else
1233  {
1234  FatalErrorIn("") << "turbModelType_ not valid!" << abort(FatalError);
1235  }
1236 
1237  // set boundary values
1238  // for decomposed domain, don't set BC if the patch is empty
1239  if (mesh_.boundaryMesh()[patchI].size() > 0)
1240  {
1241  scalar wallVal = alphat[0];
1242  forAll(alphat.boundaryFieldRef()[patchI], faceI)
1243  {
1244  alphat.boundaryFieldRef()[patchI][faceI] = wallVal; // assign uniform field
1245  }
1246  }
1247  }
1248  else
1249  {
1250  alphat.boundaryFieldRef().set(
1251  patchI,
1252  fvPatchField<scalar>::New("fixedValue", mesh_.boundary()[patchI], alphat));
1253 
1254  if (printInfo)
1255  {
1256  Info << "BCType=fixedValue" << endl;
1257  }
1258 
1259  // set boundary values
1260  // for decomposed domain, don't set BC if the patch is empty
1261  if (mesh_.boundaryMesh()[patchI].size() > 0)
1262  {
1263  forAll(alphat.boundaryFieldRef()[patchI], faceI)
1264  {
1265  alphat.boundaryFieldRef()[patchI][faceI] = 1e-14; // assign uniform field
1266  }
1267  }
1268  }
1269  }
1270  }
1271  }
1272  }
1273 }
1274 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1275 
1276 } // End namespace Foam
1277 
1278 // ************************************************************************* //
mu
volScalarField & mu
Definition: createRefsSolidDisplacement.H:6
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Foam::DAField::daIndex_
const DAIndex & daIndex_
DAIndex object.
Definition: DAField.H:57
Foam::DAStateInfo::New
static autoPtr< DAStateInfo > New(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel)
Definition: DAStateInfo.C:43
Foam::DAField::ofField2State
void ofField2State(scalar *states) const
assign the openfoam fields to the states array
Definition: DAField.C:213
Foam::DAField::setPrimalBoundaryConditions
void setPrimalBoundaryConditions(const label printInfo=1)
set the boundary conditions based on parameters defined in DAOption
Definition: DAField.C:633
Foam::DAIndex::bFacePatchI
labelList bFacePatchI
Definition: DAIndex.H:110
Foam::DAOption
Definition: DAOption.H:29
Foam::DAField::state2OFField
void state2OFField(const scalar *states) const
assign the fields in OpenFOAM based on the state array
Definition: DAField.C:126
Foam::DAField::specialBCs
wordList specialBCs
a list that contains the names of detected special boundary conditions
Definition: DAField.H:108
Foam::DAField::checkSpecialBCs
void checkSpecialBCs()
check if we need to do special treatment for boundary conditions
Definition: DAField.C:568
MRF
IOMRFZoneListDF & MRF
Definition: createRefsRhoSimple.H:18
Foam::DAField::mesh_
const fvMesh & mesh_
Foam::fvMesh object.
Definition: DAField.H:48
Foam::DAOption::getOption
classType getOption(const word key) const
get an option from subDict and key
Definition: DAOption.H:92
solverName
word solverName
Definition: createAdjoint.H:14
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::DAField::ofField2StateVec
void ofField2StateVec(Vec stateVec) const
set the state vector based on the latest fields in OpenFOAM
Definition: DAField.C:36
Foam::DAField::ofMesh2PointVec
void ofMesh2PointVec(Vec xvVec) const
assign the point vector based on the points in fvMesh of OpenFOAM
Definition: DAField.C:438
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAField::point2OFMesh
void point2OFMesh(const scalar *volCoords) const
assign the points in fvMesh of OpenFOAM based on the point array
Definition: DAField.C:390
assignValueCheckAD
#define assignValueCheckAD(valA, valB)
Definition: DAMacroFunctions.H:103
Foam::DAField::daOption_
const DAOption & daOption_
Foam::DAOption object.
Definition: DAField.H:51
Foam::DAModel
Definition: DAModel.H:57
k
dimensionedScalar & k
Definition: createRefsHeatTransfer.H:6
Foam
Definition: checkGeometry.C:32
Foam::DAField::stateVec2OFField
void stateVec2OFField(const Vec stateVec) const
assign the fields in OpenFOAM based on the state vector
Definition: DAField.C:300
allOptions
const dictionary & allOptions
Definition: createRefsPimpleDyM.H:14
Foam::IOMRFZoneListDF
Definition: IOMRFZoneListDF.H:49
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAIndex::nLocalInternalFaces
label nLocalInternalFaces
local internal face size
Definition: DAIndex.H:98
Foam::DAIndex::getLocalXvIndex
label getLocalXvIndex(const label idxPoint, const label idxCoord) const
get local Xv index for a given point index and coordinate component (x, y, or z)
Definition: DAIndex.C:733
makeStateRes
#define makeStateRes(regState, fieldType, db)
Definition: DAMacroFunctions.H:21
alphat
volScalarField & alphat
Definition: TEqnPimpleDyM.H:5
Foam::DAField::ofResField2ResVec
void ofResField2ResVec(Vec resVec) const
assign the residual vector based on the residual field in OpenFOAM
Definition: DAField.C:477
Foam::DAIndex::getLocalAdjointStateIndex
label getLocalAdjointStateIndex(const word stateName, const label idxI, const label comp=-1) const
get local adjoint index for a given state name, cell/face indxI and its component (optional,...
Definition: DAIndex.C:516
makeState
#define makeState(regState, fieldType, db)
Definition: DAMacroFunctions.H:15
DAField.H
Foam::DAField::specialBCTreatment
void specialBCTreatment()
apply special treatment for boundary conditions
Definition: DAField.C:608
Foam::DAIndex::bFaceFaceI
labelList bFaceFaceI
Definition: DAIndex.H:116
Foam::DAField::stateInfo_
HashTable< wordList > stateInfo_
the StateInfo_ list from DAStateInfo object
Definition: DAField.H:60

Generated by   doxygen 1.8.17