DAJacCon.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v3
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAJacCon.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
16 
17 defineTypeNameAndDebug(DAJacCon, 0);
18 defineRunTimeSelectionTable(DAJacCon, dictionary);
19 
20 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
21 
22 DAJacCon::DAJacCon(
23  const word modelType,
24  const fvMesh& mesh,
25  const DAOption& daOption,
26  const DAModel& daModel,
27  const DAIndex& daIndex)
28  : modelType_(modelType),
29  mesh_(mesh),
30  daOption_(daOption),
31  daModel_(daModel),
32  daIndex_(daIndex),
33  daColoring_(mesh, daOption, daModel, daIndex),
34  daField_(mesh, daOption, daModel, daIndex)
35 {
36  // initialize stateInfo_
37  word solverName = daOption.getOption<word>("solverName");
38  autoPtr<DAStateInfo> daStateInfo(DAStateInfo::New(solverName, mesh, daOption, daModel));
39  stateInfo_ = daStateInfo->getStateInfo();
40 
41  // check if there is special boundary conditions that need special treatment in jacCon_
42  this->checkSpecialBCs();
43 }
44 
45 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
46 
47 autoPtr<DAJacCon> DAJacCon::New(
48  const word modelType,
49  const fvMesh& mesh,
50  const DAOption& daOption,
51  const DAModel& daModel,
52  const DAIndex& daIndex)
53 {
54  if (daOption.getAllOptions().lookupOrDefault<label>("debug", 0))
55  {
56  Info << "Selecting " << modelType << " for DAJacCon" << endl;
57  }
58 
59  dictionaryConstructorTable::iterator cstrIter =
60  dictionaryConstructorTablePtr_->find(modelType);
61 
62  // if the solver name is not found in any child class, print an error
63  if (cstrIter == dictionaryConstructorTablePtr_->end())
64  {
65  FatalErrorIn(
66  "DAJacCon::New"
67  "("
68  " const word,"
69  " const fvMesh&,"
70  " const DAOption&,"
71  " const DAModel&,"
72  " const DAIndex&"
73  ")")
74  << "Unknown DAJacCon type "
75  << modelType << nl << nl
76  << "Valid DAJacCon types:" << endl
77  << dictionaryConstructorTablePtr_->sortedToc()
78  << exit(FatalError);
79  }
80 
81  // child class found
82  return autoPtr<DAJacCon>(
83  cstrIter()(modelType, mesh, daOption, daModel, daIndex));
84 }
85 
86 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
87 
88 void DAJacCon::initializeStateBoundaryCon()
89 {
90  /*
91  Description:
92  Initialize state boundary connectivity matrices and variables.
93  This will be used for Jacobian with respect to states, e.g., dRdW, and dFdW.
94 
95  NOTE: no need to call this function for partial derivatives that are not wrt states
96 
97  Output:
98  neiBFaceGlobalCompact_: neibough face global index for a given local boundary face
99 
100  stateBoundaryCon_: matrix to store boundary connectivity levels for state Jacobians
101 
102  stateBoundaryConID_: matrix to store boundary connectivity ID for state Jacobians
103  */
104 
105  // Calculate the boundary connectivity
106  if (daOption_.getOption<label>("debug"))
107  {
108  Info << "Generating Connectivity for Boundaries:" << endl;
109  }
110 
111  this->calcNeiBFaceGlobalCompact(neiBFaceGlobalCompact_);
112 
113  this->setupStateBoundaryCon(&stateBoundaryCon_);
114 
115  this->setupStateBoundaryConID(&stateBoundaryConID_);
116 
117  wordList writeJacobians;
118  daOption_.getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
119  if (writeJacobians.found("stateBoundaryCon"))
120  {
121  DAUtility::writeMatrixBinary(stateBoundaryCon_, "stateBoundaryCon");
122  DAUtility::writeMatrixBinary(stateBoundaryConID_, "stateBoundaryConID");
123  }
124 }
125 
126 void DAJacCon::setupJacobianConnections(
127  Mat conMat,
128  Mat connections,
129  const PetscInt idxI)
130 {
131  /*
132  Description:
133  Assign connectivity to Jacobian conMat, e.g., dRdWCon, based on the connections input Mat
134 
135  Input:
136  idxI: Row index to added, ad the column index to added is based on connections
137 
138  connections: the one row matrix with nonzero values to add to conMat
139 
140  Output:
141  conMat: the connectivity mat to add
142  */
143 
144  PetscInt nCols;
145  const PetscInt* cols;
146  const PetscScalar* vals;
147 
148  MatAssemblyBegin(connections, MAT_FINAL_ASSEMBLY);
149  MatAssemblyEnd(connections, MAT_FINAL_ASSEMBLY);
150 
151  MatGetRow(connections, 0, &nCols, &cols, &vals);
152  MatSetValues(conMat, 1, &idxI, nCols, cols, vals, INSERT_VALUES);
153 
154  // restore the row of the matrix
155  MatRestoreRow(connections, 0, &nCols, &cols, &vals);
156  MatDestroy(&connections);
157 
158  return;
159 }
160 
161 void DAJacCon::createConnectionMat(Mat* connectedStates)
162 {
163  /*
164  Description:
165  Initialize a serial connectivity matrix connectedStates,
166  basically, it is one row of connectivity in DAJacCon::jacCon_
167 
168  Input/Output:
169  connectedStates: a 1 row matrix that will be used to store the connectivity
170  */
171 
172  // create a local matrix to store this row's connectivity
173  MatCreateSeqAIJ(
174  PETSC_COMM_SELF,
175  1,
176  daIndex_.nGlobalAdjointStates,
177  2000,
178  NULL,
179  connectedStates);
180  //MatSetOption(*connectedStates, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
181  MatSetUp(*connectedStates);
182 
183  MatZeroEntries(*connectedStates);
184 
185  return;
186 }
187 
188 void DAJacCon::addStateConnections(
189  Mat connections,
190  const label cellI,
191  const label connectedLevelLocal,
192  const wordList connectedStatesLocal,
193  const List<List<word>> connectedStatesInterProc,
194  const label addFace)
195 {
196  /*
197  Description:
198  A high level interface to add the connectivity for the row matrix connections
199  Note: the connections mat is basically one row of connectivity in DAJacCon::jacCon_
200 
201  Input:
202  cellI: cell index based on which we want to add the connectivity. We can add any level of
203  connected states to this cellI
204 
205  connectedLevelLocal: level of local connectivity, this is useually obtained from
206  DAJacCon::adjStateResidualConInfo_
207 
208  connectedStatesLocal: list of connected states to add for the level: connectedLevelLocal
209 
210  connectedStatesInterProc: list of states to add for a given level of boundary connectivity
211 
212  addFace: add cell faces for the current level?
213 
214  Output:
215  connections: one row of connectivity in DAJacCon::jacCon_
216 
217  Example:
218  If the connectivity list reads:
219 
220  adjStateResidualConInfo_
221  {
222  "URes"
223  {
224  {"U", "p", "phi"}, // level 0 connectivity
225  {"U", "p", "phi"}, // level 1 connectivity
226  {"U"}, // level 2 connectivity
227  }
228  }
229 
230  and the cell topology with a inter-proc boundary cen be either of the following:
231  CASE 1:
232  ---------
233  | cellQ |
234  -----------------------
235  | cellP | cellJ | cellO | <------ proc1
236  ------------------------------------------------ <----- inter-processor boundary
237  | cellT | cellK | cellI | cellL | cellU | <------ proc0
238  -----------------------------------------
239  | cellN | cellM | cellR |
240  ------------------------
241  | cellS |
242  ---------
243 
244  CASE 2:
245  ---------
246  | cellQ | <------ proc1
247  -------------------------------------------------- <----- inter-processor boundary
248  | cellP | cellJ | cellO | <------ proc0
249  -----------------------------------------
250  | cellT | cellK | cellI | cellL | cellU |
251  -----------------------------------------
252  | cellN | cellM | cellR |
253  ------------------------
254  | cellS |
255  ---------
256 
257  Then, to add the connectivity correctly, we need to add all levels of connected
258  states for cellI.
259  Level 0 connectivity is straightforward becasue we don't need
260  to provide connectedStatesInterProc
261 
262  To add level 1 connectivity, we need to:
263  set connectedLevelLocal = 1
264  set connectedStatesLocal = {U, p}
265  set connectedStatesInterProc = {{U,p}, {U}}
266  set addFace = 1
267  NOTE: we need set level 1 and level 2 con in connectedStatesInterProc because the
268  north face of cellI is a inter-proc boundary and there are two levels of connected
269  state on the other side of the inter-proc boundary for CASE 1. This is the only chance we
270  can add all two levels of connected state across the boundary for CASE 1. For CASE 2, we won't
271  add any level 1 inter-proc states because non of the faces for cellI are inter-proc
272  faces so calling DAJacCon::addBoundaryFaceConnections for cellI won't add anything
273 
274  To add level 2 connectivity, we need to
275  set connectedLevelLocal = 2
276  set connectedStatesLocal = {U}
277  set connectedStatesInterProc = {{U}}
278  set addFace = 0
279  NOTE 1: we need only level 2 con (U) for connectedStatesInterProc because if we are in CASE 1,
280  the level 2 of inter-proc states have been added. For CASE 2, we only need to add cellQ
281  by calling DAJacCon::addBoundaryFaceConnections with cellJ
282  NOTE 2: If we didn't call two levels of connectedStatesInterProc in the previous call for
283  level 1 con, we can not add it for connectedLevelLocal = 2 becasue for CASE 2 there is no
284  inter-proc boundary for cellI
285 
286  NOTE: how to provide connectedLevelLocal, connectedStatesLocal, and connectedStatesInterProc
287  are done in DAJacCon::setupJacCon
288 
289  */
290 
291  // check if the input parameters are valid
292  if (connectedLevelLocal > 3 or connectedLevelLocal < 0)
293  {
294  FatalErrorIn("connectedLevelLocal not valid") << abort(FatalError);
295  }
296  if (addFace != 0 && addFace != 1)
297  {
298  FatalErrorIn("addFace not valid") << abort(FatalError);
299  }
300  if (cellI >= mesh_.nCells())
301  {
302  FatalErrorIn("cellI not valid") << abort(FatalError);
303  }
304  //if (connectedLevelLocal>=2 && addFace==1)
305  //{
306  // FatalErrorIn("addFace not supported for localLevel>=2")<< abort(FatalError);
307  //}
308 
309  labelList val1 = {1};
310  labelList vals2 = {1, 1};
311  labelList vals3 = {1, 1, 1};
312 
313  label interProcLevel = connectedStatesInterProc.size();
314 
315  if (connectedLevelLocal == 0)
316  {
317  // add connectedStatesLocal for level0
318  forAll(connectedStatesLocal, idxI)
319  {
320  word stateName = connectedStatesLocal[idxI];
321  label compMax = 1;
322  if (daIndex_.adjStateType[stateName] == "volVectorState")
323  {
324  compMax = 3;
325  }
326  for (label i = 0; i < compMax; i++)
327  {
328  label idxJ = daIndex_.getGlobalAdjointStateIndex(stateName, cellI, i);
329  this->setConnections(connections, idxJ);
330  }
331  }
332  // add faces for level0
333  if (addFace)
334  {
335  forAll(stateInfo_["surfaceScalarStates"], idxI)
336  {
337  word stateName = stateInfo_["surfaceScalarStates"][idxI];
338  this->addConMatCellFaces(connections, 0, cellI, stateName, 1.0);
339  }
340  }
341  }
342  else if (connectedLevelLocal == 1)
343  {
344  // add connectedStatesLocal for level1
345  forAll(connectedStatesLocal, idxI)
346  {
347  word stateName = connectedStatesLocal[idxI];
348  this->addConMatNeighbourCells(connections, 0, cellI, stateName, 1.0);
349  }
350 
351  // add faces for level1
352  if (addFace)
353  {
354  forAll(mesh_.cellCells()[cellI], cellJ)
355  {
356  label localCell = mesh_.cellCells()[cellI][cellJ];
357  forAll(stateInfo_["surfaceScalarStates"], idxI)
358  {
359  word stateName = stateInfo_["surfaceScalarStates"][idxI];
360  this->addConMatCellFaces(connections, 0, localCell, stateName, 1.0);
361  }
362  }
363  }
364  // add inter-proc connectivity for level1
365  if (interProcLevel == 0)
366  {
367  // pass, not adding anything
368  }
369  else if (interProcLevel == 1)
370  {
371  this->addBoundaryFaceConnections(
372  connections,
373  0,
374  cellI,
375  val1,
376  connectedStatesInterProc,
377  addFace);
378  }
379  else if (interProcLevel == 2)
380  {
381  this->addBoundaryFaceConnections(
382  connections,
383  0,
384  cellI,
385  vals2,
386  connectedStatesInterProc,
387  addFace);
388  }
389  else if (interProcLevel == 3)
390  {
391  this->addBoundaryFaceConnections(
392  connections,
393  0,
394  cellI,
395  vals3,
396  connectedStatesInterProc,
397  addFace);
398  }
399  else
400  {
401  FatalErrorIn("interProcLevel not valid") << abort(FatalError);
402  }
403  }
404  else if (connectedLevelLocal == 2)
405  {
406  forAll(mesh_.cellCells()[cellI], cellJ)
407  {
408  label localCell = mesh_.cellCells()[cellI][cellJ];
409 
410  // add connectedStatesLocal for level2
411  forAll(connectedStatesLocal, idxI)
412  {
413  word stateName = connectedStatesLocal[idxI];
414  this->addConMatNeighbourCells(connections, 0, localCell, stateName, 1.0);
415  }
416 
417  // add faces for level2
418  if (addFace)
419  {
420  forAll(mesh_.cellCells()[localCell], cellK)
421  {
422  label localCellK = mesh_.cellCells()[localCell][cellK];
423  forAll(stateInfo_["surfaceScalarStates"], idxI)
424  {
425  word stateName = stateInfo_["surfaceScalarStates"][idxI];
426  this->addConMatCellFaces(connections, 0, localCellK, stateName, 1.0);
427  }
428  }
429  }
430 
431  // add inter-proc connecitivty for level2
432  if (interProcLevel == 0)
433  {
434  // pass, not adding anything
435  }
436  else if (interProcLevel == 1)
437  {
438  this->addBoundaryFaceConnections(
439  connections,
440  0,
441  localCell,
442  val1,
443  connectedStatesInterProc,
444  addFace);
445  }
446  else if (interProcLevel == 2)
447  {
448  this->addBoundaryFaceConnections(
449  connections,
450  0,
451  localCell,
452  vals2,
453  connectedStatesInterProc,
454  addFace);
455  }
456  else if (interProcLevel == 3)
457  {
458  this->addBoundaryFaceConnections(
459  connections,
460  0,
461  localCell,
462  vals3,
463  connectedStatesInterProc,
464  addFace);
465  }
466  else
467  {
468  FatalErrorIn("interProcLevel not valid") << abort(FatalError);
469  }
470  }
471  }
472  else if (connectedLevelLocal == 3)
473  {
474 
475  forAll(mesh_.cellCells()[cellI], cellJ)
476  {
477  label localCell = mesh_.cellCells()[cellI][cellJ];
478  forAll(mesh_.cellCells()[localCell], cellK)
479  {
480  label localCell2 = mesh_.cellCells()[localCell][cellK];
481 
482  // add connectedStatesLocal for level3
483  forAll(connectedStatesLocal, idxI)
484  {
485  word stateName = connectedStatesLocal[idxI];
486  this->addConMatNeighbourCells(connections, 0, localCell2, stateName, 1.0);
487  }
488 
489  // add faces for level3
490  if (addFace)
491  {
492  forAll(mesh_.cellCells()[localCell2], cellL)
493  {
494  label localCellL = mesh_.cellCells()[localCell2][cellL];
495  forAll(stateInfo_["surfaceScalarStates"], idxI)
496  {
497  word stateName = stateInfo_["surfaceScalarStates"][idxI];
498  this->addConMatCellFaces(connections, 0, localCellL, stateName, 1.0);
499  }
500  }
501  }
502 
503  // add inter-proc connecitivty for level3
504  if (interProcLevel == 0)
505  {
506  // pass, not adding anything
507  }
508  else if (interProcLevel == 1)
509  {
510  this->addBoundaryFaceConnections(
511  connections,
512  0,
513  localCell2,
514  val1,
515  connectedStatesInterProc,
516  addFace);
517  }
518  else if (interProcLevel == 2)
519  {
520  this->addBoundaryFaceConnections(
521  connections,
522  0,
523  localCell2,
524  vals2,
525  connectedStatesInterProc,
526  addFace);
527  }
528  else if (interProcLevel == 3)
529  {
530  this->addBoundaryFaceConnections(
531  connections,
532  0,
533  localCell2,
534  vals3,
535  connectedStatesInterProc,
536  addFace);
537  }
538  else
539  {
540  FatalErrorIn("interProcLevel not valid") << abort(FatalError);
541  }
542  }
543  }
544  }
545  else
546  {
547  FatalErrorIn("connectedLevelLocal not valid") << abort(FatalError);
548  }
549 
550  return;
551 }
552 
553 void DAJacCon::setConnections(
554  Mat conMat,
555  const label idx) const
556 {
557 
558  /*
559  Description:
560  Set 1.0 for conMat, the column index is idx, the row index is
561  always 1 because conMat is a row matrix
562  */
563 
564  PetscInt idxI = 0;
565  PetscScalar v = 1;
566  MatSetValues(conMat, 1, &idxI, 1, &idx, &v, INSERT_VALUES);
567  return;
568 }
569 
570 void DAJacCon::calcNeiBFaceGlobalCompact(labelList& neiBFaceGlobalCompact)
571 {
572  /*
573  Description:
574  This function calculates DAJacCon::neiBFaceGlobalCompact[bFaceI]. Here neiBFaceGlobalCompact
575  stores the global coupled boundary face index for the face on the other side of the local
576  processor boundary. bFaceI is the "compact" face index. bFaceI=0 for the first boundary face
577  neiBFaceGlobalCompat.size() = nLocalBoundaryFaces
578  neiBFaceGlobalCompact[bFaceI] = -1 means it is not a coupled face
579  NOTE: neiBFaceGlobalCompact will be used to calculate the connectivity across processors
580  in DAJacCon::setupStateBoundaryCon
581 
582  Output:
583  neiBFaceGlobalCompact: the global coupled boundary face index for the face on the other
584  side of the local processor boundary
585 
586  Example:
587  On proc0, neiBFaceGlobalCompact[0] = 1024, then we have the following:
588 
589  localBFaceI = 0 <--proc0
590  --------------------------- coupled boundary face
591  globalBFaceI=1024 <--proc1
592  Taken and modified from the extended stencil code in fvMesh
593  Swap the global boundary face index
594  */
595 
596  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
597 
598  neiBFaceGlobalCompact.setSize(daIndex_.nLocalBoundaryFaces);
599 
600  // initialize the list with -1, i.e., non coupled face
601  forAll(neiBFaceGlobalCompact, idx)
602  {
603  neiBFaceGlobalCompact[idx] = -1;
604  }
605 
606  // loop over the patches and store the global indices
607  label counter = 0;
608  forAll(patches, patchI)
609  {
610  const polyPatch& pp = patches[patchI];
611 
612  // get the start index of this patch in the global face list
613  label faceIStart = pp.start();
614 
615  // check whether this face is coupled (cyclic or processor?)
616  if (pp.coupled())
617  {
618  // For coupled faces set the global face index so that it can be
619  // swaped across the interface.
620  forAll(pp, i)
621  {
622  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
623  neiBFaceGlobalCompact[bFaceI] = daIndex_.globalCoupledBFaceNumbering.toGlobal(counter);
624  faceIStart++;
625  counter++;
626  }
627  }
628  }
629 
630  // Swap the cell indices, the list now contains the global index for the
631  // U state for the cell on the other side of the processor boundary
632  syncTools::swapBoundaryFaceList(mesh_, neiBFaceGlobalCompact);
633 
634  return;
635 }
636 
637 label DAJacCon::getLocalCoupledBFaceIndex(const label localFaceI) const
638 {
639  /*
640  Description:
641  Calculate the index of the local inter-processor boundary face (bRow).
642 
643  Input:
644  localFaceI: The local face index. It is in a list of faces including all the
645  internal and boundary faces.
646 
647  Output:
648  bRow: A list of faces starts with the first inter-processor face.
649  See DAJacCon::globalBndNumbering_ for more details.
650  */
651 
652  label counter = 0;
653  forAll(mesh_.boundaryMesh(), patchI)
654  {
655  const polyPatch& pp = mesh_.boundaryMesh()[patchI];
656  // check whether this face is coupled (cyclic or processor?)
657  if (pp.coupled())
658  {
659  // get the start index of this patch in the global face
660  // list and the size of this patch.
661  label faceStart = pp.start();
662  label patchSize = pp.size();
663  label faceEnd = faceStart + patchSize;
664  if (localFaceI >= faceStart && localFaceI < faceEnd)
665  {
666  // this face is on this patch, find the exact index
667  label countDelta = localFaceI - pp.start(); //-faceStart;
668  PetscInt bRow = counter + countDelta;
669  return bRow;
670  }
671  else
672  {
673  //increment the counter by patchSize
674  counter += patchSize;
675  }
676  }
677  }
678 
679  // no match found
680  FatalErrorIn("getLocalBndFaceIndex") << abort(FatalError);
681  return -1;
682 }
683 
684 void DAJacCon::setupStateBoundaryCon(Mat* stateBoundaryCon)
685 {
686  /*
687  Description:
688  This function calculates DAJacCon::stateBoundaryCon_
689 
690  Output:
691  stateBoundaryCon stores the level of connected states (on the other side
692  across the boundary) for a given coupled boundary face. stateBoundaryCon is
693  a matrix with sizes of nGlobalCoupledBFaces by nGlobalAdjointStates
694  stateBoundaryCon is mainly used in the addBoundaryFaceConnection function
695 
696  Example:
697  Basically, if there are 2 levels of connected states across the inter-proc boundary
698 
699  |<-----------proc0, globalBFaceI=1024
700  ----------------------------------- <-coupled boundary face
701  globalAdjStateIdx=100 -> | lv1 | <------ proc1
702  |_____|
703  globalAdjStateIdx=200 -> | lv2 |
704  |_____|
705 
706 
707  The indices for row 1024 in the stateBoundaryCon matrix will be
708  stateBoundaryCon
709  rowI=1024
710  Cols: colI=0 ...... colI=100 ........ colI=200 ......... colI=nGlobalAdjointStates
711  Vals (level): 1 2
712  NOTE: globalBFaceI=1024 is owned by proc0
713 
714  */
715 
716  MatCreate(PETSC_COMM_WORLD, stateBoundaryCon);
717  MatSetSizes(
718  *stateBoundaryCon,
719  daIndex_.nLocalCoupledBFaces,
720  daIndex_.nLocalAdjointStates,
721  PETSC_DETERMINE,
722  PETSC_DETERMINE);
723  MatSetFromOptions(*stateBoundaryCon);
724  MatMPIAIJSetPreallocation(*stateBoundaryCon, 1000, NULL, 1000, NULL);
725  MatSeqAIJSetPreallocation(*stateBoundaryCon, 1000, NULL);
726  MatSetOption(*stateBoundaryCon, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
727  MatSetUp(*stateBoundaryCon);
728  MatZeroEntries(*stateBoundaryCon);
729 
730  Mat stateBoundaryConTmp;
731  MatCreate(PETSC_COMM_WORLD, &stateBoundaryConTmp);
732  MatSetSizes(
733  stateBoundaryConTmp,
734  daIndex_.nLocalCoupledBFaces,
735  daIndex_.nLocalAdjointStates,
736  PETSC_DETERMINE,
737  PETSC_DETERMINE);
738  MatSetFromOptions(stateBoundaryConTmp);
739  MatMPIAIJSetPreallocation(stateBoundaryConTmp, 1000, NULL, 1000, NULL);
740  MatSeqAIJSetPreallocation(stateBoundaryConTmp, 1000, NULL);
741  MatSetOption(stateBoundaryConTmp, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
742  MatSetUp(stateBoundaryConTmp);
743  MatZeroEntries(stateBoundaryConTmp);
744 
745  // loop over the patches and set the boundary connnectivity
746  // Add connectivity in reverse so that the nearer stencils take priority
747 
748  // NOTE: we need to start with level 3, then to 2, then to 1, and flush the matrix
749  // for each level before going to another level This is necessary because
750  // we need to make sure a proper INSERT_VALUE behavior in MatSetValues
751  // i.e., we found that if you use INSERT_VALUE to insert different values (e.g., 1, 2, and 3)
752  // to a same rowI and colI in MatSetValues, and call Mat_Assembly in the end. The, the actual
753  // value in rowI and colI is kind of random, it does not depend on which value is
754  // insert first, in this case, it can be 1, 2, or 3... This happens only in parallel and
755  // only happens after Petsc-3.8.4
756 
757  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
758  // level 3 con
759  forAll(patches, patchI)
760  {
761  const polyPatch& pp = patches[patchI];
762  const UList<label>& pFaceCells = pp.faceCells();
763  // get the start index of this patch in the global face list
764  label faceIStart = pp.start();
765 
766  // check whether this face is coupled (cyclic or processor?)
767  if (pp.coupled())
768  {
769  forAll(pp, faceI)
770  {
771  // get the necessary matrix row
772  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
773  faceIStart++;
774  label gRow = neiBFaceGlobalCompact_[bFaceI];
775 
776  // Now get the cell that borders this coupled bFace
777  label idxN = pFaceCells[faceI];
778 
779  // This cell is already a neighbour cell, so we need this plus two
780  // more levels
781  // Start with next to nearest neighbours
782  forAll(mesh_.cellCells()[idxN], cellI)
783  {
784  label localCell = mesh_.cellCells()[idxN][cellI];
785  forAll(daIndex_.adjStateNames, idxI)
786  {
787  word stateName = daIndex_.adjStateNames[idxI];
788  if (daIndex_.adjStateType[stateName] != "surfaceScalarState")
789  {
790  // Now add level 3 connectivity, add all vars except for
791  // surfaceScalarStates
792  this->addConMatNeighbourCells(
793  *stateBoundaryCon,
794  gRow,
795  localCell,
796  stateName,
797  3.0);
798  this->addConMatNeighbourCells(
799  stateBoundaryConTmp,
800  gRow,
801  localCell,
802  stateName,
803  3.0);
804  }
805  }
806  }
807  }
808  }
809  }
810  // NOTE: need to flush the value before assigning the next level
811  MatAssemblyBegin(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
812  MatAssemblyEnd(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
813  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
814  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
815 
816  // level 2 con
817  forAll(patches, patchI)
818  {
819  const polyPatch& pp = patches[patchI];
820  const UList<label>& pFaceCells = pp.faceCells();
821  // get the start index of this patch in the global face list
822  label faceIStart = pp.start();
823 
824  // check whether this face is coupled (cyclic or processor?)
825  if (pp.coupled())
826  {
827  forAll(pp, faceI)
828  {
829  // get the necessary matrix row
830  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
831  faceIStart++;
832  label gRow = neiBFaceGlobalCompact_[bFaceI];
833 
834  // Now get the cell that borders this coupled bFace
835  label idxN = pFaceCells[faceI];
836 
837  // now add the nearest neighbour cells, add all vars for level 2 except
838  // for surfaceScalarStates
839  forAll(daIndex_.adjStateNames, idxI)
840  {
841  word stateName = daIndex_.adjStateNames[idxI];
842  if (daIndex_.adjStateType[stateName] != "surfaceScalarState")
843  {
844  this->addConMatNeighbourCells(
845  *stateBoundaryCon,
846  gRow,
847  idxN,
848  stateName,
849  2.0);
850  this->addConMatNeighbourCells(
851  stateBoundaryConTmp,
852  gRow,
853  idxN,
854  stateName,
855  2.0);
856  }
857  }
858  }
859  }
860  }
861  // NOTE: need to flush the value before assigning the next level
862  MatAssemblyBegin(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
863  MatAssemblyEnd(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
864  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
865  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
866 
867  // face con
868  forAll(patches, patchI)
869  {
870  const polyPatch& pp = patches[patchI];
871  const UList<label>& pFaceCells = pp.faceCells();
872  // get the start index of this patch in the global face list
873  label faceIStart = pp.start();
874 
875  // check whether this face is coupled (cyclic or processor?)
876  if (pp.coupled())
877  {
878  forAll(pp, faceI)
879  {
880  // get the necessary matrix row
881  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
882  faceIStart++;
883  label gRow = neiBFaceGlobalCompact_[bFaceI];
884 
885  // Now get the cell that borders this coupled bFace
886  label idxN = pFaceCells[faceI];
887  // and add the surfaceScalarStates for idxN
888  forAll(stateInfo_["surfaceScalarStates"], idxI)
889  {
890  word stateName = stateInfo_["surfaceScalarStates"][idxI];
891  this->addConMatCellFaces(
892  *stateBoundaryCon,
893  gRow,
894  idxN,
895  stateName,
896  10.0); // for faces, its connectivity level is 10
897  this->addConMatCellFaces(
898  stateBoundaryConTmp,
899  gRow,
900  idxN,
901  stateName,
902  10.0);
903  }
904  }
905  }
906  }
907  // NOTE: need to flush the value before assigning the next level
908  MatAssemblyBegin(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
909  MatAssemblyEnd(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
910  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
911  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
912 
913  // level 1 con
914  forAll(patches, patchI)
915  {
916  const polyPatch& pp = patches[patchI];
917  const UList<label>& pFaceCells = pp.faceCells();
918  // get the start index of this patch in the global face list
919  label faceIStart = pp.start();
920 
921  // check whether this face is coupled (cyclic or processor?)
922  if (pp.coupled())
923  {
924  forAll(pp, faceI)
925  {
926  // get the necessary matrix row
927  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
928  faceIStart++;
929  label gRow = neiBFaceGlobalCompact_[bFaceI];
930 
931  // Now get the cell that borders this coupled bFace
932  label idxN = pFaceCells[faceI];
933  // Add all the cell states for idxN
934  forAll(daIndex_.adjStateNames, idxI)
935  {
936  word stateName = daIndex_.adjStateNames[idxI];
937  if (daIndex_.adjStateType[stateName] != "surfaceScalarState")
938  {
939  this->addConMatCell(
940  *stateBoundaryCon,
941  gRow,
942  idxN,
943  stateName,
944  1.0);
945  this->addConMatCell(
946  stateBoundaryConTmp,
947  gRow,
948  idxN,
949  stateName,
950  1.0);
951  }
952  }
953  }
954  }
955  }
956  // Now we can do the final assembly
957  MatAssemblyBegin(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
958  MatAssemblyEnd(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
959 
960  // Now repeat loop adding boundary connections from other procs using matrix
961  // created in the first loop.
962  // Add connectivity in reverse so that the nearer stencils take priority
963 
964  // level 3 con
965  forAll(patches, patchI)
966  {
967  const polyPatch& pp = patches[patchI];
968  const UList<label>& pFaceCells = pp.faceCells();
969  // get the start index of this patch in the global face list
970  label faceIStart = pp.start();
971 
972  // check whether this face is coupled (cyclic or processor?)
973  if (pp.coupled())
974  {
975  forAll(pp, faceI)
976  {
977  // get the necessary matrix row
978  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
979  faceIStart++;
980  label gRow = neiBFaceGlobalCompact_[bFaceI];
981 
982  // Now get the cell that borders this coupled bFace
983  label idxN = pFaceCells[faceI];
984 
985  // This cell is already a neighbour cell, so we need this plus two
986  // more levels
987  // Start with nearest neighbours
988  forAll(mesh_.cellCells()[idxN], cellI)
989  {
990  label localCell = mesh_.cellCells()[idxN][cellI];
991  labelList val1 = {3};
992  // pass a zero list to add all states
993  List<List<word>> connectedStates(0);
994  this->addBoundaryFaceConnections(
995  stateBoundaryConTmp,
996  gRow,
997  localCell,
998  val1,
999  connectedStates,
1000  0);
1001  }
1002  }
1003  }
1004  }
1005  // NOTE: need to flush the value before assigning the next level
1006  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1007  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1008 
1009  // level 2 and 3 con
1010  forAll(patches, patchI)
1011  {
1012  const polyPatch& pp = patches[patchI];
1013  const UList<label>& pFaceCells = pp.faceCells();
1014  // get the start index of this patch in the global face list
1015  label faceIStart = pp.start();
1016 
1017  // check whether this face is coupled (cyclic or processor?)
1018  if (pp.coupled())
1019  {
1020  forAll(pp, faceI)
1021  {
1022  // get the necessary matrix row
1023  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1024  faceIStart++;
1025  label gRow = neiBFaceGlobalCompact_[bFaceI];
1026 
1027  // Now get the cell that borders this coupled bFace
1028  label idxN = pFaceCells[faceI];
1029  // now add the neighbour cells
1030  labelList vals2 = {2, 3};
1031  // pass a zero list to add all states
1032  List<List<word>> connectedStates(0);
1033  this->addBoundaryFaceConnections(
1034  stateBoundaryConTmp,
1035  gRow,
1036  idxN,
1037  vals2,
1038  connectedStates,
1039  0);
1040  }
1041  }
1042  }
1043  // NOTE: need to flush the value before assigning the next level
1044  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1045  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1046 
1047  // level 2 again, because the previous call will mess up level 2 con
1048  forAll(patches, patchI)
1049  {
1050  const polyPatch& pp = patches[patchI];
1051  const UList<label>& pFaceCells = pp.faceCells();
1052  // get the start index of this patch in the global face list
1053  label faceIStart = pp.start();
1054 
1055  // check whether this face is coupled (cyclic or processor?)
1056  if (pp.coupled())
1057  {
1058  forAll(pp, faceI)
1059  {
1060  // get the necessary matrix row
1061  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1062  faceIStart++;
1063  label gRow = neiBFaceGlobalCompact_[bFaceI];
1064 
1065  // Now get the cell that borders this coupled bFace
1066  label idxN = pFaceCells[faceI];
1067  // now add the neighbour cells
1068  labelList vals1 = {2};
1069  // pass a zero list to add all states
1070  List<List<word>> connectedStates(0);
1071  this->addBoundaryFaceConnections(
1072  stateBoundaryConTmp,
1073  gRow,
1074  idxN,
1075  vals1,
1076  connectedStates,
1077  0);
1078  }
1079  }
1080  }
1081 
1082  MatAssemblyBegin(stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1083  MatAssemblyEnd(stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1084 
1085  // the above repeat loop is not enough to cover all the stencil, we need to do more
1086  this->combineStateBndCon(stateBoundaryCon, &stateBoundaryConTmp);
1087 
1088  return;
1089 }
1090 
1091 void DAJacCon::combineStateBndCon(
1092  Mat* stateBoundaryCon,
1093  Mat* stateBoundaryConTmp)
1094 {
1095  /*
1096  Description:
1097  1. Add additional adj state connectivities if the stateBoundaryCon stencil extends through
1098  three or more decomposed domains, something like this:
1099 
1100  -------- ---------
1101  | | |
1102  Con3 | Con2 | Con1 | R
1103  | | |
1104  --------- --------
1105 
1106  Here R is the residual, Con1 to 3 are its connectivity, and dashed lines
1107  are the inter-processor boundary
1108 
1109  2. Assign stateBoundaryConTmp to stateBoundaryCon.
1110 
1111  Input/Output:
1112  stateBoundaryCon, and stateBoundaryConTmp should come from DAJacCon::stateBoundaryCon
1113  */
1114 
1115  PetscInt nCols;
1116  const PetscInt* cols;
1117  const PetscScalar* vals;
1118 
1119  PetscInt nCols1;
1120  const PetscInt* cols1;
1121  const PetscScalar* vals1;
1122 
1123  // Destroy and initialize stateBoundaryCon with zeros
1124  MatDestroy(stateBoundaryCon);
1125  MatCreate(PETSC_COMM_WORLD, stateBoundaryCon);
1126  MatSetSizes(
1127  *stateBoundaryCon,
1128  daIndex_.nLocalCoupledBFaces,
1129  daIndex_.nLocalAdjointStates,
1130  PETSC_DETERMINE,
1131  PETSC_DETERMINE);
1132  MatSetFromOptions(*stateBoundaryCon);
1133  MatMPIAIJSetPreallocation(*stateBoundaryCon, 1000, NULL, 1000, NULL);
1134  MatSeqAIJSetPreallocation(*stateBoundaryCon, 1000, NULL);
1135  MatSetOption(*stateBoundaryCon, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1136  MatSetUp(*stateBoundaryCon);
1137  MatZeroEntries(*stateBoundaryCon); // initialize with zeros
1138 
1139  // assign stateBoundaryConTmp to stateBoundaryCon
1140  PetscInt Istart, Iend;
1141  MatGetOwnershipRange(*stateBoundaryConTmp, &Istart, &Iend);
1142  for (PetscInt i = Istart; i < Iend; i++)
1143  {
1144  MatGetRow(*stateBoundaryConTmp, i, &nCols, &cols, &vals);
1145  for (PetscInt j = 0; j < nCols; j++)
1146  {
1147  MatSetValue(*stateBoundaryCon, i, cols[j], vals[j], INSERT_VALUES);
1148  }
1149  MatRestoreRow(*stateBoundaryConTmp, i, &nCols, &cols, &vals);
1150  }
1151  MatAssemblyBegin(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1152  MatAssemblyEnd(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1153 
1154  MatDestroy(stateBoundaryConTmp);
1155 
1156  // copy ConMat to ConMatTmp for an extract loop
1157  MatConvert(*stateBoundaryCon, MATSAME, MAT_INITIAL_MATRIX, stateBoundaryConTmp);
1158 
1159  // We need to do another loop adding boundary connections from other procs using ConMat
1160  // this will add missing connectivity if the stateBoundaryCon stencil extends through
1161  // three more more processors
1162  // NOTE: we need to start with level 3, then to 2, then to 1, and flush the matrix
1163  // for each level before going to another level This is necessary because
1164  // we need to make sure a proper INSERT_VALUE behavior in MatSetValues
1165  // i.e., we found that if you use INSERT_VALUE to insert different values (e.g., 1, 2, and 3)
1166  // to a same rowI and colI in MatSetValues, and call Mat_Assembly in the end. The, the actual
1167  // value in rowI and colI is kind of random, it does not depend on which value is
1168  // insert first, in this case, it can be 1, 2, or 3... This happens only in parallel and
1169  // only happens after Petsc-3.8.4
1170 
1171  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1172  // level 3 con
1173  forAll(patches, patchI)
1174  {
1175  const polyPatch& pp = patches[patchI];
1176  const UList<label>& pFaceCells = pp.faceCells();
1177  label faceIStart = pp.start();
1178  if (pp.coupled())
1179  {
1180  forAll(pp, faceI)
1181  {
1182  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1183  faceIStart++;
1184  label gRow = neiBFaceGlobalCompact_[bFaceI];
1185  label idxN = pFaceCells[faceI];
1186 
1187  forAll(mesh_.cellCells()[idxN], cellI)
1188  {
1189  label localCell = mesh_.cellCells()[idxN][cellI];
1190  labelList val1 = {3};
1191  // pass a zero list to add all states
1192  List<List<word>> connectedStates(0);
1193  this->addBoundaryFaceConnections(
1194  *stateBoundaryConTmp,
1195  gRow,
1196  localCell,
1197  val1,
1198  connectedStates,
1199  0);
1200  }
1201  }
1202  }
1203  }
1204  MatAssemblyBegin(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1205  MatAssemblyEnd(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1206 
1207  // level 2, 3 con
1208  forAll(patches, patchI)
1209  {
1210  const polyPatch& pp = patches[patchI];
1211  const UList<label>& pFaceCells = pp.faceCells();
1212  label faceIStart = pp.start();
1213  if (pp.coupled())
1214  {
1215  forAll(pp, faceI)
1216  {
1217  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1218  faceIStart++;
1219  label gRow = neiBFaceGlobalCompact_[bFaceI];
1220  label idxN = pFaceCells[faceI];
1221  // now add the neighbour cells
1222  labelList vals2 = {2, 3};
1223  // pass a zero list to add all states
1224  List<List<word>> connectedStates(0);
1225  this->addBoundaryFaceConnections(
1226  *stateBoundaryConTmp,
1227  gRow,
1228  idxN,
1229  vals2,
1230  connectedStates,
1231  0);
1232  }
1233  }
1234  }
1235  MatAssemblyBegin(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1236  MatAssemblyEnd(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1237 
1238  // level 2 again, because the previous call will mess up level 2 con
1239  forAll(patches, patchI)
1240  {
1241  const polyPatch& pp = patches[patchI];
1242  const UList<label>& pFaceCells = pp.faceCells();
1243  label faceIStart = pp.start();
1244  if (pp.coupled())
1245  {
1246  forAll(pp, faceI)
1247  {
1248  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1249  faceIStart++;
1250  label gRow = neiBFaceGlobalCompact_[bFaceI];
1251  label idxN = pFaceCells[faceI];
1252  // now add the neighbour cells
1253  labelList vals1 = {2};
1254  // pass a zero list to add all states
1255  List<List<word>> connectedStates(0);
1256  this->addBoundaryFaceConnections(
1257  *stateBoundaryConTmp,
1258  gRow,
1259  idxN,
1260  vals1,
1261  connectedStates,
1262  0);
1263  }
1264  }
1265  }
1266  MatAssemblyBegin(*stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1267  MatAssemblyEnd(*stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1268 
1269  // Now stateBoundaryConTmp will have all the missing stencil. However, it will also mess
1270  // up the existing stencil in stateBoundaryCon. So we need to do a check to make sure that
1271  // stateBoundaryConTmp only add stencil, not replacing any existing stencil in stateBoundaryCon.
1272  // If anything in stateBoundaryCon is replaced, rollback the changes.
1273  Mat tmpMat; // create a temp mat
1274  MatCreate(PETSC_COMM_WORLD, &tmpMat);
1275  MatSetSizes(
1276  tmpMat,
1277  daIndex_.nLocalCoupledBFaces,
1278  daIndex_.nLocalAdjointStates,
1279  PETSC_DETERMINE,
1280  PETSC_DETERMINE);
1281  MatSetFromOptions(tmpMat);
1282  MatMPIAIJSetPreallocation(tmpMat, 1000, NULL, 1000, NULL);
1283  MatSeqAIJSetPreallocation(tmpMat, 1000, NULL);
1284  MatSetOption(tmpMat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1285  MatSetUp(tmpMat);
1286  MatZeroEntries(tmpMat); // initialize with zeros
1287  for (PetscInt i = Istart; i < Iend; i++)
1288  {
1289  MatGetRow(*stateBoundaryCon, i, &nCols, &cols, &vals);
1290  MatGetRow(*stateBoundaryConTmp, i, &nCols1, &cols1, &vals1);
1291  for (PetscInt j = 0; j < nCols1; j++)
1292  {
1293  // for each col in stateBoundaryConTmp, we need to check if there are any existing
1294  // values for the same col in stateBoundaryCon. If yes, assign the val from
1295  // stateBoundaryCon instead of stateBoundaryConTmp
1296  PetscScalar newVal = vals1[j];
1297  PetscInt newCol = cols1[j];
1298  for (PetscInt k = 0; k < nCols; k++)
1299  {
1300  if (int(cols[k]) == int(cols1[j]))
1301  {
1302  newVal = vals[k];
1303  newCol = cols[k];
1304  break;
1305  }
1306  }
1307  MatSetValue(tmpMat, i, newCol, newVal, INSERT_VALUES);
1308  }
1309  MatRestoreRow(*stateBoundaryCon, i, &nCols, &cols, &vals);
1310  MatRestoreRow(*stateBoundaryConTmp, i, &nCols1, &cols1, &vals1);
1311  }
1312  MatAssemblyBegin(tmpMat, MAT_FINAL_ASSEMBLY);
1313  MatAssemblyEnd(tmpMat, MAT_FINAL_ASSEMBLY);
1314 
1315  // copy ConMat to ConMatTmp
1316  MatDestroy(stateBoundaryCon);
1317  MatConvert(tmpMat, MATSAME, MAT_INITIAL_MATRIX, stateBoundaryCon);
1318  MatAssemblyBegin(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1319  MatAssemblyEnd(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1320 
1321  MatDestroy(stateBoundaryConTmp);
1322  MatDestroy(&tmpMat);
1323 
1324  return;
1325 }
1326 
1327 void DAJacCon::setupStateBoundaryConID(Mat* stateBoundaryConID)
1328 {
1329  /*
1330  Description:
1331  This function computes DAJacCon::stateBoundaryConID_.
1332 
1333  Output:
1334  stateBoundaryConID: it has the exactly same structure as DAJacCon::stateBoundaryCon_
1335  except that stateBoundaryConID stores the connected stateID instead of connected
1336  levels. stateBoundaryConID will be used in DAJacCon::addBoundaryFaceConnections
1337  */
1338 
1339  PetscInt nCols, colI;
1340  const PetscInt* cols;
1341  const PetscScalar* vals;
1342  PetscInt Istart, Iend;
1343 
1344  PetscScalar valIn;
1345 
1346  // assemble adjStateID4GlobalAdjIdx
1347  // adjStateID4GlobalAdjIdx stores the adjStateID for given a global adj index
1348  labelList adjStateID4GlobalAdjIdx;
1349  adjStateID4GlobalAdjIdx.setSize(daIndex_.nGlobalAdjointStates);
1350  daIndex_.calcAdjStateID4GlobalAdjIdx(adjStateID4GlobalAdjIdx);
1351 
1352  // initialize
1353  MatCreate(PETSC_COMM_WORLD, stateBoundaryConID);
1354  MatSetSizes(
1355  *stateBoundaryConID,
1356  daIndex_.nLocalCoupledBFaces,
1357  daIndex_.nLocalAdjointStates,
1358  PETSC_DETERMINE,
1359  PETSC_DETERMINE);
1360  MatSetFromOptions(*stateBoundaryConID);
1361  MatMPIAIJSetPreallocation(*stateBoundaryConID, 1000, NULL, 1000, NULL);
1362  MatSeqAIJSetPreallocation(*stateBoundaryConID, 1000, NULL);
1363  MatSetOption(*stateBoundaryConID, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1364  MatSetUp(*stateBoundaryConID);
1365  MatZeroEntries(*stateBoundaryConID);
1366 
1367  MatGetOwnershipRange(stateBoundaryCon_, &Istart, &Iend);
1368 
1369  // set stateBoundaryConID_ based on stateBoundaryCon_ and adjStateID4GlobalAdjIdx
1370  for (PetscInt i = Istart; i < Iend; i++)
1371  {
1372  MatGetRow(stateBoundaryCon_, i, &nCols, &cols, &vals);
1373  for (PetscInt j = 0; j < nCols; j++)
1374  {
1375  if (!DAUtility::isValueCloseToRef(vals[j], 0.0))
1376  {
1377  colI = cols[j];
1378  valIn = adjStateID4GlobalAdjIdx[colI];
1379  MatSetValue(*stateBoundaryConID, i, colI, valIn, INSERT_VALUES);
1380  }
1381  }
1382  MatRestoreRow(stateBoundaryCon_, i, &nCols, &cols, &vals);
1383  }
1384 
1385  adjStateID4GlobalAdjIdx.clear();
1386 
1387  MatAssemblyBegin(*stateBoundaryConID, MAT_FINAL_ASSEMBLY);
1388  MatAssemblyEnd(*stateBoundaryConID, MAT_FINAL_ASSEMBLY);
1389 
1390  return;
1391 }
1392 
1393 void DAJacCon::addConMatCell(
1394  Mat conMat,
1395  const label gRow,
1396  const label cellI,
1397  const word stateName,
1398  const PetscScalar val)
1399 {
1400  /*
1401  Description:
1402  Insert a value (val) to the connectivity Matrix (conMat)
1403  This value will be inserted at rowI=gRow
1404  The column index is dependent on the cellI and stateName
1405 
1406  Input:
1407  gRow: which row to insert the value for conMat
1408 
1409  cellI: the index of the cell to compute the column index to add
1410 
1411  stateName: the name of the state variable to compute the column index to add
1412 
1413  val: the value to add to conMat
1414 
1415  Output:
1416  conMat: the matrix to add value to
1417 
1418  Example:
1419 
1420  If we want to add value 1.0 to conMat for
1421  column={the U globalAdjointIndice of cellI} where cellI=5
1422  row = gRow = 100
1423  Then, call addConMatCell(conMat, 100, 5, "U", 1.0)
1424 
1425  -------
1426  | cellI | <----------- value 1.0 will be added to
1427  ------- column = {global index of U
1428  for cellI}
1429  */
1430 
1431  PetscInt idxJ, idxI;
1432 
1433  idxI = gRow;
1434 
1435  // find the global index of this state
1436  label compMax = 1;
1437  if (daIndex_.adjStateType[stateName] == "volVectorState")
1438  {
1439  compMax = 3;
1440  }
1441 
1442  for (label i = 0; i < compMax; i++)
1443  {
1444  idxJ = daIndex_.getGlobalAdjointStateIndex(stateName, cellI, i);
1445  // set it in the matrix
1446  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &val, INSERT_VALUES);
1447  }
1448 
1449  return;
1450 }
1451 
1452 void DAJacCon::addConMatNeighbourCells(
1453  Mat conMat,
1454  const label gRow,
1455  const label cellI,
1456  const word stateName,
1457  const PetscScalar val)
1458 {
1459 
1460  /*
1461  Description:
1462  Insert a value (val) to the connectivity Matrix (conMat)
1463  This value will be inserted at rowI=gRow
1464  The column index is dependent on the cellI, and cellI's neibough and stateName
1465 
1466  Input:
1467  gRow: which row to insert the value for conMat
1468 
1469  cellI: the index of the cell to compute the column index to add
1470 
1471  stateName: the name of the state variable to compute the column index to add
1472 
1473  val: the value to add to conMat
1474 
1475  Output:
1476  conMat: the matrix to add value to
1477 
1478  Example:
1479 
1480  If we want to add value 1.0 to conMat for
1481  columns={the U globalAdjointIndice of all the neiboughs of cellI} where cellI=5
1482  row = gRow = 100
1483  Then, call addConMatNeighbourCells(conMat, 100, 5, "U", 1.0)
1484 
1485  -------
1486  | cellL | <----------- value 1.0 will be added to
1487  ------- ------- ------- column = {global index of U
1488  | cellJ | cellI | cellK | for cellL}, similarly for all
1489  ------- ------- ------- the neiboughs of cellI
1490  | cellM |
1491  -------
1492  */
1493 
1494  label localCellJ;
1495  PetscInt idxJ, idxI;
1496 
1497  idxI = gRow;
1498  // Add the nearest neighbour cells for cell
1499  forAll(mesh_.cellCells()[cellI], cellJ)
1500  {
1501  // get the local neighbour cell
1502  localCellJ = mesh_.cellCells()[cellI][cellJ];
1503 
1504  // find the global index of this state
1505  label compMax = 1;
1506  if (daIndex_.adjStateType[stateName] == "volVectorState")
1507  {
1508  compMax = 3;
1509  }
1510  for (label i = 0; i < compMax; i++)
1511  {
1512  idxJ = daIndex_.getGlobalAdjointStateIndex(stateName, localCellJ, i);
1513  // set it in the matrix
1514  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &val, INSERT_VALUES);
1515  }
1516  }
1517 
1518  return;
1519 }
1520 
1521 void DAJacCon::addConMatCellFaces(
1522  Mat conMat,
1523  const label gRow,
1524  const label cellI,
1525  const word stateName,
1526  const PetscScalar val)
1527 {
1528 
1529  /*
1530  Description:
1531  Insert a value (val) to the connectivity Matrix (conMat)
1532  This value will be inserted at rowI=gRow
1533  The column index is dependent on the cellI's faces and stateName
1534 
1535  Input:
1536  gRow: which row to insert the value for conMat
1537 
1538  cellI: the index of the cell to compute the column index to add
1539 
1540  stateName: the name of the state variable to compute the column index to add
1541 
1542  val: the value to add to conMat
1543 
1544  Output:
1545  conMat: the matrix to add value to
1546 
1547  Example:
1548 
1549  If we want to add value 10.0 to conMat for
1550  columns={the phi globalAdjointIndice of cellI's faces} where cellI=5
1551  row = gRow = 100
1552  Then, call addConMatCell(conMat, 100, 5, "U", 1.0)
1553 
1554  -------
1555  | cellI | <----------- value 10.0 will be added to
1556  ------- column = {global adjoint index
1557  of all cellI's faces}
1558  */
1559 
1560  PetscInt idxJ, idxI;
1561  idxI = gRow;
1562 
1563  // get the faces connected to this cell, note these are in a single
1564  // list that includes all internal and boundary faces
1565  const labelList& faces = mesh_.cells()[cellI];
1566  forAll(faces, idx)
1567  {
1568  //get the appropriate index for this face
1569  label globalState = daIndex_.getGlobalAdjointStateIndex(stateName, faces[idx]);
1570  idxJ = globalState;
1571  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &val, INSERT_VALUES);
1572  }
1573 
1574  return;
1575 }
1576 
1577 void DAJacCon::addBoundaryFaceConnections(
1578  Mat conMat,
1579  const label gRow,
1580  const label cellI,
1581  const labelList v,
1582  const List<List<word>> connectedStates,
1583  const label addFaces)
1584 {
1585  /*
1586  Description:
1587  This function adds inter-proc connectivity into conMat.
1588  For all the inter-proc faces owned by cellI, get the global adj state indices
1589  from DAJacCon::stateBoundaryCon_ and then add them into conMat
1590  Col index to add: the same col index for a given row (bRowGlobal) in the stateBoundaryCon
1591  mat if the element value in the stateBoundaryCon mat is less than the input level,
1592  i.e., v.size().
1593 
1594  Input:
1595  gRow: Row index to add
1596 
1597  cellI: the cell index for getting the faces to add inter-proc connectivity, NOTE: depending on the level
1598  of requested connections, we may add inter-proc face that are not belonged to cellI
1599 
1600  v: an array denoting the desired values to add, the size of v denotes the maximal levels to add
1601 
1602  connectedStates: selectively add some states into the conMat for the current level. If its size is 0,
1603  add all the possible states (except for surfaceStates). The dimension of connectedStates is nLevel
1604  by nStates.
1605 
1606  addFaces: whether to add indices for face (phi) connectivity
1607 
1608  Example:
1609 
1610  labelList val2={1,2};
1611  PetscInt gRow=1024, idxN = 100, addFaces=1;
1612  wordListList connectedStates={{"U","p"},{"U"}};
1613  addBoundaryFaceConnections(stateBoundaryCon,gRow,idxN,vals2,connectedStates,addFaces);
1614  The above call will add 2 levels of connected states for all the inter-proc faces belonged to cellI=idxN
1615  The cols to added are: the level1 connected states (U, p) for all the inter-proc faces belonged to
1616  cellI=idxN. the level2 connected states (U only) for all the inter-proc faces belonged to cellI=idxN
1617  The valus "1" will be added to conMat for all the level1 connected states while the value "2" will be
1618  added for level2.
1619  Note: this function will also add all faces belonged to level1 of the inter-proc faces, see the
1620  following for reference
1621 
1622  -------
1623  | idxN|
1624  | | proc0, idxN=100, globalBFaceI=1024 for the south face of idxN
1625  ----------------- <----coupled boundary face
1626  add state U and p -----> | lv1 | proc1
1627  also add faces --------> | |
1628  -------
1629  | lv2 |
1630  add state U ----------> | |
1631  -------
1632 
1633 
1634  ****** NOTE: *******
1635  If the inter-proc boundary is like the following, calling this function will NOT add any
1636  inter-proc connection for idxN because there is no inter-proc boundary for cell idxN
1637 
1638  -------
1639  | idxN|
1640  | |
1641  ------- <--------- there is no inter-proc boundary for idxN, not adding anything
1642  | lv1 |
1643  | | proc0
1644  ------------------- <----coupled boundary face
1645  | lv2 | proc1
1646  | |
1647  -------
1648 
1649  */
1650 
1651  if (v.size() != connectedStates.size() && connectedStates.size() != 0)
1652  {
1653  FatalErrorIn("") << "size of v and connectedStates are not identical!"
1654  << abort(FatalError);
1655  }
1656 
1657  PetscInt idxJ, idxI, bRow, bRowGlobal;
1658  PetscInt nCols;
1659  const PetscInt* cols;
1660  const PetscScalar* vals;
1661 
1662  PetscInt nColsID;
1663  const PetscInt* colsID;
1664  const PetscScalar* valsID;
1665 
1666  // convert stateNames to stateIDs
1667  labelListList connectedStateIDs(connectedStates.size());
1668  forAll(connectedStates, idxI)
1669  {
1670  forAll(connectedStates[idxI], idxJ)
1671  {
1672  word stateName = connectedStates[idxI][idxJ];
1673  label stateID = daIndex_.adjStateID[stateName];
1674  connectedStateIDs[idxI].append(stateID);
1675  }
1676  }
1677 
1678  idxI = gRow;
1679  // get the faces connected to this cell, note these are in a single
1680  // list that includes all internal and boundary faces
1681  const labelList& faces = mesh_.cells()[cellI];
1682 
1683  //get the level
1684  label level = v.size();
1685 
1686  for (label lv = level; lv >= 1; lv--) // we need to start from the largest levels since they have higher priority
1687  {
1688  forAll(faces, faceI)
1689  {
1690  // Now deal with coupled faces
1691  label currFace = faces[faceI];
1692 
1693  if (daIndex_.isCoupledFace[currFace])
1694  {
1695  //this is a coupled face
1696 
1697  // use the boundary connectivity to figure out what is connected
1698  // to this face for this level
1699 
1700  // get bRow in boundaryCon for this face
1701  bRow = this->getLocalCoupledBFaceIndex(currFace);
1702 
1703  // get the global bRow index
1704  bRowGlobal = daIndex_.globalCoupledBFaceNumbering.toGlobal(bRow);
1705 
1706  // now extract the boundaryCon row
1707  MatGetRow(stateBoundaryCon_, bRowGlobal, &nCols, &cols, &vals);
1708  if (connectedStates.size() != 0)
1709  {
1710  // check if we need to get stateID
1711  MatGetRow(stateBoundaryConID_, bRowGlobal, &nColsID, &colsID, &valsID);
1712  }
1713 
1714  // now loop over the row and set any column that match this level
1715  // in conMat
1716  for (label i = 0; i < nCols; i++)
1717  {
1718  idxJ = cols[i];
1719  label val = round(vals[i]); // val is the connectivity level extracted from stateBoundaryCon_ at this col
1720  // selectively add some states into conMat
1721  label addState;
1722  label stateID = -9999;
1723  // check if we need to get stateID
1724  if (connectedStates.size() != 0)
1725  {
1726  stateID = round(valsID[i]);
1727  }
1728 
1729  if (connectedStates.size() == 0)
1730  {
1731  addState = 1;
1732  }
1733  else if (DAUtility::isInList<label>(stateID, connectedStateIDs[lv - 1]))
1734  {
1735  addState = 1;
1736  }
1737  else
1738  {
1739  addState = 0;
1740  }
1741  // if the level match and the state is what you want
1742  if (val == lv && addState)
1743  {
1744  // need to do v[lv-1] here since v is an array with starting index 0
1745  PetscScalar valIn = v[lv - 1];
1746  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &valIn, INSERT_VALUES);
1747  }
1748  if (val == 10 && addFaces)
1749  {
1750  // this is a necessary connection
1751  PetscScalar valIn = v[lv - 1];
1752  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &valIn, INSERT_VALUES);
1753  }
1754  }
1755 
1756  // restore the row of the matrix
1757  MatRestoreRow(stateBoundaryCon_, bRowGlobal, &nCols, &cols, &vals);
1758  if (connectedStates.size() != 0)
1759  {
1760  // check if we need to get stateID
1761  MatRestoreRow(stateBoundaryConID_, bRowGlobal, &nColsID, &colsID, &valsID);
1762  }
1763  }
1764  }
1765  }
1766 
1767  return;
1768 }
1769 
1770 label DAJacCon::coloringExists(const word postFix) const
1771 {
1772  /*
1773  Description:
1774  Check whether the coloring file exists
1775 
1776  Input:
1777  postFix: the post fix of the file name, e.g., the original
1778  name is dFdWColoring_1.bin, then the new name is
1779  dFdWColoring_drag_1.bin with postFix = _drag
1780 
1781  Output:
1782  return 1 if coloring files exist, otherwise, return 0
1783  */
1784 
1785  Info << "Checking if Coloring file exists.." << endl;
1786  label nProcs = Pstream::nProcs();
1787  word fileName = modelType_ + "Coloring" + postFix + "_" + Foam::name(nProcs);
1788  word checkFile = fileName + ".bin";
1789  std::ifstream fIn(checkFile);
1790  if (!fIn.fail())
1791  {
1792  Info << checkFile << " exists." << endl;
1793  return 1;
1794  }
1795  else
1796  {
1797  return 0;
1798  }
1799 }
1800 
1801 void DAJacCon::calcJacConColoring(const word postFix)
1802 {
1803  /*
1804  Description:
1805  Calculate the coloring for jacCon.
1806 
1807  Input:
1808  postFix: the post fix of the file name, e.g., the original
1809  name is dFdWColoring_1.bin, then the new name is
1810  dFdWColoring_drag_1.bin with postFix = _drag
1811 
1812  Output:
1813  jacConColors_: jacCon coloring and save to files.
1814  The naming convention for coloring vector is
1815  coloringVecName_nProcs.bin. This is necessary because
1816  using different CPU cores result in different jacCon
1817  and therefore different coloring
1818 
1819  nJacColors: number of jacCon colors
1820 
1821  */
1822 
1823  // first check if the file name exists, if yes, return and
1824  // don't compute the coloring
1825  Info << "Calculating " << modelType_ << " Coloring.." << endl;
1826  label nProcs = Pstream::nProcs();
1827  word fileName = modelType_ + "Coloring" + postFix + "_" + Foam::name(nProcs);
1828 
1829  VecZeroEntries(jacConColors_);
1830  if (daOption_.getOption<label>("adjUseColoring"))
1831  {
1832  // use parallelD2 coloring to compute colors
1833  daColoring_.parallelD2Coloring(jacCon_, jacConColors_, nJacConColors_);
1834  }
1835  else
1836  {
1837  // use brute force to compute colors, basically, we assign
1838  // color to its global index
1839  PetscScalar* jacConColorsArray;
1840  VecGetArray(jacConColors_, &jacConColorsArray);
1841  label Istart, Iend;
1842  VecGetOwnershipRange(jacConColors_, &Istart, &Iend);
1843  for (label i = Istart; i < Iend; i++)
1844  {
1845  label relIdx = i - Istart;
1846  jacConColorsArray[relIdx] = i * 1.0;
1847  }
1848  VecRestoreArray(jacConColors_, &jacConColorsArray);
1849 
1850  PetscReal maxVal;
1851  VecMax(jacConColors_, NULL, &maxVal);
1852  nJacConColors_ = maxVal + 1;
1853  }
1854 
1855  daColoring_.validateColoring(jacCon_, jacConColors_);
1856  Info << " nJacConColors: " << nJacConColors_ << endl;
1857  // write jacCon colors
1858  Info << "Writing Colors to " << fileName << endl;
1859  DAUtility::writeVectorBinary(jacConColors_, fileName);
1860 
1861  return;
1862 }
1863 
1864 void DAJacCon::readJacConColoring(const word postFix)
1865 {
1866  /*
1867  Description:
1868  Read the jacCon coloring from files and
1869  compute nJacConColors. The naming convention for
1870  coloring vector is coloringVecName_nProcs.bin
1871  This is necessary because using different CPU
1872  cores result in different jacCon and therefore
1873  different coloring
1874 
1875  Input:
1876  postFix: the post fix of the file name, e.g., the original
1877  name is dFdWColoring_1.bin, then the new name is
1878  dFdWColoring_drag_1.bin with postFix = _drag
1879 
1880  Output:
1881  jacConColors_: read from file
1882 
1883  nJacConColors: number of jacCon colors
1884  */
1885 
1886  label nProcs = Pstream::nProcs();
1887  word fileName = modelType_ + "Coloring" + postFix + "_" + Foam::name(nProcs);
1888  Info << "Reading Coloring " << fileName << endl;
1889 
1890  VecZeroEntries(jacConColors_);
1891  DAUtility::readVectorBinary(jacConColors_, fileName);
1892 
1893  daColoring_.validateColoring(jacCon_, jacConColors_);
1894 
1895  PetscReal maxVal;
1896  VecMax(jacConColors_, NULL, &maxVal);
1897  nJacConColors_ = maxVal + 1;
1898 }
1899 
1900 void DAJacCon::setupJacConPreallocation(const dictionary& options)
1901 {
1902  /*
1903  Description:
1904  Compute the preallocation vectors.
1905  NOTE: this need to be implemented in the child class, if not,
1906  print an error! For example, for dRdW, this needs to be implemented;
1907  however, for dFdW, no setupJacConPreallocation is needed so users
1908  shouldn't call this function at all!
1909  */
1910  FatalErrorIn("") << "setupJacConPreallocation not implemented " << endl
1911  << " in the child class for " << modelType_
1912  << abort(FatalError);
1913 }
1914 
1915 void DAJacCon::preallocatedRdW(
1916  Mat dRMat,
1917  const label transposed) const
1918 {
1919  /*
1920  Description:
1921  Preallocate the dRMat
1922  NOTE: this need to be implemented in the child class, if not,
1923  print an error! For example, for dRdW, this needs to be implemented;
1924  however, for dFdW, no preallocatedRdW is needed so users
1925  shouldn't call this function at all!
1926  */
1927  FatalErrorIn("") << "preallocatedRdW not implemented " << endl
1928  << " in the child class for " << modelType_
1929  << abort(FatalError);
1930 }
1931 
1932 void DAJacCon::setObjFuncVec(
1933  scalarList objFuncFaceValues,
1934  scalarList objFuncCellValues,
1935  Vec objFuncVec) const
1936 {
1937  /*
1938  Description:
1939  Set the objective function vector
1940  NOTE: this need to be implemented in the child class, if not,
1941  print an error! For example, for dFdW, this needs to be implemented;
1942  however, for dRdW, no setObjFuncVec is needed so users
1943  shouldn't call this function at all!
1944  */
1945  FatalErrorIn("") << "setObjFuncVec not implemented " << endl
1946  << " in the child class for " << modelType_
1947  << abort(FatalError);
1948 }
1949 
1950 void DAJacCon::calcColoredColumns(
1951  const label colorI,
1952  Vec coloredColumn) const
1953 {
1954  /*
1955  Description:
1956  Compute the colored column vector: coloredColumn. This vector will then
1957  be used to assign resVec to dRdW in DAPartDeriv::calcPartDeriv
1958 
1959  Input:
1960  colorI: the ith color index
1961 
1962  Output:
1963  coloredColumn: For a given colorI, coloredColumn vector contains the column
1964  index for non-zero elements in the DAJacCon::jacCon_ matrix. If there is
1965  no non-zero element for this color, set the value to -1
1966 
1967  Example:
1968 
1969  If the DAJacCon::jacCon_ matrix reads,
1970 
1971  color0 color1
1972  | |
1973  1 0 0 0
1974  jacCon = 0 1 1 0
1975  0 0 1 0
1976  0 0 0 1
1977  | |
1978  color0 color0
1979 
1980  and the coloring vector DAJacCon::jacConColors_ = {0, 0, 1, 0}.
1981 
1982  **************************
1983  ***** If colorI = 0 ******
1984  **************************
1985  Calling calcColoredColumns(0, coloredColumn) will return
1986 
1987  coloredColumn = {0, 1, -1, 3}
1988 
1989  In this case, we have three columns (0, 1, and 3) for color=0, the nonzero pattern is:
1990 
1991  color0
1992  |
1993  1 0 0 0
1994  jacCon = 0 1 0 0
1995  0 0 0 0
1996  0 0 0 1
1997  | |
1998  color0 color0
1999 
2000  So for the 0th, 1th, and 3th rows, the non-zero elements in the jacCon matrix are at the
2001  0th, 1th, and 3th columns, respectively, which gives coloredColumn = {0, 1, -1, 3}
2002 
2003  **************************
2004  ***** If colorI = 1 ******
2005  **************************
2006  Calling calcColoredColumns(1, coloredColumn) will return
2007 
2008  coloredColumn = {-1, 2, 2, -1}
2009 
2010  In this case, we have one column (2) for color=1, , the nonzero pattern is:
2011 
2012  color1
2013  |
2014  0 0 0 0
2015  jacCon = 0 0 1 0
2016  0 0 1 0
2017  0 0 0 0
2018 
2019  So for the 1th and 2th rows, the non-zero elements in the jacCon matrix are at the
2020  2th and 2th columns, respectively, which gives coloredColumn = {-1, 2, 2, -1}
2021 
2022  */
2023 
2024  if (daOption_.getOption<label>("adjUseColoring"))
2025  {
2026 
2027  Vec colorIdx;
2028  label Istart, Iend;
2029 
2030  /* for the current color, determine which row/column pairs match up. */
2031 
2032  // create a vector to hold the column indices associated with this color
2033  VecDuplicate(jacConColors_, &colorIdx);
2034  VecZeroEntries(colorIdx);
2035 
2036  // Start by looping over the color vector. Set each column index associated
2037  // with the current color to its own value in the color idx vector
2038  // get the values on this proc
2039  VecGetOwnershipRange(jacConColors_, &Istart, &Iend);
2040 
2041  // create the arrays to access them directly
2042  const PetscScalar* jacConColorsArray;
2043  PetscScalar* colorIdxArray;
2044  VecGetArrayRead(jacConColors_, &jacConColorsArray);
2045  VecGetArray(colorIdx, &colorIdxArray);
2046 
2047  // loop over the entries to find the ones that match this color
2048  for (label j = Istart; j < Iend; j++)
2049  {
2050  label idx = j - Istart;
2051  if (DAUtility::isValueCloseToRef(jacConColorsArray[idx], colorI * 1.0))
2052  {
2053  // use 1 based indexing here and then subtract 1 from all values in
2054  // the mat mult. This will handle the zero index case in the first row
2055  colorIdxArray[idx] = j + 1;
2056  }
2057  }
2058  VecRestoreArrayRead(jacConColors_, &jacConColorsArray);
2059  VecRestoreArray(colorIdx, &colorIdxArray);
2060 
2061  //VecAssemblyBegin(colorIdx);
2062  //VecAssemblyEnd(colorIdx);
2063 
2064  //Set coloredColumn to -1 to account for the 1 based indexing in the above loop
2065  VecSet(coloredColumn, -1);
2066  // Now do a MatVec with the conMat to get the row coloredColumn pairs.
2067  MatMultAdd(jacCon_, colorIdx, coloredColumn, coloredColumn);
2068 
2069  // destroy the temporary vector
2070  VecDestroy(&colorIdx);
2071  }
2072  else
2073  {
2074  // uncolored case, we just set all elements to colorI
2075  PetscScalar val = colorI * 1.0;
2076  VecSet(coloredColumn, val);
2077 
2078  VecAssemblyBegin(coloredColumn);
2079  VecAssemblyEnd(coloredColumn);
2080  }
2081 }
2082 
2083 void DAJacCon::checkSpecialBCs()
2084 {
2085  /*
2086  Description:
2087  Check if there is special boundary conditions that need special treatment in jacCon_
2088  */
2089 
2090  // *******************************************************************
2091  // pressureInletVelocity
2092  // *******************************************************************
2093  if (DAUtility::isInList<word>("pressureInletVelocity", daField_.specialBCs))
2094  {
2095  // we need special treatment for connectivity
2096  Info << "pressureInletVelocity detected, applying special treatment for connectivity." << endl;
2097  // initialize and compute isPIVBCState_
2098  VecCreate(PETSC_COMM_WORLD, &isPIVBCState_);
2099  VecSetSizes(isPIVBCState_, daIndex_.nLocalAdjointStates, PETSC_DECIDE);
2100  VecSetFromOptions(isPIVBCState_);
2101  VecZeroEntries(isPIVBCState_);
2102 
2103  // compute isPIVBCState_
2104  // Note we need to read the U field, instead of getting it from db
2105  // this is because coloringSolver does not read U
2106  Info << "Calculating pressureInletVelocity state vec..." << endl;
2107  volVectorField U(
2108  IOobject(
2109  "U",
2110  mesh_.time().timeName(),
2111  mesh_,
2112  IOobject::READ_IF_PRESENT,
2113  IOobject::NO_WRITE,
2114  false),
2115  mesh_,
2116  dimensionedVector("U", dimensionSet(0, 1, -1, 0, 0, 0, 0), vector::zero),
2117  zeroGradientFvPatchField<vector>::typeName);
2118 
2119  forAll(U.boundaryField(), patchI)
2120  {
2121  if (U.boundaryFieldRef()[patchI].type() == "pressureInletVelocity")
2122  {
2123  const UList<label>& pFaceCells = mesh_.boundaryMesh()[patchI].faceCells();
2124  forAll(U.boundaryFieldRef()[patchI], faceI)
2125  {
2126 
2127  label faceCellI = pFaceCells[faceI]; // lv 1 face neibouring cells
2128  this->setPIVVec(isPIVBCState_, faceCellI);
2129 
2130  forAll(mesh_.cellCells()[faceCellI], cellJ)
2131  {
2132 
2133  label faceCellJ = mesh_.cellCells()[faceCellI][cellJ]; // lv 2 face neibouring cells
2134  this->setPIVVec(isPIVBCState_, faceCellJ);
2135 
2136  forAll(mesh_.cellCells()[faceCellJ], cellK)
2137  {
2138  label faceCellK = mesh_.cellCells()[faceCellI][cellJ]; // lv 3 face neibouring cells
2139  this->setPIVVec(isPIVBCState_, faceCellK);
2140  }
2141  }
2142  }
2143  }
2144  }
2145 
2146  VecAssemblyBegin(isPIVBCState_);
2147  VecAssemblyEnd(isPIVBCState_);
2148 
2149  wordList writeJacobians;
2150  daOption_.getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
2151  if (writeJacobians.found("isPIVVec"))
2152  {
2153  DAUtility::writeVectorASCII(isPIVBCState_, "isPIVVec");
2154  }
2155  }
2156 
2157  // *******************************************************************
2158  // append more special BCs
2159  // *******************************************************************
2160 
2161  return;
2162 }
2163 
2164 void DAJacCon::setPIVVec(
2165  Vec isPIV,
2166  const label cellI)
2167 {
2168  forAll(daIndex_.adjStateNames, idxI)
2169  {
2170  // get stateName and residual names
2171  word stateName = daIndex_.adjStateNames[idxI];
2172 
2173  // check if this state is a cell state, we do surfaceScalarState separately
2174  if (daIndex_.adjStateType[stateName] == "surfaceScalarState")
2175  {
2176  forAll(mesh_.cells()[cellI], faceI)
2177  {
2178  label cellFaceI = mesh_.cells()[cellI][faceI];
2179  label rowI = daIndex_.getGlobalAdjointStateIndex(stateName, cellFaceI);
2180  VecSetValue(isPIV, rowI, 1.0, INSERT_VALUES);
2181  }
2182  }
2183  else
2184  {
2185  // if it is a vectorState, set compMax=3
2186  label compMax = 1;
2187  if (daIndex_.adjStateType[stateName] == "volVectorState")
2188  {
2189  compMax = 3;
2190  }
2191 
2192  for (label comp = 0; comp < compMax; comp++)
2193  {
2194  label rowI = daIndex_.getGlobalAdjointStateIndex(stateName, cellI, comp);
2195  VecSetValue(isPIV, rowI, 1.0, INSERT_VALUES);
2196  }
2197  }
2198  }
2199 }
2200 
2201 label DAJacCon::addPhi4PIV(
2202  const word stateName,
2203  const label idxI,
2204  const label comp)
2205 {
2206 
2207  label localI = daIndex_.getLocalAdjointStateIndex(stateName, idxI, comp);
2208  PetscScalar* isPIVArray;
2209  VecGetArray(isPIVBCState_, &isPIVArray);
2210  if (DAUtility::isValueCloseToRef(isPIVArray[localI], 1.0))
2211  {
2212  VecRestoreArray(isPIVBCState_, &isPIVArray);
2213  return 1;
2214  }
2215  else
2216  {
2217  VecRestoreArray(isPIVBCState_, &isPIVArray);
2218  return 0;
2219  }
2220 
2221  VecRestoreArray(isPIVBCState_, &isPIVArray);
2222  return 0;
2223 }
2224 
2225 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2226 
2227 } // End namespace Foam
2228 
2229 // ************************************************************************* //
Foam::DAJacCon::modelType_
const word modelType_
the name of the jacCon matrix
Definition: DAJacCon.H:46
Foam::DAJacCon::daColoring_
DAColoring daColoring_
DAColoring object.
Definition: DAJacCon.H:61
Foam::DAJacCon::setPIVVec
void setPIVVec(Vec iSPIV, const label idxI)
function used to add connectivity for pressureInletVelocity
Definition: DAJacCon.C:2164
U
U
Definition: pEqnPimpleDyM.H:60
Foam::DAUtility::readVectorBinary
static void readVectorBinary(Vec vecIn, const word prefix)
read petsc vector in binary format
Definition: DAUtility.C:226
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::DAJacCon::checkSpecialBCs
void checkSpecialBCs()
check if there is special boundary conditions that need special treatment in jacCon_
Definition: DAJacCon.C:2083
Foam::DAJacCon::calcColoredColumns
void calcColoredColumns(const label colorI, Vec coloredColumn) const
calculate the colored column vector
Definition: DAJacCon.C:1950
Foam::DAJacCon::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAJacCon.H:49
Foam::DAJacCon::jacCon_
Mat jacCon_
Jacobian connectivity mat.
Definition: DAJacCon.H:79
Foam::DAJacCon::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAJacCon.H:52
Foam::DAOption
Definition: DAOption.H:29
Foam::DAJacCon::addPhi4PIV
label addPhi4PIV(const word stateName, const label idxI, const label comp=-1)
add connectivity phi for pressureInletVelocity
Definition: DAJacCon.C:2201
Foam::DAIndex::getGlobalAdjointStateIndex
label getGlobalAdjointStateIndex(const word stateName, const label idxI, const label comp=-1) const
get global adjoint index for a given state name, cell/face indxI and its component (optional,...
Definition: DAIndex.C:661
Foam::DAJacCon::setObjFuncVec
virtual void setObjFuncVec(scalarList objFuncFaceValues, scalarList objFuncCellValues, Vec objFuncVec) const
assign values for the objective function vector based on the face and cell value lists
Definition: DAJacCon.C:1932
daOption
DAOption daOption(mesh, pyOptions_)
Foam::DAJacCon::setupStateBoundaryConID
void setupStateBoundaryConID(Mat *stateBoundaryConID)
calculate DAJacCon::stateBoundaryConID_
Definition: DAJacCon.C:1327
Foam::DAColoring::validateColoring
void validateColoring(Mat conMat, Vec colors) const
validate if there is coloring conflict
Definition: DAColoring.C:931
Foam::DAField::specialBCs
wordList specialBCs
a list that contains the names of detected special boundary conditions
Definition: DAField.H:125
Foam::DAColoring::parallelD2Coloring
void parallelD2Coloring(const Mat conMat, Vec colors, label &nColors) const
a parallel distance-2 graph coloring function
Definition: DAColoring.C:32
Foam::DAJacCon::combineStateBndCon
void combineStateBndCon(Mat *stateBoundaryCon, Mat *stateBoundaryConTmp)
combine stateBoundaryCon and stateBoundaryConTmp, this is to ensure including all connected states fo...
Definition: DAJacCon.C:1091
Foam::DAJacCon::getLocalCoupledBFaceIndex
label getLocalCoupledBFaceIndex(const label localFaceI) const
given a local face index, return the local index of the coupled boundary face
Definition: DAJacCon.C:637
Foam::DAOption::getOption
classType getOption(const word key) const
get an option from subDict and key
Definition: DAOption.H:92
Foam::DAIndex::adjStateType
HashTable< word > adjStateType
hash table of adjoint state types, e.g., volVectorState for a given state name
Definition: DAIndex.H:77
DAJacCon.H
Foam::DAJacCon::initializeStateBoundaryCon
void initializeStateBoundaryCon()
initialize state boundary connection
Definition: DAJacCon.C:88
Foam::DAIndex::isCoupledFace
labelList isCoupledFace
for a given face index, return whether this face is a coupled boundary face
Definition: DAIndex.H:186
Foam::DAJacCon::setupJacobianConnections
void setupJacobianConnections(Mat conMat, Mat connections, const label idxI)
assign values in connections to a specific row idxI in conMat
Definition: DAJacCon.C:126
Foam::DAUtility::isValueCloseToRef
static label isValueCloseToRef(const scalar val, const scalar refVal, const scalar tol=1.0e-6)
check whether a value is close to a reference value by a tolerance
Definition: DAUtility.C:651
Foam::DAJacCon::calcNeiBFaceGlobalCompact
void calcNeiBFaceGlobalCompact(labelList &neiBFaceGlobalCompact)
calculate DAJacCon::neiBFaceGlobalCompact_
Definition: DAJacCon.C:570
Foam::DAOption::getAllOptions
const dictionary & getAllOptions() const
return a reference of allOptions_ dictionary
Definition: DAOption.H:56
Foam::DAIndex::nLocalBoundaryFaces
label nLocalBoundaryFaces
local boundary face size
Definition: DAIndex.H:95
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::defineRunTimeSelectionTable
defineRunTimeSelectionTable(DAFvSource, dictionary)
solverName
word solverName
Definition: createAdjointCompressible.H:30
Foam::DAJacCon::addConMatCell
void addConMatCell(Mat conMat, const label gRow, const label cellI, const word stateName, const PetscScalar val)
add val to the gRow row in conMat, the column indice are the state (given by stateName) for a given c...
Definition: DAJacCon.C:1393
Foam::DAJacCon::addConMatNeighbourCells
void addConMatNeighbourCells(Mat conMat, const label gRow, const label cellI, const word stateName, const PetscScalar val)
add val to gRow row in conMat, column indice are the neighbouring states (given by stateName) for a g...
Definition: DAJacCon.C:1452
Foam::DAIndex::globalCoupledBFaceNumbering
globalIndex globalCoupledBFaceNumbering
Definition: DAIndex.H:171
Foam::DAUtility::writeVectorBinary
static void writeVectorBinary(const Vec vecIn, const word prefix)
write petsc vector in binary format
Definition: DAUtility.C:259
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAIndex::adjStateNames
List< word > adjStateNames
list of adjoint state names for a specific solver
Definition: DAIndex.H:74
dafoam_plot3dtransform.vals
vals
Definition: dafoam_plot3dtransform.py:39
Foam::DAJacCon::preallocatedRdW
virtual void preallocatedRdW(Mat dRMat, const label transposed) const
preallocate dRdW matrix using the preallocVec
Definition: DAJacCon.C:1915
Foam::DAModel
Definition: DAModel.H:59
Foam::DAJacCon::isPIVBCState_
Vec isPIVBCState_
a vector to show whether a state is connected to a pressureInletVelocity boundary face (3 level max)
Definition: DAJacCon.H:179
k
dimensionedScalar & k
Definition: createRefsHeatTransfer.H:6
Foam
Definition: multiFreqScalarFvPatchField.C:144
Foam::DAJacCon::stateBoundaryCon_
Mat stateBoundaryCon_
matrix to store boundary connectivity levels for state Jacobians
Definition: DAJacCon.H:70
Foam::DAJacCon::daField_
DAField daField_
DAField object.
Definition: DAJacCon.H:64
Foam::DAJacCon::daIndex_
const DAIndex & daIndex_
DAIndex object.
Definition: DAJacCon.H:58
Foam::DAIndex::adjStateID
HashTable< label > adjStateID
a unique number ID for adjoint states, it depends on the sequence of adjStateNames
Definition: DAIndex.H:125
Foam::DAIndex::nLocalInternalFaces
label nLocalInternalFaces
local internal face size
Definition: DAIndex.H:98
Foam::DAIndex::calcAdjStateID4GlobalAdjIdx
void calcAdjStateID4GlobalAdjIdx(labelList &adjStateID4GlobalAdjIdx) const
compute global list adjStateID4GlobalAdjIdx that stores the stateID for a given globalAdjIndx
Definition: DAIndex.C:912
Foam::DAJacCon::coloringExists
label coloringExists(const word postFix="") const
whether the coloring file exists
Definition: DAJacCon.C:1770
Foam::DAUtility::writeMatrixBinary
static void writeMatrixBinary(const Mat matIn, const word prefix)
write petsc matrix in binary format
Definition: DAUtility.C:355
Foam::DAJacCon::setupStateBoundaryCon
void setupStateBoundaryCon(Mat *stateBoundaryCon)
calculate DAJacCon::stateBoundaryCon_
Definition: DAJacCon.C:684
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFvSource, 0)
Foam::DAUtility::writeVectorASCII
static void writeVectorASCII(const Vec vecIn, const word prefix)
write petsc vector in ascii format
Definition: DAUtility.C:291
Foam::DAJacCon::nJacConColors_
label nJacConColors_
number of jacCon colors
Definition: DAJacCon.H:85
Foam::DAJacCon::neiBFaceGlobalCompact_
labelList neiBFaceGlobalCompact_
neibough face global index for a given local boundary face
Definition: DAJacCon.H:76
daModel
DAModel daModel(mesh, daOption)
Foam::DAJacCon::createConnectionMat
void createConnectionMat(Mat *connectedStates)
allocate connectedState matrix
Definition: DAJacCon.C:161
Foam::DAJacCon::setupJacConPreallocation
virtual void setupJacConPreallocation(const dictionary &options)
calculate the preallocation vector for initializing the JacCon mat, if necessary
Definition: DAJacCon.C:1900
Foam::DAJacCon::setConnections
void setConnections(Mat conMat, const label idx) const
add value 1 for the colume idx to conMat
Definition: DAJacCon.C:553
Foam::DAJacCon::calcJacConColoring
void calcJacConColoring(const word postFix="")
compute graph coloring for Jacobian connectivity matrix
Definition: DAJacCon.C:1801
Foam::DAJacCon::stateBoundaryConID_
Mat stateBoundaryConID_
matrix to store boundary connectivity ID for state Jacobians
Definition: DAJacCon.H:73
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
daIndex
DAIndex daIndex(mesh, daOption, daModel)
Foam::DAJacCon::jacConColors_
Vec jacConColors_
jacCon matrix colors
Definition: DAJacCon.H:82
Foam::DAJacCon::addConMatCellFaces
void addConMatCellFaces(Mat conMat, const label gRow, const label cellI, const word stateName, const PetscScalar val)
add val to gRow row in conMat, column indice are the faces (given by stateName) for a given cellI
Definition: DAJacCon.C:1521
Foam::DAIndex::nGlobalAdjointStates
label nGlobalAdjointStates
number of global adjoint states (including all cells and faces)
Definition: DAIndex.H:148
Foam::DAIndex::nLocalCoupledBFaces
label nLocalCoupledBFaces
local coupled boundary patch size
Definition: DAIndex.H:104
Foam::DAJacCon::addBoundaryFaceConnections
void addBoundaryFaceConnections(Mat conMat, const label gRow, const label cellI, const labelList v, const List< List< word >> connectedStates, const label addFaces)
add the column index of the (iner-proc) connected states and faces to conMat, given a local face inde...
Definition: DAJacCon.C:1577
Foam::DAJacCon::readJacConColoring
void readJacConColoring(const word postFix="")
read colors for JacCon
Definition: DAJacCon.C:1864
Foam::DAJacCon::addStateConnections
void addStateConnections(Mat connections, const label cellI, const label connectedLevelLocal, const wordList connectedStatesLocal, const List< List< word >> connectedStateInterProc, const label addFace)
a high-level function to add connected state column indices to the connectivity matrix
Definition: DAJacCon.C:188
Foam::DAJacCon::stateInfo_
HashTable< wordList > stateInfo_
the regState_ list from DAStateInfo object
Definition: DAJacCon.H:67
Foam::DAIndex::nLocalAdjointStates
label nLocalAdjointStates
number of local adjoint states (including all cells and faces)
Definition: DAIndex.H:145
Foam::DAJacCon::New
static autoPtr< DAJacCon > New(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAJacCon.C:47

Generated by   doxygen 1.8.17