DAField.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v3
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::resetOFSeeds()
37 {
38  /*
39  Description:
40  Reset the seeds for both state and mesh variables in OpenFOAM by setting their gradient() to zeros
41  Then, we will propagate the zero seeds to other variables by calling correctBC or movePoints
42  NOTE. we need to zero the seeds for all boundary values as well. For example, in aerothermal coupling
43  we need to set seeds for boundary values and compute the matrix-vector product. We need to reset
44  the seeds for the boundary values after the product is done.
45  TODO. We only reset seeds for boundary types fixedValue and fixedGradient for now. If we need to couple
46  more boundary types in the future, make sure we add codes to reset their seeds here.
47 
48  Output:
49  OpenFoam field variables's gradients will be reset to zeros
50  */
51 
52 #ifdef CODI_AD_REVERSE
53 
54  const objectRegistry& db = mesh_.thisDb();
55 
56  forAll(stateInfo_["volVectorStates"], idxI)
57  {
58  // lookup state from meshDb
59  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
60 
61  forAll(state, cellI)
62  {
63  for (label comp = 0; comp < 3; comp++)
64  {
65  state[cellI][comp].setGradient(0.0);
66  }
67  }
68 
69  forAll(state.boundaryField(), patchI)
70  {
71  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
72  {
73  forAll(state.boundaryField()[patchI], faceI)
74  {
75  for (label comp = 0; comp < 3; comp++)
76  {
77  state.boundaryFieldRef()[patchI][faceI][comp].setGradient(0.0);
78  }
79  }
80  }
81  else if (state.boundaryFieldRef()[patchI].type() == "fixedGradient")
82  {
83  fixedGradientFvPatchField<vector>& patchBC =
84  refCast<fixedGradientFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
85  vectorField& grad = const_cast<vectorField&>(patchBC.gradient());
86  forAll(grad, idxI)
87  {
88  for (label comp = 0; comp < 3; comp++)
89  {
90  grad[idxI][comp].setGradient(0.0);
91  }
92  }
93  }
94  }
95  }
96 
97  forAll(stateInfo_["volScalarStates"], idxI)
98  {
99  // lookup state from meshDb
100  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
101 
102  forAll(state, cellI)
103  {
104  state[cellI].setGradient(0.0);
105  }
106 
107  forAll(state.boundaryField(), patchI)
108  {
109  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
110  {
111  forAll(state.boundaryField()[patchI], faceI)
112  {
113  state.boundaryFieldRef()[patchI][faceI].setGradient(0.0);
114  }
115  }
116  else if (state.boundaryFieldRef()[patchI].type() == "fixedGradient")
117  {
118  fixedGradientFvPatchField<scalar>& patchBC =
119  refCast<fixedGradientFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
120  scalarField& grad = const_cast<scalarField&>(patchBC.gradient());
121  forAll(grad, idxI)
122  {
123  grad[idxI].setGradient(0.0);
124  }
125  }
126  }
127  }
128 
129  forAll(stateInfo_["modelStates"], idxI)
130  {
131  // lookup state from meshDb
132  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
133 
134  forAll(state, cellI)
135  {
136  state[cellI].setGradient(0.0);
137  }
138 
139  forAll(state.boundaryField(), patchI)
140  {
141  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
142  {
143  forAll(state.boundaryField()[patchI], faceI)
144  {
145  state.boundaryFieldRef()[patchI][faceI].setGradient(0.0);
146  }
147  }
148  else if (state.boundaryFieldRef()[patchI].type() == "fixedGradient")
149  {
150  fixedGradientFvPatchField<scalar>& patchBC =
151  refCast<fixedGradientFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
152  scalarField& grad = const_cast<scalarField&>(patchBC.gradient());
153  forAll(grad, idxI)
154  {
155  grad[idxI].setGradient(0.0);
156  }
157  }
158  }
159  }
160 
161  forAll(stateInfo_["surfaceScalarStates"], idxI)
162  {
163  // lookup state from meshDb
164  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
165 
166  forAll(state, faceI)
167  {
168  state[faceI].setGradient(0.0);
169  }
170 
171  forAll(state.boundaryField(), patchI)
172  {
173  forAll(state.boundaryField()[patchI], faceI)
174  {
175  state.boundaryFieldRef()[patchI][faceI].setGradient(0.0);
176  }
177  }
178  }
179 
180  pointField meshPoints(mesh_.points());
181 
182  forAll(mesh_.points(), pointI)
183  {
184  for (label comp = 0; comp < 3; comp++)
185  {
186  meshPoints[pointI][comp].setGradient(0.0);
187  }
188  }
189 
190  // movePoints update the mesh metrics such as volume, surface area and cell centers
191  fvMesh& mesh = const_cast<fvMesh&>(mesh_);
192  mesh.movePoints(meshPoints);
193  mesh.moving(false);
194 
195 #endif
196 }
197 
198 void DAField::ofField2StateVec(Vec stateVec) const
199 {
200  /*
201  Description:
202  Assign values for the state variable vector based on the
203  latest OpenFOAM field values
204 
205  Input:
206  OpenFOAM field variables
207 
208  Output:
209  stateVec: state variable vector
210 
211  Example:
212  Image we have two state variables (p,T) and five cells, running on two CPU
213  processors, the proc0 owns two cells and the proc1 owns three cells,
214  then calling this function gives the state vector (state-by-state ordering):
215 
216  stateVec = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
217  0 1 2 3 | 4 5 6 7 8 9 <- global state vec index
218  ---- proc0 -----|--------- proc1 -------
219  */
220 
221  const objectRegistry& db = mesh_.thisDb();
222  PetscScalar* stateVecArray;
223  VecGetArray(stateVec, &stateVecArray);
224 
225  forAll(stateInfo_["volVectorStates"], idxI)
226  {
227  // lookup state from meshDb
228  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
229 
230  forAll(mesh_.cells(), cellI)
231  {
232  for (label comp = 0; comp < 3; comp++)
233  {
234  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
235  assignValueCheckAD(stateVecArray[localIdx], state[cellI][comp]);
236  }
237  }
238  }
239 
240  forAll(stateInfo_["volScalarStates"], idxI)
241  {
242  // lookup state from meshDb
243  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
244 
245  forAll(mesh_.cells(), cellI)
246  {
247  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
248  assignValueCheckAD(stateVecArray[localIdx], state[cellI]);
249  }
250  }
251 
252  forAll(stateInfo_["modelStates"], idxI)
253  {
254  // lookup state from meshDb
255  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
256 
257  forAll(mesh_.cells(), cellI)
258  {
259  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
260  assignValueCheckAD(stateVecArray[localIdx], state[cellI]);
261  }
262  }
263 
264  forAll(stateInfo_["surfaceScalarStates"], idxI)
265  {
266  // lookup state from meshDb
267  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
268 
269  forAll(mesh_.faces(), faceI)
270  {
271  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
272  if (faceI < daIndex_.nLocalInternalFaces)
273  {
274  assignValueCheckAD(stateVecArray[localIdx], state[faceI]);
275  }
276  else
277  {
278  label relIdx = faceI - daIndex_.nLocalInternalFaces;
279  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
280  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
281  assignValueCheckAD(stateVecArray[localIdx], state.boundaryField()[patchIdx][faceIdx]);
282  }
283  }
284  }
285  VecRestoreArray(stateVec, &stateVecArray);
286 }
287 
288 void DAField::state2OFField(const scalar* states) const
289 {
290  /*
291  Description:
292  Assign values OpenFOAM field values based on the state variable array
293 
294  Input:
295  states: state variable array
296 
297  Output:
298  OpenFoam field variables
299 
300  Example:
301  Image we have two state variables (p,T) and five cells, running on two CPU
302  processors, the proc0 owns two cells and the proc1 owns three cells,
303  then calling this function will assign the p, and T based on the the state
304  vector (state-by-state ordering):
305 
306  stateVec = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
307  0 1 2 3 | 4 5 6 7 8 9 <- global state vec index
308  ---- proc0 -----|--------- proc1 -------
309  */
310 
311  const objectRegistry& db = mesh_.thisDb();
312 
313  forAll(stateInfo_["volVectorStates"], idxI)
314  {
315  // lookup state from meshDb
316  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
317 
318  forAll(mesh_.cells(), cellI)
319  {
320  for (label comp = 0; comp < 3; comp++)
321  {
322  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
323  state[cellI][comp] = states[localIdx];
324  }
325  }
326  }
327 
328  forAll(stateInfo_["volScalarStates"], idxI)
329  {
330  // lookup state from meshDb
331  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
332 
333  forAll(mesh_.cells(), cellI)
334  {
335  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
336  state[cellI] = states[localIdx];
337  }
338  }
339 
340  forAll(stateInfo_["modelStates"], idxI)
341  {
342  // lookup state from meshDb
343  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
344 
345  forAll(mesh_.cells(), cellI)
346  {
347  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
348  state[cellI] = states[localIdx];
349  }
350  }
351 
352  forAll(stateInfo_["surfaceScalarStates"], idxI)
353  {
354  // lookup state from meshDb
355  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
356 
357  forAll(mesh_.faces(), faceI)
358  {
359  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
360  if (faceI < daIndex_.nLocalInternalFaces)
361  {
362  state[faceI] = states[localIdx];
363  }
364  else
365  {
366  label relIdx = faceI - daIndex_.nLocalInternalFaces;
367  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
368  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
369  state.boundaryFieldRef()[patchIdx][faceIdx] = states[localIdx];
370  }
371  }
372  }
373 }
374 
375 void DAField::stateVec2OFField(const Vec stateVec) const
376 {
377  /*
378  Description:
379  Assign values OpenFOAM field values based on the state variable vector
380 
381  Input:
382  stateVec: state variable vector
383 
384  Output:
385  OpenFoam field variables
386 
387  Example:
388  Image we have two state variables (p,T) and five cells, running on two CPU
389  processors, the proc0 owns two cells and the proc1 owns three cells,
390  then calling this function will assign the p, and T based on the the state
391  vector (state-by-state ordering):
392 
393  stateVec = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
394  0 1 2 3 | 4 5 6 7 8 9 <- global state vec index
395  ---- proc0 -----|--------- proc1 -------
396  */
397 
398  const objectRegistry& db = mesh_.thisDb();
399  const PetscScalar* stateVecArray;
400  VecGetArrayRead(stateVec, &stateVecArray);
401 
402  forAll(stateInfo_["volVectorStates"], idxI)
403  {
404  // lookup state from meshDb
405  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
406 
407  forAll(mesh_.cells(), cellI)
408  {
409  for (label comp = 0; comp < 3; comp++)
410  {
411  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
412  state[cellI][comp] = stateVecArray[localIdx];
413  }
414  }
415  }
416 
417  forAll(stateInfo_["volScalarStates"], idxI)
418  {
419  // lookup state from meshDb
420  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
421 
422  forAll(mesh_.cells(), cellI)
423  {
424  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
425  state[cellI] = stateVecArray[localIdx];
426  }
427  }
428 
429  forAll(stateInfo_["modelStates"], idxI)
430  {
431  // lookup state from meshDb
432  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
433 
434  forAll(mesh_.cells(), cellI)
435  {
436  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
437  state[cellI] = stateVecArray[localIdx];
438  }
439  }
440 
441  forAll(stateInfo_["surfaceScalarStates"], idxI)
442  {
443  // lookup state from meshDb
444  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
445 
446  forAll(mesh_.faces(), faceI)
447  {
448  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
449  if (faceI < daIndex_.nLocalInternalFaces)
450  {
451  state[faceI] = stateVecArray[localIdx];
452  }
453  else
454  {
455  label relIdx = faceI - daIndex_.nLocalInternalFaces;
456  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
457  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
458  state.boundaryFieldRef()[patchIdx][faceIdx] = stateVecArray[localIdx];
459  }
460  }
461  }
462  VecRestoreArrayRead(stateVec, &stateVecArray);
463 }
464 
465 void DAField::point2OFMesh(const scalar* volCoords) const
466 {
467  /*
468  Description:
469  Assign the points in fvMesh of OpenFOAM based on the volCoords array
470 
471  Input:
472  volCoords: a vector that stores the x, y, and z coordinates for all
473  points in the fvMesh mesh
474 
475  Output:
476  New mesh metrics in fvMesh, effectively by calling mesh.movePoints
477 
478  Example:
479  Image we have three points in fvMesh, running on two CPU
480  processors, the proc0 owns one point and the proc1 owns two points,
481  then calling this function will assign xvVec based on the the points
482  coordinates in fvMesh
483 
484  xvVec = [x0, y0, z0 | x0, y0, z0, x1, y1, z1] <- x0 means x coordinate for the 0th point on local processor
485  0 1 2 | 3 4 5 6 7 8 <- global point vec index
486  --- proc0 --|--------- proc1 -------
487  */
488 
489  pointField meshPoints(mesh_.points());
490 
491  forAll(mesh_.points(), pointI)
492  {
493  for (label comp = 0; comp < 3; comp++)
494  {
495  label localIdx = daIndex_.getLocalXvIndex(pointI, comp);
496  meshPoints[pointI][comp] = volCoords[localIdx];
497  }
498  }
499 
500  // movePoints update the mesh metrics such as volume, surface area and cell centers
501  fvMesh& mesh = const_cast<fvMesh&>(mesh_);
502  mesh.movePoints(meshPoints);
503  mesh.moving(false);
504 }
505 
506 void DAField::pointVec2OFMesh(const Vec xvVec) const
507 {
508  /*
509  Description:
510  Assign the points in fvMesh of OpenFOAM based on the point vector
511 
512  Input:
513  xvVec: a vector that stores the x, y, and z coordinates for all
514  points in the fvMesh mesh
515 
516  Output:
517  New mesh metrics in fvMesh, effectively by calling mesh.movePoints
518 
519  Example:
520  Image we have three points in fvMesh, running on two CPU
521  processors, the proc0 owns one point and the proc1 owns two points,
522  then calling this function will assign xvVec based on the the points
523  coordinates in fvMesh
524 
525  xvVec = [x0, y0, z0 | x0, y0, z0, x1, y1, z1] <- x0 means x coordinate for the 0th point on local processor
526  0 1 2 | 3 4 5 6 7 8 <- global point vec index
527  --- proc0 --|--------- proc1 -------
528  */
529 
530  const PetscScalar* xvVecArray;
531  VecGetArrayRead(xvVec, &xvVecArray);
532 
533  pointField meshPoints(mesh_.points());
534 
535  forAll(mesh_.points(), pointI)
536  {
537  for (label comp = 0; comp < 3; comp++)
538  {
539  label localIdx = daIndex_.getLocalXvIndex(pointI, comp);
540  meshPoints[pointI][comp] = xvVecArray[localIdx];
541  }
542  }
543 
544  VecRestoreArrayRead(xvVec, &xvVecArray);
545 
546  // movePoints update the mesh metrics such as volume, surface area and cell centers
547  fvMesh& mesh = const_cast<fvMesh&>(mesh_);
548  mesh.movePoints(meshPoints);
549  mesh.moving(false);
550 }
551 
552 void DAField::ofMesh2PointVec(Vec xvVec) const
553 {
554  /*
555  Description:
556  Assign the point vector based on the points in fvMesh of OpenFOAM
557 
558  Input:
559  Mesh coordinates in fvMesh
560 
561  Output:
562  xvVec: a vector that stores the x, y, and z coordinates for all
563  points in the fvMesh mesh
564 
565  Example:
566  Image we have three points in fvMesh, running on two CPU
567  processors, the proc0 owns one point and the proc1 owns two points,
568  then calling this function will assign xvVec based on the the points
569  coordinates in fvMesh
570 
571  xvVec = [x0, y0, z0 | x0, y0, z0, x1, y1, z1] <- x0 means x coordinate for the 0th point on local processor
572  0 1 2 | 3 4 5 6 7 8 <- global point vec index
573  --- proc0 --|--------- proc1 -------
574  */
575 
576  PetscScalar* xvVecArray;
577  VecGetArray(xvVec, &xvVecArray);
578 
579  forAll(mesh_.points(), pointI)
580  {
581  for (label comp = 0; comp < 3; comp++)
582  {
583  label localIdx = daIndex_.getLocalXvIndex(pointI, comp);
584  assignValueCheckAD(xvVecArray[localIdx], mesh_.points()[pointI][comp]);
585  }
586  }
587 
588  VecRestoreArray(xvVec, &xvVecArray);
589 }
590 
591 void DAField::ofResField2ResVec(Vec resVec) const
592 {
593  /*
594  Description:
595  Assign values for the residual vector based on the
596  latest OpenFOAM residual field values
597 
598  Input:
599  OpenFOAM residual field variables
600 
601  Output:
602  resVec: state residual vector
603 
604  Example:
605  Image we have two state residuals (pRes,TRes) and five cells, running on two CPU
606  processors, the proc0 owns two cells and the proc1 owns three cells,
607  then calling this function gives the residual vector (state-by-state ordering):
608 
609  resVec = [pRes0, pRes1, TRes0, TRes1 | pRes0, pRes1, pRes2, TRes0, TRes1, TRes2]
610  0 1 2 3 | 4 5 6 7 8 9 <- global residual vec index
611  ---------- proc0 ---------|------------- proc1 ----------------------
612  NOTE: pRes0 means p residual for the 0th cell on local processor
613  */
614 
615  const objectRegistry& db = mesh_.thisDb();
616  PetscScalar* stateResVecArray;
617  VecGetArray(resVec, &stateResVecArray);
618 
619  forAll(stateInfo_["volVectorStates"], idxI)
620  {
621  // lookup state from meshDb
622  makeStateRes(stateInfo_["volVectorStates"][idxI], volVectorField, db);
623 
624  forAll(mesh_.cells(), cellI)
625  {
626  for (label comp = 0; comp < 3; comp++)
627  {
628  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
629  assignValueCheckAD(stateResVecArray[localIdx], stateRes[cellI][comp]);
630  }
631  }
632  }
633 
634  forAll(stateInfo_["volScalarStates"], idxI)
635  {
636  // lookup state from meshDb
637  makeStateRes(stateInfo_["volScalarStates"][idxI], volScalarField, db);
638 
639  forAll(mesh_.cells(), cellI)
640  {
641  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
642  assignValueCheckAD(stateResVecArray[localIdx], stateRes[cellI]);
643  }
644  }
645 
646  forAll(stateInfo_["modelStates"], idxI)
647  {
648  // lookup state from meshDb
649  makeStateRes(stateInfo_["modelStates"][idxI], volScalarField, db);
650 
651  forAll(mesh_.cells(), cellI)
652  {
653  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
654  assignValueCheckAD(stateResVecArray[localIdx], stateRes[cellI]);
655  }
656  }
657 
658  forAll(stateInfo_["surfaceScalarStates"], idxI)
659  {
660  // lookup state from meshDb
661  makeStateRes(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
662 
663  forAll(mesh_.faces(), faceI)
664  {
665  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
666  if (faceI < daIndex_.nLocalInternalFaces)
667  {
668  assignValueCheckAD(stateResVecArray[localIdx], stateRes[faceI]);
669  }
670  else
671  {
672  label relIdx = faceI - daIndex_.nLocalInternalFaces;
673  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
674  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
675  assignValueCheckAD(stateResVecArray[localIdx], stateRes.boundaryField()[patchIdx][faceIdx]);
676  }
677  }
678  }
679  VecRestoreArray(resVec, &stateResVecArray);
680 }
681 
682 void DAField::resVec2OFResField(const Vec resVec) const
683 {
684  /*
685  Description:
686  Assign OpenFOAM residual values based on the residual vector
687 
688  Input:
689  resVec: residual vector
690 
691  Output:
692  OpenFoam field variables
693 
694  Example:
695  Image we have two state residuals (pRes,TRes) and five cells, running on two CPU
696  processors, the proc0 owns two cells and the proc1 owns three cells,
697  then calling this function gives the residual vector (state-by-state ordering):
698 
699  resVec = [pRes0, pRes1, TRes0, TRes1 | pRes0, pRes1, pRes2, TRes0, TRes1, TRes2]
700  0 1 2 3 | 4 5 6 7 8 9 <- global residual vec index
701  ---------- proc0 ---------|------------- proc1 ----------------------
702  NOTE: pRes0 means p residual for the 0th cell on local processor
703  */
704 
705  const objectRegistry& db = mesh_.thisDb();
706  const PetscScalar* stateResVecArray;
707  VecGetArrayRead(resVec, &stateResVecArray);
708 
709  forAll(stateInfo_["volVectorStates"], idxI)
710  {
711  // lookup state from meshDb
712  makeStateRes(stateInfo_["volVectorStates"][idxI], volVectorField, db);
713 
714  forAll(mesh_.cells(), cellI)
715  {
716  for (label comp = 0; comp < 3; comp++)
717  {
718  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
719  stateRes[cellI][comp] = stateResVecArray[localIdx];
720  }
721  }
722  }
723 
724  forAll(stateInfo_["volScalarStates"], idxI)
725  {
726  // lookup state from meshDb
727  makeStateRes(stateInfo_["volScalarStates"][idxI], volScalarField, db);
728 
729  forAll(mesh_.cells(), cellI)
730  {
731  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
732  stateRes[cellI] = stateResVecArray[localIdx];
733  }
734  }
735 
736  forAll(stateInfo_["modelStates"], idxI)
737  {
738  // lookup state from meshDb
739  makeStateRes(stateInfo_["modelStates"][idxI], volScalarField, db);
740 
741  forAll(mesh_.cells(), cellI)
742  {
743  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
744  stateRes[cellI] = stateResVecArray[localIdx];
745  }
746  }
747 
748  forAll(stateInfo_["surfaceScalarStates"], idxI)
749  {
750  // lookup state from meshDb
751  makeStateRes(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
752 
753  forAll(mesh_.faces(), faceI)
754  {
755  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
756  if (faceI < daIndex_.nLocalInternalFaces)
757  {
758  stateRes[faceI] = stateResVecArray[localIdx];
759  }
760  else
761  {
762  label relIdx = faceI - daIndex_.nLocalInternalFaces;
763  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
764  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
765  stateRes.boundaryFieldRef()[patchIdx][faceIdx] = stateResVecArray[localIdx];
766  }
767  }
768  }
769  VecRestoreArrayRead(resVec, &stateResVecArray);
770 }
771 
772 void DAField::checkSpecialBCs()
773 {
774  /*
775  Description:
776  Check if we need to do special treatment for boundary conditions
777  If any special BC is detected, append their names to specialBCs list
778  */
779 
780  // *******************************************************************
781  // pressureInletVelocity
782  // *******************************************************************
783  // Note we need to read the U field, instead of getting it from db
784  // this is because coloringSolver does not read U
785  // Also we need to set read_if_present for solid solvers in which
786  // there is no U field
787  volVectorField U(
788  IOobject(
789  "U",
790  mesh_.time().timeName(),
791  mesh_,
792  IOobject::READ_IF_PRESENT,
793  IOobject::NO_WRITE,
794  false),
795  mesh_,
796  dimensionedVector("U", dimensionSet(0, 1, -1, 0, 0, 0, 0), vector::zero),
797  zeroGradientFvPatchField<vector>::typeName);
798  forAll(U.boundaryField(), patchI)
799  {
800  if (U.boundaryFieldRef()[patchI].type() == "pressureInletVelocity")
801  {
802  specialBCs.append("pressureInletVelocity");
803  break;
804  }
805  }
806 
807  // *******************************************************************
808  // append more special BCs
809  // *******************************************************************
810 }
811 
812 void DAField::specialBCTreatment()
813 {
814  /*
815  Description:
816  Apply special treatment for boundary conditions
817  */
818 
819  // *******************************************************************
820  // pressureInletVelocity
821  // *******************************************************************
822  // for pressureInletVelocity, the inlet U depends on
823  // rho and phi, so we need to call U.correctBoundaryConditions again
824  if (DAUtility::isInList<word>("pressureInletVelocity", specialBCs))
825  {
826 
827  volVectorField& U(const_cast<volVectorField&>(
828  mesh_.thisDb().lookupObject<volVectorField>("U")));
829  U.correctBoundaryConditions();
830  }
831 
832  // *******************************************************************
833  // append more special BCs
834  // *******************************************************************
835 }
836 
837 void DAField::setPrimalBoundaryConditions(const label printInfo)
838 {
839  /*
840  Description:
841  A general function to read the inlet/outlet values from DAOption, and set
842  the corresponding values to the boundary field. It also setup turbulence
843  wall boundary condition
844  Note: this function should be called before running the primal solver
845  If nothing is set, the BC will remain unchanged
846  Example
847  "primalBC":
848  {
849  "bc0":
850  {
851  "patches": ["inlet", "side"],
852  "variable": "U",
853  "value": [10.0 0.0, 0.0],
854  },
855  "bc1":
856  {
857  "patches": ["outlet"],
858  "variable": "p",
859  "value": [101325.0],
860  },
861  "useWallFunction": True,
862  "MRF": 1000.0,
863  "fvSource": {"value": 0.01, "comp": 0}
864  }
865  */
866 
867  const objectRegistry& db = mesh_.thisDb();
868 
869  label setTurbWallBCs = 0;
870  label useWallFunction = 0;
871 
872  const dictionary& allOptions = daOption_.getAllOptions();
873  dictionary bcDict = allOptions.subDict("primalBC");
874 
875  forAll(bcDict.toc(), idxI)
876  {
877  word bcKey = bcDict.toc()[idxI];
878 
879  if (bcKey == "useWallFunction")
880  {
881  setTurbWallBCs = 1;
882  useWallFunction = bcDict.getLabel("useWallFunction");
883  continue;
884  }
885  else if (bcKey == "MRF")
886  {
887  // change the rotation speed in MRF
888  scalar omegaNew = bcDict.getScalar("MRF");
889  const IOMRFZoneListDF& MRF = db.lookupObject<IOMRFZoneListDF>("MRFProperties");
890  scalar& omega = const_cast<scalar&>(MRF.getOmegaRef());
891  omega = omegaNew;
892 
893  if (printInfo)
894  {
895  Info << "Setting MRF omega to " << omegaNew << endl;
896  }
897 
898  continue;
899  }
900  else if (bcKey == "fvSource")
901  {
902  // assign the fvSource field with an uniform value
903  dictionary subDict = bcDict.subDict("fvSource");
904  scalar fvSourceVal = subDict.getScalar("value");
905  label compI = subDict.getLabel("comp");
906 
907  volVectorField& fvSource = const_cast<volVectorField&>(
908  db.lookupObject<volVectorField>("fvSource"));
909 
910  forAll(fvSource, cellI)
911  {
912  fvSource[cellI][compI] = fvSourceVal;
913  }
914  forAll(fvSource.boundaryField(), patchI)
915  {
916  forAll(fvSource.boundaryField()[patchI], faceI)
917  {
918  fvSource.boundaryFieldRef()[patchI][faceI][compI] = fvSourceVal;
919  }
920  }
921  fvSource.correctBoundaryConditions();
922 
923  if (printInfo)
924  {
925  Info << "Setting fvSource to " << fvSourceVal << endl;
926  }
927 
928  continue;
929  }
930  else if (bcKey == "transport:nu")
931  {
932  // change the nu field
933  scalar nu = bcDict.getScalar("transport:nu");
934  volScalarField& nuField = const_cast<volScalarField&>(
935  db.lookupObject<volScalarField>("nu"));
936  forAll(nuField, cellI)
937  {
938  nuField[cellI] = nu;
939  }
940  forAll(nuField.boundaryField(), patchI)
941  {
942  forAll(nuField.boundaryField()[patchI], faceI)
943  {
944  nuField.boundaryFieldRef()[patchI][faceI] = nu;
945  }
946  }
947  nuField.correctBoundaryConditions();
948 
949  if (printInfo)
950  {
951  Info << "Setting transportProperties nu to " << nu << endl;
952  }
953 
954  continue;
955  }
956  else if (bcKey == "thermo:mu")
957  {
958  // change the nu field
959  scalar mu = bcDict.getScalar("thermo:mu");
960  volScalarField& muField = const_cast<volScalarField&>(
961  db.lookupObject<volScalarField>("thermo:mu"));
962  forAll(muField, cellI)
963  {
964  muField[cellI] = mu;
965  }
966  forAll(muField.boundaryField(), patchI)
967  {
968  forAll(muField.boundaryField()[patchI], faceI)
969  {
970  muField.boundaryFieldRef()[patchI][faceI] = mu;
971  }
972  }
973  muField.correctBoundaryConditions();
974 
975  if (printInfo)
976  {
977  Info << "Setting thermalphysicalProperties mu to " << mu << endl;
978  }
979 
980  continue;
981  }
982 
983  dictionary bcSubDict = bcDict.subDict(bcKey);
984 
985  wordList patches;
986  bcSubDict.readEntry<wordList>("patches", patches);
987  word variable = bcSubDict.getWord("variable");
988  scalarList value;
989  bcSubDict.readEntry<scalarList>("value", value);
990 
991  // loop over all patches and set values
992  forAll(patches, idxI)
993  {
994  word patch = patches[idxI];
995 
996  // it should be a scalar
997  if (value.size() == 1)
998  {
999  if (!db.foundObject<volScalarField>(variable))
1000  {
1001  if (printInfo)
1002  {
1003  Info << variable << " not found, skip it." << endl;
1004  }
1005  continue;
1006  }
1007  // it is a scalar
1008  volScalarField& state(const_cast<volScalarField&>(
1009  db.lookupObject<volScalarField>(variable)));
1010 
1011  if (printInfo)
1012  {
1013  Info << "Setting primal boundary conditions..." << endl;
1014  Info << "Setting " << variable << " = " << value[0] << " at " << patch << endl;
1015  }
1016 
1017  label patchI = mesh_.boundaryMesh().findPatchID(patch);
1018 
1019  // for decomposed domain, don't set BC if the patch is empty
1020  if (mesh_.boundaryMesh()[patchI].size() > 0)
1021  {
1022  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
1023  {
1024  forAll(state.boundaryFieldRef()[patchI], faceI)
1025  {
1026  state.boundaryFieldRef()[patchI][faceI] = value[0];
1027  }
1028  }
1029  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
1030  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
1031  {
1032  // set value
1033  forAll(state.boundaryFieldRef()[patchI], faceI)
1034  {
1035  state.boundaryFieldRef()[patchI][faceI] = value[0];
1036  }
1037  // set inletValue
1038  mixedFvPatchField<scalar>& inletOutletPatch =
1039  refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
1040 
1041  inletOutletPatch.refValue() = value[0];
1042  }
1043  else if (state.boundaryFieldRef()[patchI].type() == "fixedGradient")
1044  {
1045  fixedGradientFvPatchField<scalar>& patchBC =
1046  refCast<fixedGradientFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
1047  scalarField& grad = const_cast<scalarField&>(patchBC.gradient());
1048  forAll(grad, idxI)
1049  {
1050  grad[idxI] = value[0];
1051  }
1052  }
1053  else
1054  {
1055  FatalErrorIn("") << "only support fixedValues, inletOutlet, "
1056  << "outletInlet, fixedGradient!" << abort(FatalError);
1057  }
1058  }
1059  }
1060  else if (value.size() == 3)
1061  {
1062  if (!db.foundObject<volVectorField>(variable))
1063  {
1064  if (printInfo)
1065  {
1066  Info << variable << " not found, skip it." << endl;
1067  }
1068  continue;
1069  }
1070  // it is a vector
1071  volVectorField& state(const_cast<volVectorField&>(
1072  db.lookupObject<volVectorField>(variable)));
1073 
1074  vector valVec = {value[0], value[1], value[2]};
1075  if (printInfo)
1076  {
1077  Info << "Setting primal boundary conditions..." << endl;
1078  Info << "Setting " << variable << " = (" << value[0] << " "
1079  << value[1] << " " << value[2] << ") at " << patch << endl;
1080  }
1081 
1082  label patchI = mesh_.boundaryMesh().findPatchID(patch);
1083 
1084  // for decomposed domain, don't set BC if the patch is empty
1085  if (mesh_.boundaryMesh()[patchI].size() > 0)
1086  {
1087  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
1088  {
1089  forAll(state.boundaryFieldRef()[patchI], faceI)
1090  {
1091  state.boundaryFieldRef()[patchI][faceI] = valVec;
1092  }
1093  }
1094  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
1095  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
1096  {
1097  // set value
1098  forAll(state.boundaryFieldRef()[patchI], faceI)
1099  {
1100  state.boundaryFieldRef()[patchI][faceI] = valVec;
1101  }
1102  // set inletValue
1103  mixedFvPatchField<vector>& inletOutletPatch =
1104  refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
1105  inletOutletPatch.refValue() = valVec;
1106  }
1107  else if (state.boundaryFieldRef()[patchI].type() == "fixedGradient")
1108  {
1109  fixedGradientFvPatchField<vector>& patchBC =
1110  refCast<fixedGradientFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
1111  vectorField& grad = const_cast<vectorField&>(patchBC.gradient());
1112  forAll(grad, idxI)
1113  {
1114  grad[idxI][0] = value[0];
1115  grad[idxI][1] = value[1];
1116  grad[idxI][2] = value[2];
1117  }
1118  }
1119  else
1120  {
1121  FatalErrorIn("") << "only support fixedValues, inletOutlet, "
1122  << "fixedGradient, and tractionDisplacement!"
1123  << abort(FatalError);
1124  }
1125  }
1126  }
1127  else
1128  {
1129  FatalErrorIn("") << "value should be a list of either 1 (scalar) "
1130  << "or 3 (vector) elements" << abort(FatalError);
1131  }
1132  }
1133  }
1134 
1135  // we also set wall boundary conditions for turbulence variables
1136  if (setTurbWallBCs)
1137  {
1138 
1139  IOdictionary turbDict(
1140  IOobject(
1141  "turbulenceProperties",
1142  mesh_.time().constant(),
1143  mesh_,
1144  IOobject::MUST_READ,
1145  IOobject::NO_WRITE,
1146  false));
1147 
1148  dictionary coeffDict(turbDict.subDict("RAS"));
1149  word turbModelType = word(coeffDict["RASModel"]);
1150  //Info<<"turbModelType: "<<turbModelType<<endl;
1151 
1152  // create word regular expression for SA model
1153  wordReList SAModelWordReList {
1154  {"SpalartAllmaras.*", wordRe::REGEX}};
1155  wordRes SAModelWordRes(SAModelWordReList);
1156 
1157  // ------ nut ----------
1158  if (db.foundObject<volScalarField>("nut"))
1159  {
1160 
1161  volScalarField& nut(const_cast<volScalarField&>(
1162  db.lookupObject<volScalarField>("nut")));
1163 
1164  forAll(nut.boundaryField(), patchI)
1165  {
1166  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1167  {
1168  if (printInfo)
1169  {
1170  Info << "Setting nut wall BC for "
1171  << mesh_.boundaryMesh()[patchI].name() << ". ";
1172  }
1173 
1174  if (useWallFunction)
1175  {
1176  // wall function for SA
1177  if (SAModelWordRes(turbModelType))
1178  {
1179  nut.boundaryFieldRef().set(
1180  patchI,
1181  fvPatchField<scalar>::New(
1182  "nutUSpaldingWallFunction",
1183  mesh_.boundary()[patchI],
1184  nut));
1185 
1186  if (printInfo)
1187  {
1188  Info << "BCType=nutUSpaldingWallFunction" << endl;
1189  }
1190  }
1191  else // wall function for kOmega and kEpsilon
1192  {
1193  nut.boundaryFieldRef().set(
1194  patchI,
1195  fvPatchField<scalar>::New(
1196  "nutkWallFunction",
1197  mesh_.boundary()[patchI],
1198  nut));
1199 
1200  if (printInfo)
1201  {
1202  Info << "BCType=nutkWallFunction" << endl;
1203  }
1204  }
1205 
1206  // set boundary values
1207  // for decomposed domain, don't set BC if the patch is empty
1208  if (mesh_.boundaryMesh()[patchI].size() > 0)
1209  {
1210  scalar wallVal = nut[0];
1211  forAll(nut.boundaryFieldRef()[patchI], faceI)
1212  {
1213  // assign uniform field
1214  nut.boundaryFieldRef()[patchI][faceI] = wallVal;
1215  }
1216  }
1217  }
1218  else
1219  {
1220  nut.boundaryFieldRef().set(
1221  patchI,
1222  fvPatchField<scalar>::New(
1223  "nutLowReWallFunction",
1224  mesh_.boundary()[patchI],
1225  nut));
1226 
1227  if (printInfo)
1228  {
1229  Info << "BCType=nutLowReWallFunction" << endl;
1230  }
1231 
1232  // set boundary values
1233  // for decomposed domain, don't set BC if the patch is empty
1234  if (mesh_.boundaryMesh()[patchI].size() > 0)
1235  {
1236  forAll(nut.boundaryFieldRef()[patchI], faceI)
1237  {
1238  // assign uniform field
1239  nut.boundaryFieldRef()[patchI][faceI] = 1e-14;
1240  }
1241  }
1242  }
1243  }
1244  }
1245  }
1246 
1247  // ------ k ----------
1248  if (db.foundObject<volScalarField>("k"))
1249  {
1250 
1251  volScalarField& k(const_cast<volScalarField&>(
1252  db.lookupObject<volScalarField>("k")));
1253 
1254  forAll(k.boundaryField(), patchI)
1255  {
1256  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1257  {
1258  if (printInfo)
1259  {
1260  Info << "Setting k wall BC for "
1261  << mesh_.boundaryMesh()[patchI].name() << ". ";
1262  }
1263 
1264  if (useWallFunction)
1265  {
1266  // wall function for SA
1267  k.boundaryFieldRef().set(
1268  patchI,
1269  fvPatchField<scalar>::New("kqRWallFunction", mesh_.boundary()[patchI], k));
1270 
1271  if (printInfo)
1272  {
1273  Info << "BCType=kqRWallFunction" << endl;
1274  }
1275 
1276  // set boundary values
1277  // for decomposed domain, don't set BC if the patch is empty
1278  if (mesh_.boundaryMesh()[patchI].size() > 0)
1279  {
1280  scalar wallVal = k[0];
1281  forAll(k.boundaryFieldRef()[patchI], faceI)
1282  {
1283  k.boundaryFieldRef()[patchI][faceI] = wallVal; // assign uniform field
1284  }
1285  }
1286  }
1287  else
1288  {
1289  k.boundaryFieldRef().set(
1290  patchI,
1291  fvPatchField<scalar>::New("fixedValue", mesh_.boundary()[patchI], k));
1292 
1293  if (printInfo)
1294  {
1295  Info << "BCType=fixedValue" << endl;
1296  }
1297 
1298  // set boundary values
1299  // for decomposed domain, don't set BC if the patch is empty
1300  if (mesh_.boundaryMesh()[patchI].size() > 0)
1301  {
1302  forAll(k.boundaryFieldRef()[patchI], faceI)
1303  {
1304  k.boundaryFieldRef()[patchI][faceI] = 1e-14; // assign uniform field
1305  }
1306  }
1307  }
1308  }
1309  }
1310  }
1311 
1312  // ------ omega ----------
1313  if (db.foundObject<volScalarField>("omega"))
1314  {
1315 
1316  volScalarField& omega(const_cast<volScalarField&>(
1317  db.lookupObject<volScalarField>("omega")));
1318 
1319  forAll(omega.boundaryField(), patchI)
1320  {
1321  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1322  {
1323  if (printInfo)
1324  {
1325  Info << "Setting omega wall BC for "
1326  << mesh_.boundaryMesh()[patchI].name() << ". ";
1327  }
1328 
1329  // always use omegaWallFunction
1330  omega.boundaryFieldRef().set(
1331  patchI,
1332  fvPatchField<scalar>::New("omegaWallFunction", mesh_.boundary()[patchI], omega));
1333 
1334  if (printInfo)
1335  {
1336  Info << "BCType=omegaWallFunction" << endl;
1337  }
1338 
1339  // set boundary values
1340  // for decomposed domain, don't set BC if the patch is empty
1341  if (mesh_.boundaryMesh()[patchI].size() > 0)
1342  {
1343  scalar wallVal = omega[0];
1344  forAll(omega.boundaryFieldRef()[patchI], faceI)
1345  {
1346  omega.boundaryFieldRef()[patchI][faceI] = wallVal; // assign uniform field
1347  }
1348  }
1349  }
1350  }
1351  }
1352 
1353  // ------ epsilon ----------
1354  if (db.foundObject<volScalarField>("epsilon"))
1355  {
1356 
1357  volScalarField& epsilon(const_cast<volScalarField&>(
1358  db.lookupObject<volScalarField>("epsilon")));
1359 
1360  forAll(epsilon.boundaryField(), patchI)
1361  {
1362  if (mesh_.boundaryMesh()[patchI].type() == "wall")
1363  {
1364 
1365  if (printInfo)
1366  {
1367  Info << "Setting epsilon wall BC for "
1368  << mesh_.boundaryMesh()[patchI].name() << ". ";
1369  }
1370 
1371  if (useWallFunction)
1372  {
1373  epsilon.boundaryFieldRef().set(
1374  patchI,
1375  fvPatchField<scalar>::New("epsilonWallFunction", mesh_.boundary()[patchI], epsilon));
1376 
1377  if (printInfo)
1378  {
1379  Info << "BCType=epsilonWallFunction" << endl;
1380  }
1381 
1382  // set boundary values
1383  // for decomposed domain, don't set BC if the patch is empty
1384  if (mesh_.boundaryMesh()[patchI].size() > 0)
1385  {
1386  scalar wallVal = epsilon[0];
1387  forAll(epsilon.boundaryFieldRef()[patchI], faceI)
1388  {
1389  epsilon.boundaryFieldRef()[patchI][faceI] = wallVal; // assign uniform field
1390  }
1391  }
1392  }
1393  else
1394  {
1395  epsilon.boundaryFieldRef().set(
1396  patchI,
1397  fvPatchField<scalar>::New("fixedValue", mesh_.boundary()[patchI], epsilon));
1398 
1399  if (printInfo)
1400  {
1401  Info << "BCType=fixedValue" << endl;
1402  }
1403 
1404  // set boundary values
1405  // for decomposed domain, don't set BC if the patch is empty
1406  if (mesh_.boundaryMesh()[patchI].size() > 0)
1407  {
1408  forAll(epsilon.boundaryFieldRef()[patchI], faceI)
1409  {
1410  epsilon.boundaryFieldRef()[patchI][faceI] = 1e-14; // assign uniform field
1411  }
1412  }
1413  }
1414  }
1415  }
1416  }
1417  }
1418 }
1419 
1420 void DAField::ofField2List(
1421  scalarList& stateList,
1422  scalarList& stateBoundaryList) const
1423 {
1424  /*
1425  Description:
1426  Assign values for the scalar list of states based on the latest OpenFOAM field values.
1427  This function is similar to DAField::ofField2StateVec except that the output are
1428  scalarLists instead of a Petsc vector
1429 
1430  Input:
1431  OpenFOAM field variables
1432 
1433  Output:
1434  stateList: scalar list of states
1435  stateBoundaryList: scalar list of boundary states
1436 
1437  Example:
1438  Image we have two state variables (p,T) and five cells, running on two CPU
1439  processors, the proc0 owns two cells and the proc1 owns three cells,
1440  then calling this function gives the scalar list of states (state-by-state ordering):
1441 
1442  scalarList = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
1443  0 1 2 3 | 4 5 6 7 8 9 <- global state index
1444  ---- proc0 -----|--------- proc1 -------
1445  */
1446 
1447  const objectRegistry& db = mesh_.thisDb();
1448 
1449  label localBFaceI = 0;
1450 
1451  forAll(stateInfo_["volVectorStates"], idxI)
1452  {
1453  // lookup state from meshDb
1454  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
1455 
1456  forAll(mesh_.cells(), cellI)
1457  {
1458  for (label comp = 0; comp < 3; comp++)
1459  {
1460  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
1461  stateList[localIdx] = state[cellI][comp];
1462  }
1463  }
1464 
1465  forAll(state.boundaryField(), patchI)
1466  {
1467  if (state.boundaryField()[patchI].size() > 0)
1468  {
1469  forAll(state.boundaryField()[patchI], faceI)
1470  {
1471  for (label comp = 0; comp < 3; comp++)
1472  {
1473  stateBoundaryList[localBFaceI] = state.boundaryField()[patchI][faceI][comp];
1474  localBFaceI++;
1475  }
1476  }
1477  }
1478  }
1479  }
1480 
1481  forAll(stateInfo_["volScalarStates"], idxI)
1482  {
1483  // lookup state from meshDb
1484  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
1485 
1486  forAll(mesh_.cells(), cellI)
1487  {
1488  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
1489  stateList[localIdx] = state[cellI];
1490  }
1491 
1492  forAll(state.boundaryField(), patchI)
1493  {
1494  if (state.boundaryField()[patchI].size() > 0)
1495  {
1496  forAll(state.boundaryField()[patchI], faceI)
1497  {
1498  stateBoundaryList[localBFaceI] = state.boundaryField()[patchI][faceI];
1499  localBFaceI++;
1500  }
1501  }
1502  }
1503  }
1504 
1505  forAll(stateInfo_["modelStates"], idxI)
1506  {
1507  // lookup state from meshDb
1508  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
1509 
1510  forAll(mesh_.cells(), cellI)
1511  {
1512  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
1513  stateList[localIdx] = state[cellI];
1514  }
1515 
1516  forAll(state.boundaryField(), patchI)
1517  {
1518  if (state.boundaryField()[patchI].size() > 0)
1519  {
1520  forAll(state.boundaryField()[patchI], faceI)
1521  {
1522  stateBoundaryList[localBFaceI] = state.boundaryField()[patchI][faceI];
1523  localBFaceI++;
1524  }
1525  }
1526  }
1527  }
1528 
1529  forAll(stateInfo_["surfaceScalarStates"], idxI)
1530  {
1531  // lookup state from meshDb
1532  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
1533 
1534  forAll(mesh_.faces(), faceI)
1535  {
1536  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
1537  if (faceI < daIndex_.nLocalInternalFaces)
1538  {
1539  stateList[localIdx] = state[faceI];
1540  }
1541  else
1542  {
1543  label relIdx = faceI - daIndex_.nLocalInternalFaces;
1544  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
1545  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
1546  stateList[localIdx] = state.boundaryField()[patchIdx][faceIdx];
1547  }
1548  }
1549  }
1550 }
1551 
1552 void DAField::list2OFField(
1553  const scalarList& stateList,
1554  const scalarList& stateBoundaryList,
1555  const label oldTimeLevel) const
1556 {
1557  /*
1558  Description:
1559  Assign values OpenFOAM field values based on the scalar list of states
1560 
1561  Input:
1562  stateList: scalar list of states
1563  stateBoundaryList: scalar list of boundary states
1564  oldTimeLevel: assign to oldTime field instead of the original field, this will
1565  be used in time-accurate adjoint
1566 
1567  Output:
1568  OpenFoam field variables
1569 
1570  Example:
1571  Image we have two state variables (p,T) and five cells, running on two CPU
1572  processors, the proc0 owns two cells and the proc1 owns three cells,
1573  then calling this function gives the scalar list of states (state-by-state ordering):
1574 
1575  scalarList = [p0, p1, T0, T1 | p0, p1, p2, T0, T1, T2] <- p0 means p for the 0th cell on local processor
1576  0 1 2 3 | 4 5 6 7 8 9 <- global state index
1577  ---- proc0 -----|--------- proc1 -------
1578  */
1579 
1580  const objectRegistry& db = mesh_.thisDb();
1581 
1582  label localBFaceI = 0;
1583 
1584  if (oldTimeLevel > 2 || oldTimeLevel < 0)
1585  {
1586  FatalErrorIn("") << "oldTimeLevel not valid!"
1587  << abort(FatalError);
1588  }
1589 
1590  forAll(stateInfo_["volVectorStates"], idxI)
1591  {
1592  // lookup state from meshDb
1593  makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db);
1594 
1595  label maxOldTimes = state.nOldTimes();
1596 
1597  if (maxOldTimes >= oldTimeLevel)
1598  {
1599  forAll(mesh_.cells(), cellI)
1600  {
1601  for (label comp = 0; comp < 3; comp++)
1602  {
1603  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp);
1604  if (oldTimeLevel == 0)
1605  {
1606  state[cellI][comp] = stateList[localIdx];
1607  }
1608  else if (oldTimeLevel == 1)
1609  {
1610  state.oldTime()[cellI][comp] = stateList[localIdx];
1611  }
1612  else if (oldTimeLevel == 2)
1613  {
1614  state.oldTime().oldTime()[cellI][comp] = stateList[localIdx];
1615  }
1616  }
1617  }
1618 
1619  forAll(state.boundaryField(), patchI)
1620  {
1621  if (state.boundaryField()[patchI].size() > 0)
1622  {
1623  forAll(state.boundaryField()[patchI], faceI)
1624  {
1625  for (label comp = 0; comp < 3; comp++)
1626  {
1627  if (oldTimeLevel == 0)
1628  {
1629  state.boundaryFieldRef()[patchI][faceI][comp] =
1630  stateBoundaryList[localBFaceI];
1631  }
1632  else if (oldTimeLevel == 1)
1633  {
1634  state.oldTime().boundaryFieldRef()[patchI][faceI][comp] =
1635  stateBoundaryList[localBFaceI];
1636  }
1637  else if (oldTimeLevel == 2)
1638  {
1639  state.oldTime().oldTime().boundaryFieldRef()[patchI][faceI][comp] =
1640  stateBoundaryList[localBFaceI];
1641  }
1642  localBFaceI++;
1643  }
1644  }
1645  }
1646  }
1647  }
1648  }
1649 
1650  forAll(stateInfo_["volScalarStates"], idxI)
1651  {
1652  // lookup state from meshDb
1653  makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db);
1654 
1655  label maxOldTimes = state.nOldTimes();
1656 
1657  if (maxOldTimes >= oldTimeLevel)
1658  {
1659 
1660  forAll(mesh_.cells(), cellI)
1661  {
1662  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
1663  if (oldTimeLevel == 0)
1664  {
1665  state[cellI] = stateList[localIdx];
1666  }
1667  else if (oldTimeLevel == 1)
1668  {
1669  state.oldTime()[cellI] = stateList[localIdx];
1670  }
1671  else if (oldTimeLevel == 2)
1672  {
1673  state.oldTime().oldTime()[cellI] = stateList[localIdx];
1674  }
1675  }
1676 
1677  forAll(state.boundaryField(), patchI)
1678  {
1679  if (state.boundaryField()[patchI].size() > 0)
1680  {
1681  forAll(state.boundaryField()[patchI], faceI)
1682  {
1683  if (oldTimeLevel == 0)
1684  {
1685  state.boundaryFieldRef()[patchI][faceI] =
1686  stateBoundaryList[localBFaceI];
1687  }
1688  else if (oldTimeLevel == 1)
1689  {
1690  state.oldTime().boundaryFieldRef()[patchI][faceI] =
1691  stateBoundaryList[localBFaceI];
1692  }
1693  else if (oldTimeLevel == 2)
1694  {
1695  state.oldTime().oldTime().boundaryFieldRef()[patchI][faceI] =
1696  stateBoundaryList[localBFaceI];
1697  }
1698  localBFaceI++;
1699  }
1700  }
1701  }
1702  }
1703  }
1704 
1705  forAll(stateInfo_["modelStates"], idxI)
1706  {
1707  // lookup state from meshDb
1708  makeState(stateInfo_["modelStates"][idxI], volScalarField, db);
1709 
1710  label maxOldTimes = state.nOldTimes();
1711 
1712  if (maxOldTimes >= oldTimeLevel)
1713  {
1714 
1715  forAll(mesh_.cells(), cellI)
1716  {
1717  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI);
1718  if (oldTimeLevel == 0)
1719  {
1720  state[cellI] = stateList[localIdx];
1721  }
1722  else if (oldTimeLevel == 1)
1723  {
1724  state.oldTime()[cellI] = stateList[localIdx];
1725  }
1726  else if (oldTimeLevel == 2)
1727  {
1728  state.oldTime().oldTime()[cellI] = stateList[localIdx];
1729  }
1730  }
1731 
1732  forAll(state.boundaryField(), patchI)
1733  {
1734  if (state.boundaryField()[patchI].size() > 0)
1735  {
1736  forAll(state.boundaryField()[patchI], faceI)
1737  {
1738  if (oldTimeLevel == 0)
1739  {
1740  state.boundaryFieldRef()[patchI][faceI] =
1741  stateBoundaryList[localBFaceI];
1742  }
1743  else if (oldTimeLevel == 1)
1744  {
1745  state.oldTime().boundaryFieldRef()[patchI][faceI] =
1746  stateBoundaryList[localBFaceI];
1747  }
1748  else if (oldTimeLevel == 2)
1749  {
1750  state.oldTime().oldTime().boundaryFieldRef()[patchI][faceI] =
1751  stateBoundaryList[localBFaceI];
1752  }
1753  localBFaceI++;
1754  }
1755  }
1756  }
1757  }
1758  }
1759 
1760  forAll(stateInfo_["surfaceScalarStates"], idxI)
1761  {
1762  // lookup state from meshDb
1763  makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db);
1764 
1765  label maxOldTimes = state.nOldTimes();
1766 
1767  if (maxOldTimes >= oldTimeLevel)
1768  {
1769 
1770  forAll(mesh_.faces(), faceI)
1771  {
1772  label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI);
1773  if (faceI < daIndex_.nLocalInternalFaces)
1774  {
1775 
1776  if (oldTimeLevel == 0)
1777  {
1778  state[faceI] = stateList[localIdx];
1779  }
1780  else if (oldTimeLevel == 1)
1781  {
1782  state.oldTime()[faceI] = stateList[localIdx];
1783  }
1784  else if (oldTimeLevel == 2)
1785  {
1786  state.oldTime().oldTime()[faceI] = stateList[localIdx];
1787  }
1788  }
1789  else
1790  {
1791  label relIdx = faceI - daIndex_.nLocalInternalFaces;
1792  const label& patchIdx = daIndex_.bFacePatchI[relIdx];
1793  const label& faceIdx = daIndex_.bFaceFaceI[relIdx];
1794  if (oldTimeLevel == 0)
1795  {
1796  state.boundaryFieldRef()[patchIdx][faceIdx] =
1797  stateList[localIdx];
1798  }
1799  else if (oldTimeLevel == 1)
1800  {
1801  state.oldTime().boundaryFieldRef()[patchIdx][faceIdx] =
1802  stateList[localIdx];
1803  }
1804  else if (oldTimeLevel == 2)
1805  {
1806  state.oldTime().oldTime().boundaryFieldRef()[patchIdx][faceIdx] =
1807  stateList[localIdx];
1808  }
1809  }
1810  }
1811  }
1812  }
1813 }
1814 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1815 
1816 } // End namespace Foam
1817 
1818 // ************************************************************************* //
mu
volScalarField & mu
Definition: createRefsSolidDisplacement.H:6
Foam::DAField::list2OFField
void list2OFField(const scalarList &stateList, const scalarList &stateBounaryList, const label oldTimeLevel) const
assign the fields in OpenFOAM based on the scalar list of states
Definition: DAField.C:1552
allOptions
const dictionary & allOptions
Definition: createRefsRhoSimpleC.H:15
U
U
Definition: pEqnPimpleDyM.H:60
Foam::DAField::daIndex_
const DAIndex & daIndex_
DAIndex object.
Definition: DAField.H:57
forAll
forAll(pseudoP.boundaryField(), patchI)
Definition: solvePseudoPEqn.H:10
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::setPrimalBoundaryConditions
void setPrimalBoundaryConditions(const label printInfo=1)
set the boundary conditions based on parameters defined in DAOption
Definition: DAField.C:837
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:288
daOption
DAOption daOption(mesh, pyOptions_)
Foam::DAField::specialBCs
wordList specialBCs
a list that contains the names of detected special boundary conditions
Definition: DAField.H:125
Foam::DAField::checkSpecialBCs
void checkSpecialBCs()
check if we need to do special treatment for boundary conditions
Definition: DAField.C:772
MRF
IOMRFZoneListDF & MRF
Definition: createRefsRhoSimple.H:18
Foam::DAField::resetOFSeeds
void resetOFSeeds()
reset the seeds for both state and mesh variables in OpenFOAM by setting their gradient() to zeros
Definition: DAField.C:36
Foam::DAField::mesh_
const fvMesh & mesh_
Foam::fvMesh object.
Definition: DAField.H:48
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
solverName
word solverName
Definition: createAdjointCompressible.H:30
Foam::DAField::ofField2StateVec
void ofField2StateVec(Vec stateVec) const
set the state vector based on the latest fields in OpenFOAM
Definition: DAField.C:198
Foam::DAField::ofMesh2PointVec
void ofMesh2PointVec(Vec xvVec) const
assign the point vector based on the points in fvMesh of OpenFOAM
Definition: DAField.C:552
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:465
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:59
k
dimensionedScalar & k
Definition: createRefsHeatTransfer.H:6
Foam
Definition: multiFreqScalarFvPatchField.C:144
Foam::DAField::ofField2List
void ofField2List(scalarList &stateList, scalarList &stateBoundaryList) const
set the scalar list of states based on the latest fields in OpenFOAM
Definition: DAField.C:1420
Foam::DAField::stateVec2OFField
void stateVec2OFField(const Vec stateVec) const
assign the fields in OpenFOAM based on the state vector
Definition: DAField.C:375
Foam::IOMRFZoneListDF
Definition: IOMRFZoneListDF.H:49
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
Foam::DAField::resVec2OFResField
void resVec2OFResField(const Vec resVec) const
assign the residual field in OpenFOAM based on the residual vector
Definition: DAField.C:682
Foam::DAField::ofResField2ResVec
void ofResField2ResVec(Vec resVec) const
assign the residual vector based on the residual field in OpenFOAM
Definition: DAField.C:591
daModel
DAModel daModel(mesh, daOption)
Foam::DAField::pointVec2OFMesh
void pointVec2OFMesh(const Vec xvVec) const
assign the points in fvMesh of OpenFOAM based on the point vector
Definition: DAField.C:506
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
daIndex
DAIndex daIndex(mesh, daOption, daModel)
DAField.H
Foam::DAField::specialBCTreatment
void specialBCTreatment()
apply special treatment for boundary conditions
Definition: DAField.C:812
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