DAJacCon.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAJacCon.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
16 
17 DAJacCon::DAJacCon(
18  const word modelType,
19  const fvMesh& mesh,
20  const DAOption& daOption,
21  const DAModel& daModel,
22  const DAIndex& daIndex)
23  : modelType_(modelType),
24  mesh_(mesh),
25  daOption_(daOption),
26  daModel_(daModel),
27  daIndex_(daIndex),
28  daColoring_(mesh, daOption, daModel, daIndex),
29  daField_(mesh, daOption, daModel, daIndex)
30 {
31  // initialize stateInfo_
32  word solverName = daOption.getOption<word>("solverName");
33  autoPtr<DAStateInfo> daStateInfo(DAStateInfo::New(solverName, mesh, daOption, daModel));
34  stateInfo_ = daStateInfo->getStateInfo();
35 
36  // check if there is special boundary conditions that need special treatment in jacCon_
37  this->checkSpecialBCs();
38 
39  this->initializeStateBoundaryCon();
40  this->initializePetscVecs();
41 }
42 
43 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
44 
45 void DAJacCon::initializeStateBoundaryCon()
46 {
47  /*
48  Description:
49  Initialize state boundary connectivity matrices and variables.
50  This will be used for Jacobian with respect to states, e.g., dRdW, and dFdW.
51 
52  NOTE: no need to call this function for partial derivatives that are not wrt states
53 
54  Output:
55  neiBFaceGlobalCompact_: neibough face global index for a given local boundary face
56 
57  stateBoundaryCon_: matrix to store boundary connectivity levels for state Jacobians
58 
59  stateBoundaryConID_: matrix to store boundary connectivity ID for state Jacobians
60  */
61 
62  // Calculate the boundary connectivity
63  if (daOption_.getOption<label>("debug"))
64  {
65  Info << "Generating Connectivity for Boundaries:" << endl;
66  }
67 
68  this->calcNeiBFaceGlobalCompact(neiBFaceGlobalCompact_);
69 
70  this->setupStateBoundaryCon(&stateBoundaryCon_);
71 
72  this->setupStateBoundaryConID(&stateBoundaryConID_);
73 
74  wordList writeJacobians;
75  daOption_.getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
76  if (writeJacobians.found("stateBoundaryCon"))
77  {
78  DAUtility::writeMatrixBinary(stateBoundaryCon_, "stateBoundaryCon");
79  DAUtility::writeMatrixBinary(stateBoundaryConID_, "stateBoundaryConID");
80  }
81 }
82 
83 void DAJacCon::initializePetscVecs()
84 {
85  /*
86  Description:
87  Initialize dRdWTPreallocOn_, dRdWTPreallocOff_, dRdWPreallocOn_,
88  dRdWPreallocOff_, and jacConColors_
89 
90  */
91 
92  // initialize the preallocation vecs
93  VecCreate(PETSC_COMM_WORLD, &dRdWTPreallocOn_);
94  VecSetSizes(dRdWTPreallocOn_, daIndex_.nLocalAdjointStates, PETSC_DECIDE);
95  VecSetFromOptions(dRdWTPreallocOn_);
96 
97  VecDuplicate(dRdWTPreallocOn_, &dRdWTPreallocOff_);
98  VecDuplicate(dRdWTPreallocOn_, &dRdWPreallocOn_);
99  VecDuplicate(dRdWTPreallocOn_, &dRdWPreallocOff_);
100 
101  // initialize coloring vectors
102 
103  // dRdW Colors
104  VecCreate(PETSC_COMM_WORLD, &jacConColors_);
105  VecSetSizes(jacConColors_, daIndex_.nLocalAdjointStates, PETSC_DECIDE);
106  VecSetFromOptions(jacConColors_);
107 
108  return;
109 }
110 
111 void DAJacCon::setupJacobianConnections(
112  Mat conMat,
113  Mat connections,
114  const PetscInt idxI)
115 {
116  /*
117  Description:
118  Assign connectivity to Jacobian conMat, e.g., dRdWCon, based on the connections input Mat
119 
120  Input:
121  idxI: Row index to added, ad the column index to added is based on connections
122 
123  connections: the one row matrix with nonzero values to add to conMat
124 
125  Output:
126  conMat: the connectivity mat to add
127  */
128 
129  PetscInt nCols;
130  const PetscInt* cols;
131  const PetscScalar* vals;
132 
133  MatAssemblyBegin(connections, MAT_FINAL_ASSEMBLY);
134  MatAssemblyEnd(connections, MAT_FINAL_ASSEMBLY);
135 
136  MatGetRow(connections, 0, &nCols, &cols, &vals);
137  MatSetValues(conMat, 1, &idxI, nCols, cols, vals, INSERT_VALUES);
138 
139  // restore the row of the matrix
140  MatRestoreRow(connections, 0, &nCols, &cols, &vals);
141  MatDestroy(&connections);
142 
143  return;
144 }
145 
146 void DAJacCon::createConnectionMat(Mat* connectedStates)
147 {
148  /*
149  Description:
150  Initialize a serial connectivity matrix connectedStates,
151  basically, it is one row of connectivity in DAJacCon::jacCon_
152 
153  Input/Output:
154  connectedStates: a 1 row matrix that will be used to store the connectivity
155  */
156 
157  // create a local matrix to store this row's connectivity
158  MatCreateSeqAIJ(
159  PETSC_COMM_SELF,
160  1,
161  daIndex_.nGlobalAdjointStates,
162  2000,
163  NULL,
164  connectedStates);
165  //MatSetOption(*connectedStates, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
166  MatSetUp(*connectedStates);
167 
168  MatZeroEntries(*connectedStates);
169 
170  return;
171 }
172 
173 void DAJacCon::preallocateJacobianMatrix(
174  Mat dRMat,
175  const Vec preallocOnProc,
176  const Vec preallocOffProc) const
177 {
178  /*
179  Description:
180  Preallocate memory for dRMat.
181 
182  Input:
183  preallocOnProc, preallocOffProc: the on and off diagonal nonzeros for
184  dRMat
185 
186  Output:
187  dRMat: matrix to preallocate memory for
188  */
189 
190  PetscScalar normOn, normOff;
191  VecNorm(preallocOnProc, NORM_2, &normOn);
192  VecNorm(preallocOffProc, NORM_2, &normOff);
193  PetscScalar normSum = normOn + normOff;
194  if (normSum < 1.0e-10)
195  {
196  FatalErrorIn("preallocateJacobianMatrix")
197  << "preallocOnProc and preallocOffProc are not allocated!"
198  << abort(FatalError);
199  }
200 
201  PetscScalar *onVec, *offVec;
202  PetscInt onSize[daIndex_.nLocalAdjointStates], offSize[daIndex_.nLocalAdjointStates];
203 
204  VecGetArray(preallocOnProc, &onVec);
205  VecGetArray(preallocOffProc, &offVec);
206  for (label i = 0; i < daIndex_.nLocalAdjointStates; i++)
207  {
208  onSize[i] = round(onVec[i]);
209  if (onSize[i] > daIndex_.nLocalAdjointStates)
210  {
211  onSize[i] = daIndex_.nLocalAdjointStates;
212  }
213  offSize[i] = round(offVec[i]) + 5; // reserve 5 more?
214  }
215 
216  VecRestoreArray(preallocOnProc, &onVec);
217  VecRestoreArray(preallocOffProc, &offVec);
218 
219  // MatMPIAIJSetPreallocation(dRMat,NULL,preallocOnProc,NULL,preallocOffProc);
220  // MatSeqAIJSetPreallocation(dRMat,NULL,preallocOnProc);
221 
222  MatMPIAIJSetPreallocation(dRMat, NULL, onSize, NULL, offSize);
223  MatSeqAIJSetPreallocation(dRMat, NULL, onSize);
224 
225  return;
226 }
227 
228 void DAJacCon::preallocatedRdW(
229  Mat dRMat,
230  const label transposed) const
231 {
232  /*
233  Description:
234  Call the DAJacCondRdW::preallocateJacobianMatrix function with the
235  correct vectors, depending on transposed
236 
237  Input:
238  transposed: whether the state Jacobian mat is transposed, i.e., it
239  is for dRdW or dRdWT (transposed)
240 
241  Output:
242  dRMat: the matrix to preallocate
243  */
244  if (transposed)
245  {
246  this->preallocateJacobianMatrix(dRMat, dRdWTPreallocOn_, dRdWTPreallocOff_);
247  }
248  else
249  {
250  this->preallocateJacobianMatrix(dRMat, dRdWPreallocOn_, dRdWPreallocOff_);
251  }
252 }
253 
254 void DAJacCon::initializeJacCon(const dictionary& options)
255 {
256  /*
257  Description:
258  Initialize the connectivity matrix and preallocate memory
259 
260  Input:
261  options: it is not used.
262  */
263 
264  MatCreate(PETSC_COMM_WORLD, &jacCon_);
265  MatSetSizes(
266  jacCon_,
267  daIndex_.nLocalAdjointStates,
268  daIndex_.nLocalAdjointStates,
269  PETSC_DETERMINE,
270  PETSC_DETERMINE);
271  MatSetFromOptions(jacCon_);
272 
273  this->preallocateJacobianMatrix(
274  jacCon_,
275  dRdWPreallocOn_,
276  dRdWPreallocOff_);
277  //MatSetOption(jacCon_, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
278  MatSetUp(jacCon_);
279 
280  if (daOption_.getOption<label>("debug"))
281  {
282  Info << "Connectivity matrix initialized." << endl;
283  }
284 }
285 
286 void DAJacCon::setupJacCon(const dictionary& options)
287 {
288  /*
289  Description:
290  Setup DAJacCon::jacCon_
291 
292  Input:
293  options.stateResConInfo: a hashtable that contains the connectivity
294  information for dRdW, usually obtained from Foam::DAStateInfo
295  */
296 
297  HashTable<List<List<word>>> stateResConInfo;
298  options.readEntry<HashTable<List<List<word>>>>("stateResConInfo", stateResConInfo);
299 
300  label isPrealloc = 0;
301  this->setupdRdWCon(stateResConInfo, isPrealloc);
302 }
303 
304 void DAJacCon::addStateConnections(
305  Mat connections,
306  const label cellI,
307  const label connectedLevelLocal,
308  const wordList connectedStatesLocal,
309  const List<List<word>> connectedStatesInterProc,
310  const label addFace)
311 {
312  /*
313  Description:
314  A high level interface to add the connectivity for the row matrix connections
315  Note: the connections mat is basically one row of connectivity in DAJacCon::jacCon_
316 
317  Input:
318  cellI: cell index based on which we want to add the connectivity. We can add any level of
319  connected states to this cellI
320 
321  connectedLevelLocal: level of local connectivity, this is useually obtained from
322  DAJacCon::adjStateResidualConInfo_
323 
324  connectedStatesLocal: list of connected states to add for the level: connectedLevelLocal
325 
326  connectedStatesInterProc: list of states to add for a given level of boundary connectivity
327 
328  addFace: add cell faces for the current level?
329 
330  Output:
331  connections: one row of connectivity in DAJacCon::jacCon_
332 
333  Example:
334  If the connectivity list reads:
335 
336  adjStateResidualConInfo_
337  {
338  "URes"
339  {
340  {"U", "p", "phi"}, // level 0 connectivity
341  {"U", "p", "phi"}, // level 1 connectivity
342  {"U"}, // level 2 connectivity
343  }
344  }
345 
346  and the cell topology with a inter-proc boundary cen be either of the following:
347  CASE 1:
348  ---------
349  | cellQ |
350  -----------------------
351  | cellP | cellJ | cellO | <------ proc1
352  ------------------------------------------------ <----- inter-processor boundary
353  | cellT | cellK | cellI | cellL | cellU | <------ proc0
354  -----------------------------------------
355  | cellN | cellM | cellR |
356  ------------------------
357  | cellS |
358  ---------
359 
360  CASE 2:
361  ---------
362  | cellQ | <------ proc1
363  -------------------------------------------------- <----- inter-processor boundary
364  | cellP | cellJ | cellO | <------ proc0
365  -----------------------------------------
366  | cellT | cellK | cellI | cellL | cellU |
367  -----------------------------------------
368  | cellN | cellM | cellR |
369  ------------------------
370  | cellS |
371  ---------
372 
373  Then, to add the connectivity correctly, we need to add all levels of connected
374  states for cellI.
375  Level 0 connectivity is straightforward becasue we don't need
376  to provide connectedStatesInterProc
377 
378  To add level 1 connectivity, we need to:
379  set connectedLevelLocal = 1
380  set connectedStatesLocal = {U, p}
381  set connectedStatesInterProc = {{U,p}, {U}}
382  set addFace = 1
383  NOTE: we need set level 1 and level 2 con in connectedStatesInterProc because the
384  north face of cellI is a inter-proc boundary and there are two levels of connected
385  state on the other side of the inter-proc boundary for CASE 1. This is the only chance we
386  can add all two levels of connected state across the boundary for CASE 1. For CASE 2, we won't
387  add any level 1 inter-proc states because non of the faces for cellI are inter-proc
388  faces so calling DAJacCon::addBoundaryFaceConnections for cellI won't add anything
389 
390  To add level 2 connectivity, we need to
391  set connectedLevelLocal = 2
392  set connectedStatesLocal = {U}
393  set connectedStatesInterProc = {{U}}
394  set addFace = 0
395  NOTE 1: we need only level 2 con (U) for connectedStatesInterProc because if we are in CASE 1,
396  the level 2 of inter-proc states have been added. For CASE 2, we only need to add cellQ
397  by calling DAJacCon::addBoundaryFaceConnections with cellJ
398  NOTE 2: If we didn't call two levels of connectedStatesInterProc in the previous call for
399  level 1 con, we can not add it for connectedLevelLocal = 2 becasue for CASE 2 there is no
400  inter-proc boundary for cellI
401 
402  NOTE: how to provide connectedLevelLocal, connectedStatesLocal, and connectedStatesInterProc
403  are done in DAJacCon::setupJacCon
404 
405  */
406 
407  // check if the input parameters are valid
408  if (connectedLevelLocal > 3 or connectedLevelLocal < 0)
409  {
410  FatalErrorIn("connectedLevelLocal not valid") << abort(FatalError);
411  }
412  if (addFace != 0 && addFace != 1)
413  {
414  FatalErrorIn("addFace not valid") << abort(FatalError);
415  }
416  if (cellI >= mesh_.nCells())
417  {
418  FatalErrorIn("cellI not valid") << abort(FatalError);
419  }
420  //if (connectedLevelLocal>=2 && addFace==1)
421  //{
422  // FatalErrorIn("addFace not supported for localLevel>=2")<< abort(FatalError);
423  //}
424 
425  labelList val1 = {1};
426  labelList vals2 = {1, 1};
427  labelList vals3 = {1, 1, 1};
428 
429  label interProcLevel = connectedStatesInterProc.size();
430 
431  if (connectedLevelLocal == 0)
432  {
433  // add connectedStatesLocal for level0
434  forAll(connectedStatesLocal, idxI)
435  {
436  word stateName = connectedStatesLocal[idxI];
437  label compMax = 1;
438  if (daIndex_.adjStateType[stateName] == "volVectorState")
439  {
440  compMax = 3;
441  }
442  for (label i = 0; i < compMax; i++)
443  {
444  label idxJ = daIndex_.getGlobalAdjointStateIndex(stateName, cellI, i);
445  this->setConnections(connections, idxJ);
446  }
447  }
448  // add faces for level0
449  if (addFace)
450  {
451  forAll(stateInfo_["surfaceScalarStates"], idxI)
452  {
453  word stateName = stateInfo_["surfaceScalarStates"][idxI];
454  this->addConMatCellFaces(connections, 0, cellI, stateName, 1.0);
455  }
456  }
457  }
458  else if (connectedLevelLocal == 1)
459  {
460  // add connectedStatesLocal for level1
461  forAll(connectedStatesLocal, idxI)
462  {
463  word stateName = connectedStatesLocal[idxI];
464  this->addConMatNeighbourCells(connections, 0, cellI, stateName, 1.0);
465  }
466 
467  // add faces for level1
468  if (addFace)
469  {
470  forAll(mesh_.cellCells()[cellI], cellJ)
471  {
472  label localCell = mesh_.cellCells()[cellI][cellJ];
473  forAll(stateInfo_["surfaceScalarStates"], idxI)
474  {
475  word stateName = stateInfo_["surfaceScalarStates"][idxI];
476  this->addConMatCellFaces(connections, 0, localCell, stateName, 1.0);
477  }
478  }
479  }
480  // add inter-proc connectivity for level1
481  if (interProcLevel == 0)
482  {
483  // pass, not adding anything
484  }
485  else if (interProcLevel == 1)
486  {
487  this->addBoundaryFaceConnections(
488  connections,
489  0,
490  cellI,
491  val1,
492  connectedStatesInterProc,
493  addFace);
494  }
495  else if (interProcLevel == 2)
496  {
497  this->addBoundaryFaceConnections(
498  connections,
499  0,
500  cellI,
501  vals2,
502  connectedStatesInterProc,
503  addFace);
504  }
505  else if (interProcLevel == 3)
506  {
507  this->addBoundaryFaceConnections(
508  connections,
509  0,
510  cellI,
511  vals3,
512  connectedStatesInterProc,
513  addFace);
514  }
515  else
516  {
517  FatalErrorIn("interProcLevel not valid") << abort(FatalError);
518  }
519  }
520  else if (connectedLevelLocal == 2)
521  {
522  forAll(mesh_.cellCells()[cellI], cellJ)
523  {
524  label localCell = mesh_.cellCells()[cellI][cellJ];
525 
526  // add connectedStatesLocal for level2
527  forAll(connectedStatesLocal, idxI)
528  {
529  word stateName = connectedStatesLocal[idxI];
530  this->addConMatNeighbourCells(connections, 0, localCell, stateName, 1.0);
531  }
532 
533  // add faces for level2
534  if (addFace)
535  {
536  forAll(mesh_.cellCells()[localCell], cellK)
537  {
538  label localCellK = mesh_.cellCells()[localCell][cellK];
539  forAll(stateInfo_["surfaceScalarStates"], idxI)
540  {
541  word stateName = stateInfo_["surfaceScalarStates"][idxI];
542  this->addConMatCellFaces(connections, 0, localCellK, stateName, 1.0);
543  }
544  }
545  }
546 
547  // add inter-proc connecitivty for level2
548  if (interProcLevel == 0)
549  {
550  // pass, not adding anything
551  }
552  else if (interProcLevel == 1)
553  {
554  this->addBoundaryFaceConnections(
555  connections,
556  0,
557  localCell,
558  val1,
559  connectedStatesInterProc,
560  addFace);
561  }
562  else if (interProcLevel == 2)
563  {
564  this->addBoundaryFaceConnections(
565  connections,
566  0,
567  localCell,
568  vals2,
569  connectedStatesInterProc,
570  addFace);
571  }
572  else if (interProcLevel == 3)
573  {
574  this->addBoundaryFaceConnections(
575  connections,
576  0,
577  localCell,
578  vals3,
579  connectedStatesInterProc,
580  addFace);
581  }
582  else
583  {
584  FatalErrorIn("interProcLevel not valid") << abort(FatalError);
585  }
586  }
587  }
588  else if (connectedLevelLocal == 3)
589  {
590 
591  forAll(mesh_.cellCells()[cellI], cellJ)
592  {
593  label localCell = mesh_.cellCells()[cellI][cellJ];
594  forAll(mesh_.cellCells()[localCell], cellK)
595  {
596  label localCell2 = mesh_.cellCells()[localCell][cellK];
597 
598  // add connectedStatesLocal for level3
599  forAll(connectedStatesLocal, idxI)
600  {
601  word stateName = connectedStatesLocal[idxI];
602  this->addConMatNeighbourCells(connections, 0, localCell2, stateName, 1.0);
603  }
604 
605  // add faces for level3
606  if (addFace)
607  {
608  forAll(mesh_.cellCells()[localCell2], cellL)
609  {
610  label localCellL = mesh_.cellCells()[localCell2][cellL];
611  forAll(stateInfo_["surfaceScalarStates"], idxI)
612  {
613  word stateName = stateInfo_["surfaceScalarStates"][idxI];
614  this->addConMatCellFaces(connections, 0, localCellL, stateName, 1.0);
615  }
616  }
617  }
618 
619  // add inter-proc connecitivty for level3
620  if (interProcLevel == 0)
621  {
622  // pass, not adding anything
623  }
624  else if (interProcLevel == 1)
625  {
626  this->addBoundaryFaceConnections(
627  connections,
628  0,
629  localCell2,
630  val1,
631  connectedStatesInterProc,
632  addFace);
633  }
634  else if (interProcLevel == 2)
635  {
636  this->addBoundaryFaceConnections(
637  connections,
638  0,
639  localCell2,
640  vals2,
641  connectedStatesInterProc,
642  addFace);
643  }
644  else if (interProcLevel == 3)
645  {
646  this->addBoundaryFaceConnections(
647  connections,
648  0,
649  localCell2,
650  vals3,
651  connectedStatesInterProc,
652  addFace);
653  }
654  else
655  {
656  FatalErrorIn("interProcLevel not valid") << abort(FatalError);
657  }
658  }
659  }
660  }
661  else
662  {
663  FatalErrorIn("connectedLevelLocal not valid") << abort(FatalError);
664  }
665 
666  return;
667 }
668 
669 void DAJacCon::setConnections(
670  Mat conMat,
671  const label idx) const
672 {
673 
674  /*
675  Description:
676  Set 1.0 for conMat, the column index is idx, the row index is
677  always 1 because conMat is a row matrix
678  */
679 
680  PetscInt idxI = 0;
681  PetscScalar v = 1;
682  MatSetValues(conMat, 1, &idxI, 1, &idx, &v, INSERT_VALUES);
683  return;
684 }
685 
686 void DAJacCon::calcNeiBFaceGlobalCompact(labelList& neiBFaceGlobalCompact)
687 {
688  /*
689  Description:
690  This function calculates DAJacCon::neiBFaceGlobalCompact[bFaceI]. Here neiBFaceGlobalCompact
691  stores the global coupled boundary face index for the face on the other side of the local
692  processor boundary. bFaceI is the "compact" face index. bFaceI=0 for the first boundary face
693  neiBFaceGlobalCompat.size() = nLocalBoundaryFaces
694  neiBFaceGlobalCompact[bFaceI] = -1 means it is not a coupled face
695  NOTE: neiBFaceGlobalCompact will be used to calculate the connectivity across processors
696  in DAJacCon::setupStateBoundaryCon
697 
698  Output:
699  neiBFaceGlobalCompact: the global coupled boundary face index for the face on the other
700  side of the local processor boundary
701 
702  Example:
703  On proc0, neiBFaceGlobalCompact[0] = 1024, then we have the following:
704 
705  localBFaceI = 0 <--proc0
706  --------------------------- coupled boundary face
707  globalBFaceI=1024 <--proc1
708  Taken and modified from the extended stencil code in fvMesh
709  Swap the global boundary face index
710  */
711 
712  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
713 
714  neiBFaceGlobalCompact.setSize(daIndex_.nLocalBoundaryFaces);
715 
716  // initialize the list with -1, i.e., non coupled face
717  forAll(neiBFaceGlobalCompact, idx)
718  {
719  neiBFaceGlobalCompact[idx] = -1;
720  }
721 
722  // loop over the patches and store the global indices
723  label counter = 0;
724  forAll(patches, patchI)
725  {
726  const polyPatch& pp = patches[patchI];
727 
728  // get the start index of this patch in the global face list
729  label faceIStart = pp.start();
730 
731  // check whether this face is coupled (cyclic or processor?)
732  if (pp.coupled())
733  {
734  // For coupled faces set the global face index so that it can be
735  // swaped across the interface.
736  forAll(pp, i)
737  {
738  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
739  neiBFaceGlobalCompact[bFaceI] = daIndex_.globalCoupledBFaceNumbering.toGlobal(counter);
740  faceIStart++;
741  counter++;
742  }
743  }
744  }
745 
746  // Swap the cell indices, the list now contains the global index for the
747  // U state for the cell on the other side of the processor boundary
748  syncTools::swapBoundaryFaceList(mesh_, neiBFaceGlobalCompact);
749 
750  return;
751 }
752 
753 label DAJacCon::getLocalCoupledBFaceIndex(const label localFaceI) const
754 {
755  /*
756  Description:
757  Calculate the index of the local inter-processor boundary face (bRow).
758 
759  Input:
760  localFaceI: The local face index. It is in a list of faces including all the
761  internal and boundary faces.
762 
763  Output:
764  bRow: A list of faces starts with the first inter-processor face.
765  See DAJacCon::globalBndNumbering_ for more details.
766  */
767 
768  label counter = 0;
769  forAll(mesh_.boundaryMesh(), patchI)
770  {
771  const polyPatch& pp = mesh_.boundaryMesh()[patchI];
772  // check whether this face is coupled (cyclic or processor?)
773  if (pp.coupled())
774  {
775  // get the start index of this patch in the global face
776  // list and the size of this patch.
777  label faceStart = pp.start();
778  label patchSize = pp.size();
779  label faceEnd = faceStart + patchSize;
780  if (localFaceI >= faceStart && localFaceI < faceEnd)
781  {
782  // this face is on this patch, find the exact index
783  label countDelta = localFaceI - pp.start(); //-faceStart;
784  PetscInt bRow = counter + countDelta;
785  return bRow;
786  }
787  else
788  {
789  //increment the counter by patchSize
790  counter += patchSize;
791  }
792  }
793  }
794 
795  // no match found
796  FatalErrorIn("getLocalBndFaceIndex") << abort(FatalError);
797  return -1;
798 }
799 
800 void DAJacCon::setupStateBoundaryCon(Mat* stateBoundaryCon)
801 {
802  /*
803  Description:
804  This function calculates DAJacCon::stateBoundaryCon_
805 
806  Output:
807  stateBoundaryCon stores the level of connected states (on the other side
808  across the boundary) for a given coupled boundary face. stateBoundaryCon is
809  a matrix with sizes of nGlobalCoupledBFaces by nGlobalAdjointStates
810  stateBoundaryCon is mainly used in the addBoundaryFaceConnection function
811 
812  Example:
813  Basically, if there are 2 levels of connected states across the inter-proc boundary
814 
815  |<-----------proc0, globalBFaceI=1024
816  ----------------------------------- <-coupled boundary face
817  globalAdjStateIdx=100 -> | lv1 | <------ proc1
818  |_____|
819  globalAdjStateIdx=200 -> | lv2 |
820  |_____|
821 
822 
823  The indices for row 1024 in the stateBoundaryCon matrix will be
824  stateBoundaryCon
825  rowI=1024
826  Cols: colI=0 ...... colI=100 ........ colI=200 ......... colI=nGlobalAdjointStates
827  Vals (level): 1 2
828  NOTE: globalBFaceI=1024 is owned by proc0
829 
830  */
831 
832  MatCreate(PETSC_COMM_WORLD, stateBoundaryCon);
833  MatSetSizes(
834  *stateBoundaryCon,
835  daIndex_.nLocalCoupledBFaces,
836  daIndex_.nLocalAdjointStates,
837  PETSC_DETERMINE,
838  PETSC_DETERMINE);
839  MatSetFromOptions(*stateBoundaryCon);
840  MatMPIAIJSetPreallocation(*stateBoundaryCon, 1000, NULL, 1000, NULL);
841  MatSeqAIJSetPreallocation(*stateBoundaryCon, 1000, NULL);
842  MatSetOption(*stateBoundaryCon, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
843  MatSetUp(*stateBoundaryCon);
844  MatZeroEntries(*stateBoundaryCon);
845 
846  Mat stateBoundaryConTmp;
847  MatCreate(PETSC_COMM_WORLD, &stateBoundaryConTmp);
848  MatSetSizes(
849  stateBoundaryConTmp,
850  daIndex_.nLocalCoupledBFaces,
851  daIndex_.nLocalAdjointStates,
852  PETSC_DETERMINE,
853  PETSC_DETERMINE);
854  MatSetFromOptions(stateBoundaryConTmp);
855  MatMPIAIJSetPreallocation(stateBoundaryConTmp, 1000, NULL, 1000, NULL);
856  MatSeqAIJSetPreallocation(stateBoundaryConTmp, 1000, NULL);
857  MatSetOption(stateBoundaryConTmp, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
858  MatSetUp(stateBoundaryConTmp);
859  MatZeroEntries(stateBoundaryConTmp);
860 
861  // loop over the patches and set the boundary connnectivity
862  // Add connectivity in reverse so that the nearer stencils take priority
863 
864  // NOTE: we need to start with level 3, then to 2, then to 1, and flush the matrix
865  // for each level before going to another level This is necessary because
866  // we need to make sure a proper INSERT_VALUE behavior in MatSetValues
867  // i.e., we found that if you use INSERT_VALUE to insert different values (e.g., 1, 2, and 3)
868  // to a same rowI and colI in MatSetValues, and call Mat_Assembly in the end. The, the actual
869  // value in rowI and colI is kind of random, it does not depend on which value is
870  // insert first, in this case, it can be 1, 2, or 3... This happens only in parallel and
871  // only happens after Petsc-3.8.4
872 
873  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
874  // level 3 con
875  forAll(patches, patchI)
876  {
877  const polyPatch& pp = patches[patchI];
878  const UList<label>& pFaceCells = pp.faceCells();
879  // get the start index of this patch in the global face list
880  label faceIStart = pp.start();
881 
882  // check whether this face is coupled (cyclic or processor?)
883  if (pp.coupled())
884  {
885  forAll(pp, faceI)
886  {
887  // get the necessary matrix row
888  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
889  faceIStart++;
890  label gRow = neiBFaceGlobalCompact_[bFaceI];
891 
892  // Now get the cell that borders this coupled bFace
893  label idxN = pFaceCells[faceI];
894 
895  // This cell is already a neighbour cell, so we need this plus two
896  // more levels
897  // Start with next to nearest neighbours
898  forAll(mesh_.cellCells()[idxN], cellI)
899  {
900  label localCell = mesh_.cellCells()[idxN][cellI];
901  forAll(daIndex_.adjStateNames, idxI)
902  {
903  word stateName = daIndex_.adjStateNames[idxI];
904  if (daIndex_.adjStateType[stateName] != "surfaceScalarState")
905  {
906  // Now add level 3 connectivity, add all vars except for
907  // surfaceScalarStates
908  this->addConMatNeighbourCells(
909  *stateBoundaryCon,
910  gRow,
911  localCell,
912  stateName,
913  3.0);
914  this->addConMatNeighbourCells(
915  stateBoundaryConTmp,
916  gRow,
917  localCell,
918  stateName,
919  3.0);
920  }
921  }
922  }
923  }
924  }
925  }
926  // NOTE: need to flush the value before assigning the next level
927  MatAssemblyBegin(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
928  MatAssemblyEnd(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
929  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
930  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
931 
932  // level 2 con
933  forAll(patches, patchI)
934  {
935  const polyPatch& pp = patches[patchI];
936  const UList<label>& pFaceCells = pp.faceCells();
937  // get the start index of this patch in the global face list
938  label faceIStart = pp.start();
939 
940  // check whether this face is coupled (cyclic or processor?)
941  if (pp.coupled())
942  {
943  forAll(pp, faceI)
944  {
945  // get the necessary matrix row
946  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
947  faceIStart++;
948  label gRow = neiBFaceGlobalCompact_[bFaceI];
949 
950  // Now get the cell that borders this coupled bFace
951  label idxN = pFaceCells[faceI];
952 
953  // now add the nearest neighbour cells, add all vars for level 2 except
954  // for surfaceScalarStates
955  forAll(daIndex_.adjStateNames, idxI)
956  {
957  word stateName = daIndex_.adjStateNames[idxI];
958  if (daIndex_.adjStateType[stateName] != "surfaceScalarState")
959  {
960  this->addConMatNeighbourCells(
961  *stateBoundaryCon,
962  gRow,
963  idxN,
964  stateName,
965  2.0);
966  this->addConMatNeighbourCells(
967  stateBoundaryConTmp,
968  gRow,
969  idxN,
970  stateName,
971  2.0);
972  }
973  }
974  }
975  }
976  }
977  // NOTE: need to flush the value before assigning the next level
978  MatAssemblyBegin(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
979  MatAssemblyEnd(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
980  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
981  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
982 
983  // face con
984  forAll(patches, patchI)
985  {
986  const polyPatch& pp = patches[patchI];
987  const UList<label>& pFaceCells = pp.faceCells();
988  // get the start index of this patch in the global face list
989  label faceIStart = pp.start();
990 
991  // check whether this face is coupled (cyclic or processor?)
992  if (pp.coupled())
993  {
994  forAll(pp, faceI)
995  {
996  // get the necessary matrix row
997  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
998  faceIStart++;
999  label gRow = neiBFaceGlobalCompact_[bFaceI];
1000 
1001  // Now get the cell that borders this coupled bFace
1002  label idxN = pFaceCells[faceI];
1003  // and add the surfaceScalarStates for idxN
1004  forAll(stateInfo_["surfaceScalarStates"], idxI)
1005  {
1006  word stateName = stateInfo_["surfaceScalarStates"][idxI];
1007  this->addConMatCellFaces(
1008  *stateBoundaryCon,
1009  gRow,
1010  idxN,
1011  stateName,
1012  10.0); // for faces, its connectivity level is 10
1013  this->addConMatCellFaces(
1014  stateBoundaryConTmp,
1015  gRow,
1016  idxN,
1017  stateName,
1018  10.0);
1019  }
1020  }
1021  }
1022  }
1023  // NOTE: need to flush the value before assigning the next level
1024  MatAssemblyBegin(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
1025  MatAssemblyEnd(*stateBoundaryCon, MAT_FLUSH_ASSEMBLY);
1026  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1027  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1028 
1029  // level 1 con
1030  forAll(patches, patchI)
1031  {
1032  const polyPatch& pp = patches[patchI];
1033  const UList<label>& pFaceCells = pp.faceCells();
1034  // get the start index of this patch in the global face list
1035  label faceIStart = pp.start();
1036 
1037  // check whether this face is coupled (cyclic or processor?)
1038  if (pp.coupled())
1039  {
1040  forAll(pp, faceI)
1041  {
1042  // get the necessary matrix row
1043  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1044  faceIStart++;
1045  label gRow = neiBFaceGlobalCompact_[bFaceI];
1046 
1047  // Now get the cell that borders this coupled bFace
1048  label idxN = pFaceCells[faceI];
1049  // Add all the cell states for idxN
1050  forAll(daIndex_.adjStateNames, idxI)
1051  {
1052  word stateName = daIndex_.adjStateNames[idxI];
1053  if (daIndex_.adjStateType[stateName] != "surfaceScalarState")
1054  {
1055  this->addConMatCell(
1056  *stateBoundaryCon,
1057  gRow,
1058  idxN,
1059  stateName,
1060  1.0);
1061  this->addConMatCell(
1062  stateBoundaryConTmp,
1063  gRow,
1064  idxN,
1065  stateName,
1066  1.0);
1067  }
1068  }
1069  }
1070  }
1071  }
1072  // Now we can do the final assembly
1073  MatAssemblyBegin(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1074  MatAssemblyEnd(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1075 
1076  // Now repeat loop adding boundary connections from other procs using matrix
1077  // created in the first loop.
1078  // Add connectivity in reverse so that the nearer stencils take priority
1079 
1080  // level 3 con
1081  forAll(patches, patchI)
1082  {
1083  const polyPatch& pp = patches[patchI];
1084  const UList<label>& pFaceCells = pp.faceCells();
1085  // get the start index of this patch in the global face list
1086  label faceIStart = pp.start();
1087 
1088  // check whether this face is coupled (cyclic or processor?)
1089  if (pp.coupled())
1090  {
1091  forAll(pp, faceI)
1092  {
1093  // get the necessary matrix row
1094  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1095  faceIStart++;
1096  label gRow = neiBFaceGlobalCompact_[bFaceI];
1097 
1098  // Now get the cell that borders this coupled bFace
1099  label idxN = pFaceCells[faceI];
1100 
1101  // This cell is already a neighbour cell, so we need this plus two
1102  // more levels
1103  // Start with nearest neighbours
1104  forAll(mesh_.cellCells()[idxN], cellI)
1105  {
1106  label localCell = mesh_.cellCells()[idxN][cellI];
1107  labelList val1 = {3};
1108  // pass a zero list to add all states
1109  List<List<word>> connectedStates(0);
1110  this->addBoundaryFaceConnections(
1111  stateBoundaryConTmp,
1112  gRow,
1113  localCell,
1114  val1,
1115  connectedStates,
1116  0);
1117  }
1118  }
1119  }
1120  }
1121  // NOTE: need to flush the value before assigning the next level
1122  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1123  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1124 
1125  // level 2 and 3 con
1126  forAll(patches, patchI)
1127  {
1128  const polyPatch& pp = patches[patchI];
1129  const UList<label>& pFaceCells = pp.faceCells();
1130  // get the start index of this patch in the global face list
1131  label faceIStart = pp.start();
1132 
1133  // check whether this face is coupled (cyclic or processor?)
1134  if (pp.coupled())
1135  {
1136  forAll(pp, faceI)
1137  {
1138  // get the necessary matrix row
1139  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1140  faceIStart++;
1141  label gRow = neiBFaceGlobalCompact_[bFaceI];
1142 
1143  // Now get the cell that borders this coupled bFace
1144  label idxN = pFaceCells[faceI];
1145  // now add the neighbour cells
1146  labelList vals2 = {2, 3};
1147  // pass a zero list to add all states
1148  List<List<word>> connectedStates(0);
1149  this->addBoundaryFaceConnections(
1150  stateBoundaryConTmp,
1151  gRow,
1152  idxN,
1153  vals2,
1154  connectedStates,
1155  0);
1156  }
1157  }
1158  }
1159  // NOTE: need to flush the value before assigning the next level
1160  MatAssemblyBegin(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1161  MatAssemblyEnd(stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1162 
1163  // level 2 again, because the previous call will mess up level 2 con
1164  forAll(patches, patchI)
1165  {
1166  const polyPatch& pp = patches[patchI];
1167  const UList<label>& pFaceCells = pp.faceCells();
1168  // get the start index of this patch in the global face list
1169  label faceIStart = pp.start();
1170 
1171  // check whether this face is coupled (cyclic or processor?)
1172  if (pp.coupled())
1173  {
1174  forAll(pp, faceI)
1175  {
1176  // get the necessary matrix row
1177  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1178  faceIStart++;
1179  label gRow = neiBFaceGlobalCompact_[bFaceI];
1180 
1181  // Now get the cell that borders this coupled bFace
1182  label idxN = pFaceCells[faceI];
1183  // now add the neighbour cells
1184  labelList vals1 = {2};
1185  // pass a zero list to add all states
1186  List<List<word>> connectedStates(0);
1187  this->addBoundaryFaceConnections(
1188  stateBoundaryConTmp,
1189  gRow,
1190  idxN,
1191  vals1,
1192  connectedStates,
1193  0);
1194  }
1195  }
1196  }
1197 
1198  MatAssemblyBegin(stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1199  MatAssemblyEnd(stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1200 
1201  // the above repeat loop is not enough to cover all the stencil, we need to do more
1202  this->combineStateBndCon(stateBoundaryCon, &stateBoundaryConTmp);
1203 
1204  return;
1205 }
1206 
1207 void DAJacCon::combineStateBndCon(
1208  Mat* stateBoundaryCon,
1209  Mat* stateBoundaryConTmp)
1210 {
1211  /*
1212  Description:
1213  1. Add additional adj state connectivities if the stateBoundaryCon stencil extends through
1214  three or more decomposed domains, something like this:
1215 
1216  -------- ---------
1217  | | |
1218  Con3 | Con2 | Con1 | R
1219  | | |
1220  --------- --------
1221 
1222  Here R is the residual, Con1 to 3 are its connectivity, and dashed lines
1223  are the inter-processor boundary
1224 
1225  2. Assign stateBoundaryConTmp to stateBoundaryCon.
1226 
1227  Input/Output:
1228  stateBoundaryCon, and stateBoundaryConTmp should come from DAJacCon::stateBoundaryCon
1229  */
1230 
1231  PetscInt nCols;
1232  const PetscInt* cols;
1233  const PetscScalar* vals;
1234 
1235  PetscInt nCols1;
1236  const PetscInt* cols1;
1237  const PetscScalar* vals1;
1238 
1239  // Destroy and initialize stateBoundaryCon with zeros
1240  MatDestroy(stateBoundaryCon);
1241  MatCreate(PETSC_COMM_WORLD, stateBoundaryCon);
1242  MatSetSizes(
1243  *stateBoundaryCon,
1244  daIndex_.nLocalCoupledBFaces,
1245  daIndex_.nLocalAdjointStates,
1246  PETSC_DETERMINE,
1247  PETSC_DETERMINE);
1248  MatSetFromOptions(*stateBoundaryCon);
1249  MatMPIAIJSetPreallocation(*stateBoundaryCon, 1000, NULL, 1000, NULL);
1250  MatSeqAIJSetPreallocation(*stateBoundaryCon, 1000, NULL);
1251  MatSetOption(*stateBoundaryCon, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1252  MatSetUp(*stateBoundaryCon);
1253  MatZeroEntries(*stateBoundaryCon); // initialize with zeros
1254 
1255  // assign stateBoundaryConTmp to stateBoundaryCon
1256  PetscInt Istart, Iend;
1257  MatGetOwnershipRange(*stateBoundaryConTmp, &Istart, &Iend);
1258  for (PetscInt i = Istart; i < Iend; i++)
1259  {
1260  MatGetRow(*stateBoundaryConTmp, i, &nCols, &cols, &vals);
1261  for (PetscInt j = 0; j < nCols; j++)
1262  {
1263  MatSetValue(*stateBoundaryCon, i, cols[j], vals[j], INSERT_VALUES);
1264  }
1265  MatRestoreRow(*stateBoundaryConTmp, i, &nCols, &cols, &vals);
1266  }
1267  MatAssemblyBegin(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1268  MatAssemblyEnd(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1269 
1270  MatDestroy(stateBoundaryConTmp);
1271 
1272  // copy ConMat to ConMatTmp for an extract loop
1273  MatConvert(*stateBoundaryCon, MATSAME, MAT_INITIAL_MATRIX, stateBoundaryConTmp);
1274 
1275  // We need to do another loop adding boundary connections from other procs using ConMat
1276  // this will add missing connectivity if the stateBoundaryCon stencil extends through
1277  // three more more processors
1278  // NOTE: we need to start with level 3, then to 2, then to 1, and flush the matrix
1279  // for each level before going to another level This is necessary because
1280  // we need to make sure a proper INSERT_VALUE behavior in MatSetValues
1281  // i.e., we found that if you use INSERT_VALUE to insert different values (e.g., 1, 2, and 3)
1282  // to a same rowI and colI in MatSetValues, and call Mat_Assembly in the end. The, the actual
1283  // value in rowI and colI is kind of random, it does not depend on which value is
1284  // insert first, in this case, it can be 1, 2, or 3... This happens only in parallel and
1285  // only happens after Petsc-3.8.4
1286 
1287  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1288  // level 3 con
1289  forAll(patches, patchI)
1290  {
1291  const polyPatch& pp = patches[patchI];
1292  const UList<label>& pFaceCells = pp.faceCells();
1293  label faceIStart = pp.start();
1294  if (pp.coupled())
1295  {
1296  forAll(pp, faceI)
1297  {
1298  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1299  faceIStart++;
1300  label gRow = neiBFaceGlobalCompact_[bFaceI];
1301  label idxN = pFaceCells[faceI];
1302 
1303  forAll(mesh_.cellCells()[idxN], cellI)
1304  {
1305  label localCell = mesh_.cellCells()[idxN][cellI];
1306  labelList val1 = {3};
1307  // pass a zero list to add all states
1308  List<List<word>> connectedStates(0);
1309  this->addBoundaryFaceConnections(
1310  *stateBoundaryConTmp,
1311  gRow,
1312  localCell,
1313  val1,
1314  connectedStates,
1315  0);
1316  }
1317  }
1318  }
1319  }
1320  MatAssemblyBegin(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1321  MatAssemblyEnd(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1322 
1323  // level 2, 3 con
1324  forAll(patches, patchI)
1325  {
1326  const polyPatch& pp = patches[patchI];
1327  const UList<label>& pFaceCells = pp.faceCells();
1328  label faceIStart = pp.start();
1329  if (pp.coupled())
1330  {
1331  forAll(pp, faceI)
1332  {
1333  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1334  faceIStart++;
1335  label gRow = neiBFaceGlobalCompact_[bFaceI];
1336  label idxN = pFaceCells[faceI];
1337  // now add the neighbour cells
1338  labelList vals2 = {2, 3};
1339  // pass a zero list to add all states
1340  List<List<word>> connectedStates(0);
1341  this->addBoundaryFaceConnections(
1342  *stateBoundaryConTmp,
1343  gRow,
1344  idxN,
1345  vals2,
1346  connectedStates,
1347  0);
1348  }
1349  }
1350  }
1351  MatAssemblyBegin(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1352  MatAssemblyEnd(*stateBoundaryConTmp, MAT_FLUSH_ASSEMBLY);
1353 
1354  // level 2 again, because the previous call will mess up level 2 con
1355  forAll(patches, patchI)
1356  {
1357  const polyPatch& pp = patches[patchI];
1358  const UList<label>& pFaceCells = pp.faceCells();
1359  label faceIStart = pp.start();
1360  if (pp.coupled())
1361  {
1362  forAll(pp, faceI)
1363  {
1364  label bFaceI = faceIStart - daIndex_.nLocalInternalFaces;
1365  faceIStart++;
1366  label gRow = neiBFaceGlobalCompact_[bFaceI];
1367  label idxN = pFaceCells[faceI];
1368  // now add the neighbour cells
1369  labelList vals1 = {2};
1370  // pass a zero list to add all states
1371  List<List<word>> connectedStates(0);
1372  this->addBoundaryFaceConnections(
1373  *stateBoundaryConTmp,
1374  gRow,
1375  idxN,
1376  vals1,
1377  connectedStates,
1378  0);
1379  }
1380  }
1381  }
1382  MatAssemblyBegin(*stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1383  MatAssemblyEnd(*stateBoundaryConTmp, MAT_FINAL_ASSEMBLY);
1384 
1385  // Now stateBoundaryConTmp will have all the missing stencil. However, it will also mess
1386  // up the existing stencil in stateBoundaryCon. So we need to do a check to make sure that
1387  // stateBoundaryConTmp only add stencil, not replacing any existing stencil in stateBoundaryCon.
1388  // If anything in stateBoundaryCon is replaced, rollback the changes.
1389  Mat tmpMat; // create a temp mat
1390  MatCreate(PETSC_COMM_WORLD, &tmpMat);
1391  MatSetSizes(
1392  tmpMat,
1393  daIndex_.nLocalCoupledBFaces,
1394  daIndex_.nLocalAdjointStates,
1395  PETSC_DETERMINE,
1396  PETSC_DETERMINE);
1397  MatSetFromOptions(tmpMat);
1398  MatMPIAIJSetPreallocation(tmpMat, 1000, NULL, 1000, NULL);
1399  MatSeqAIJSetPreallocation(tmpMat, 1000, NULL);
1400  MatSetOption(tmpMat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1401  MatSetUp(tmpMat);
1402  MatZeroEntries(tmpMat); // initialize with zeros
1403  for (PetscInt i = Istart; i < Iend; i++)
1404  {
1405  MatGetRow(*stateBoundaryCon, i, &nCols, &cols, &vals);
1406  MatGetRow(*stateBoundaryConTmp, i, &nCols1, &cols1, &vals1);
1407  for (PetscInt j = 0; j < nCols1; j++)
1408  {
1409  // for each col in stateBoundaryConTmp, we need to check if there are any existing
1410  // values for the same col in stateBoundaryCon. If yes, assign the val from
1411  // stateBoundaryCon instead of stateBoundaryConTmp
1412  PetscScalar newVal = vals1[j];
1413  PetscInt newCol = cols1[j];
1414  for (PetscInt k = 0; k < nCols; k++)
1415  {
1416  if (int(cols[k]) == int(cols1[j]))
1417  {
1418  newVal = vals[k];
1419  newCol = cols[k];
1420  break;
1421  }
1422  }
1423  MatSetValue(tmpMat, i, newCol, newVal, INSERT_VALUES);
1424  }
1425  MatRestoreRow(*stateBoundaryCon, i, &nCols, &cols, &vals);
1426  MatRestoreRow(*stateBoundaryConTmp, i, &nCols1, &cols1, &vals1);
1427  }
1428  MatAssemblyBegin(tmpMat, MAT_FINAL_ASSEMBLY);
1429  MatAssemblyEnd(tmpMat, MAT_FINAL_ASSEMBLY);
1430 
1431  // copy ConMat to ConMatTmp
1432  MatDestroy(stateBoundaryCon);
1433  MatConvert(tmpMat, MATSAME, MAT_INITIAL_MATRIX, stateBoundaryCon);
1434  MatAssemblyBegin(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1435  MatAssemblyEnd(*stateBoundaryCon, MAT_FINAL_ASSEMBLY);
1436 
1437  MatDestroy(stateBoundaryConTmp);
1438  MatDestroy(&tmpMat);
1439 
1440  return;
1441 }
1442 
1443 void DAJacCon::setupStateBoundaryConID(Mat* stateBoundaryConID)
1444 {
1445  /*
1446  Description:
1447  This function computes DAJacCon::stateBoundaryConID_.
1448 
1449  Output:
1450  stateBoundaryConID: it has the exactly same structure as DAJacCon::stateBoundaryCon_
1451  except that stateBoundaryConID stores the connected stateID instead of connected
1452  levels. stateBoundaryConID will be used in DAJacCon::addBoundaryFaceConnections
1453  */
1454 
1455  PetscInt nCols, colI;
1456  const PetscInt* cols;
1457  const PetscScalar* vals;
1458  PetscInt Istart, Iend;
1459 
1460  PetscScalar valIn;
1461 
1462  // assemble adjStateID4GlobalAdjIdx
1463  // adjStateID4GlobalAdjIdx stores the adjStateID for given a global adj index
1464  labelList adjStateID4GlobalAdjIdx;
1465  adjStateID4GlobalAdjIdx.setSize(daIndex_.nGlobalAdjointStates);
1466  daIndex_.calcAdjStateID4GlobalAdjIdx(adjStateID4GlobalAdjIdx);
1467 
1468  // initialize
1469  MatCreate(PETSC_COMM_WORLD, stateBoundaryConID);
1470  MatSetSizes(
1471  *stateBoundaryConID,
1472  daIndex_.nLocalCoupledBFaces,
1473  daIndex_.nLocalAdjointStates,
1474  PETSC_DETERMINE,
1475  PETSC_DETERMINE);
1476  MatSetFromOptions(*stateBoundaryConID);
1477  MatMPIAIJSetPreallocation(*stateBoundaryConID, 1000, NULL, 1000, NULL);
1478  MatSeqAIJSetPreallocation(*stateBoundaryConID, 1000, NULL);
1479  MatSetOption(*stateBoundaryConID, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
1480  MatSetUp(*stateBoundaryConID);
1481  MatZeroEntries(*stateBoundaryConID);
1482 
1483  MatGetOwnershipRange(stateBoundaryCon_, &Istart, &Iend);
1484 
1485  // set stateBoundaryConID_ based on stateBoundaryCon_ and adjStateID4GlobalAdjIdx
1486  for (PetscInt i = Istart; i < Iend; i++)
1487  {
1488  MatGetRow(stateBoundaryCon_, i, &nCols, &cols, &vals);
1489  for (PetscInt j = 0; j < nCols; j++)
1490  {
1491  if (!DAUtility::isValueCloseToRef(vals[j], 0.0))
1492  {
1493  colI = cols[j];
1494  valIn = adjStateID4GlobalAdjIdx[colI];
1495  MatSetValue(*stateBoundaryConID, i, colI, valIn, INSERT_VALUES);
1496  }
1497  }
1498  MatRestoreRow(stateBoundaryCon_, i, &nCols, &cols, &vals);
1499  }
1500 
1501  adjStateID4GlobalAdjIdx.clear();
1502 
1503  MatAssemblyBegin(*stateBoundaryConID, MAT_FINAL_ASSEMBLY);
1504  MatAssemblyEnd(*stateBoundaryConID, MAT_FINAL_ASSEMBLY);
1505 
1506  return;
1507 }
1508 
1509 void DAJacCon::addConMatCell(
1510  Mat conMat,
1511  const label gRow,
1512  const label cellI,
1513  const word stateName,
1514  const PetscScalar val)
1515 {
1516  /*
1517  Description:
1518  Insert a value (val) to the connectivity Matrix (conMat)
1519  This value will be inserted at rowI=gRow
1520  The column index is dependent on the cellI and stateName
1521 
1522  Input:
1523  gRow: which row to insert the value for conMat
1524 
1525  cellI: the index of the cell to compute the column index to add
1526 
1527  stateName: the name of the state variable to compute the column index to add
1528 
1529  val: the value to add to conMat
1530 
1531  Output:
1532  conMat: the matrix to add value to
1533 
1534  Example:
1535 
1536  If we want to add value 1.0 to conMat for
1537  column={the U globalAdjointIndice of cellI} where cellI=5
1538  row = gRow = 100
1539  Then, call addConMatCell(conMat, 100, 5, "U", 1.0)
1540 
1541  -------
1542  | cellI | <----------- value 1.0 will be added to
1543  ------- column = {global index of U
1544  for cellI}
1545  */
1546 
1547  PetscInt idxJ, idxI;
1548 
1549  idxI = gRow;
1550 
1551  // find the global index of this state
1552  label compMax = 1;
1553  if (daIndex_.adjStateType[stateName] == "volVectorState")
1554  {
1555  compMax = 3;
1556  }
1557 
1558  for (label i = 0; i < compMax; i++)
1559  {
1560  idxJ = daIndex_.getGlobalAdjointStateIndex(stateName, cellI, i);
1561  // set it in the matrix
1562  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &val, INSERT_VALUES);
1563  }
1564 
1565  return;
1566 }
1567 
1568 void DAJacCon::addConMatNeighbourCells(
1569  Mat conMat,
1570  const label gRow,
1571  const label cellI,
1572  const word stateName,
1573  const PetscScalar val)
1574 {
1575 
1576  /*
1577  Description:
1578  Insert a value (val) to the connectivity Matrix (conMat)
1579  This value will be inserted at rowI=gRow
1580  The column index is dependent on the cellI, and cellI's neibough and stateName
1581 
1582  Input:
1583  gRow: which row to insert the value for conMat
1584 
1585  cellI: the index of the cell to compute the column index to add
1586 
1587  stateName: the name of the state variable to compute the column index to add
1588 
1589  val: the value to add to conMat
1590 
1591  Output:
1592  conMat: the matrix to add value to
1593 
1594  Example:
1595 
1596  If we want to add value 1.0 to conMat for
1597  columns={the U globalAdjointIndice of all the neiboughs of cellI} where cellI=5
1598  row = gRow = 100
1599  Then, call addConMatNeighbourCells(conMat, 100, 5, "U", 1.0)
1600 
1601  -------
1602  | cellL | <----------- value 1.0 will be added to
1603  ------- ------- ------- column = {global index of U
1604  | cellJ | cellI | cellK | for cellL}, similarly for all
1605  ------- ------- ------- the neiboughs of cellI
1606  | cellM |
1607  -------
1608  */
1609 
1610  label localCellJ;
1611  PetscInt idxJ, idxI;
1612 
1613  idxI = gRow;
1614  // Add the nearest neighbour cells for cell
1615  forAll(mesh_.cellCells()[cellI], cellJ)
1616  {
1617  // get the local neighbour cell
1618  localCellJ = mesh_.cellCells()[cellI][cellJ];
1619 
1620  // find the global index of this state
1621  label compMax = 1;
1622  if (daIndex_.adjStateType[stateName] == "volVectorState")
1623  {
1624  compMax = 3;
1625  }
1626  for (label i = 0; i < compMax; i++)
1627  {
1628  idxJ = daIndex_.getGlobalAdjointStateIndex(stateName, localCellJ, i);
1629  // set it in the matrix
1630  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &val, INSERT_VALUES);
1631  }
1632  }
1633 
1634  return;
1635 }
1636 
1637 void DAJacCon::addConMatCellFaces(
1638  Mat conMat,
1639  const label gRow,
1640  const label cellI,
1641  const word stateName,
1642  const PetscScalar val)
1643 {
1644 
1645  /*
1646  Description:
1647  Insert a value (val) to the connectivity Matrix (conMat)
1648  This value will be inserted at rowI=gRow
1649  The column index is dependent on the cellI's faces and stateName
1650 
1651  Input:
1652  gRow: which row to insert the value for conMat
1653 
1654  cellI: the index of the cell to compute the column index to add
1655 
1656  stateName: the name of the state variable to compute the column index to add
1657 
1658  val: the value to add to conMat
1659 
1660  Output:
1661  conMat: the matrix to add value to
1662 
1663  Example:
1664 
1665  If we want to add value 10.0 to conMat for
1666  columns={the phi globalAdjointIndice of cellI's faces} where cellI=5
1667  row = gRow = 100
1668  Then, call addConMatCell(conMat, 100, 5, "U", 1.0)
1669 
1670  -------
1671  | cellI | <----------- value 10.0 will be added to
1672  ------- column = {global adjoint index
1673  of all cellI's faces}
1674  */
1675 
1676  PetscInt idxJ, idxI;
1677  idxI = gRow;
1678 
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  forAll(faces, idx)
1683  {
1684  //get the appropriate index for this face
1685  label globalState = daIndex_.getGlobalAdjointStateIndex(stateName, faces[idx]);
1686  idxJ = globalState;
1687  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &val, INSERT_VALUES);
1688  }
1689 
1690  return;
1691 }
1692 
1693 void DAJacCon::addBoundaryFaceConnections(
1694  Mat conMat,
1695  const label gRow,
1696  const label cellI,
1697  const labelList v,
1698  const List<List<word>> connectedStates,
1699  const label addFaces)
1700 {
1701  /*
1702  Description:
1703  This function adds inter-proc connectivity into conMat.
1704  For all the inter-proc faces owned by cellI, get the global adj state indices
1705  from DAJacCon::stateBoundaryCon_ and then add them into conMat
1706  Col index to add: the same col index for a given row (bRowGlobal) in the stateBoundaryCon
1707  mat if the element value in the stateBoundaryCon mat is less than the input level,
1708  i.e., v.size().
1709 
1710  Input:
1711  gRow: Row index to add
1712 
1713  cellI: the cell index for getting the faces to add inter-proc connectivity, NOTE: depending on the level
1714  of requested connections, we may add inter-proc face that are not belonged to cellI
1715 
1716  v: an array denoting the desired values to add, the size of v denotes the maximal levels to add
1717 
1718  connectedStates: selectively add some states into the conMat for the current level. If its size is 0,
1719  add all the possible states (except for surfaceStates). The dimension of connectedStates is nLevel
1720  by nStates.
1721 
1722  addFaces: whether to add indices for face (phi) connectivity
1723 
1724  Example:
1725 
1726  labelList val2={1,2};
1727  PetscInt gRow=1024, idxN = 100, addFaces=1;
1728  wordListList connectedStates={{"U","p"},{"U"}};
1729  addBoundaryFaceConnections(stateBoundaryCon,gRow,idxN,vals2,connectedStates,addFaces);
1730  The above call will add 2 levels of connected states for all the inter-proc faces belonged to cellI=idxN
1731  The cols to added are: the level1 connected states (U, p) for all the inter-proc faces belonged to
1732  cellI=idxN. the level2 connected states (U only) for all the inter-proc faces belonged to cellI=idxN
1733  The valus "1" will be added to conMat for all the level1 connected states while the value "2" will be
1734  added for level2.
1735  Note: this function will also add all faces belonged to level1 of the inter-proc faces, see the
1736  following for reference
1737 
1738  -------
1739  | idxN|
1740  | | proc0, idxN=100, globalBFaceI=1024 for the south face of idxN
1741  ----------------- <----coupled boundary face
1742  add state U and p -----> | lv1 | proc1
1743  also add faces --------> | |
1744  -------
1745  | lv2 |
1746  add state U ----------> | |
1747  -------
1748 
1749 
1750  ****** NOTE: *******
1751  If the inter-proc boundary is like the following, calling this function will NOT add any
1752  inter-proc connection for idxN because there is no inter-proc boundary for cell idxN
1753 
1754  -------
1755  | idxN|
1756  | |
1757  ------- <--------- there is no inter-proc boundary for idxN, not adding anything
1758  | lv1 |
1759  | | proc0
1760  ------------------- <----coupled boundary face
1761  | lv2 | proc1
1762  | |
1763  -------
1764 
1765  */
1766 
1767  if (v.size() != connectedStates.size() && connectedStates.size() != 0)
1768  {
1769  FatalErrorIn("") << "size of v and connectedStates are not identical!"
1770  << abort(FatalError);
1771  }
1772 
1773  PetscInt idxJ, idxI, bRow, bRowGlobal;
1774  PetscInt nCols;
1775  const PetscInt* cols;
1776  const PetscScalar* vals;
1777 
1778  PetscInt nColsID;
1779  const PetscInt* colsID;
1780  const PetscScalar* valsID;
1781 
1782  // convert stateNames to stateIDs
1783  labelListList connectedStateIDs(connectedStates.size());
1784  forAll(connectedStates, idxI)
1785  {
1786  forAll(connectedStates[idxI], idxJ)
1787  {
1788  word stateName = connectedStates[idxI][idxJ];
1789  label stateID = daIndex_.adjStateID[stateName];
1790  connectedStateIDs[idxI].append(stateID);
1791  }
1792  }
1793 
1794  idxI = gRow;
1795  // get the faces connected to this cell, note these are in a single
1796  // list that includes all internal and boundary faces
1797  const labelList& faces = mesh_.cells()[cellI];
1798 
1799  //get the level
1800  label level = v.size();
1801 
1802  for (label lv = level; lv >= 1; lv--) // we need to start from the largest levels since they have higher priority
1803  {
1804  forAll(faces, faceI)
1805  {
1806  // Now deal with coupled faces
1807  label currFace = faces[faceI];
1808 
1809  if (daIndex_.isCoupledFace[currFace])
1810  {
1811  //this is a coupled face
1812 
1813  // use the boundary connectivity to figure out what is connected
1814  // to this face for this level
1815 
1816  // get bRow in boundaryCon for this face
1817  bRow = this->getLocalCoupledBFaceIndex(currFace);
1818 
1819  // get the global bRow index
1820  bRowGlobal = daIndex_.globalCoupledBFaceNumbering.toGlobal(bRow);
1821 
1822  // now extract the boundaryCon row
1823  MatGetRow(stateBoundaryCon_, bRowGlobal, &nCols, &cols, &vals);
1824  if (connectedStates.size() != 0)
1825  {
1826  // check if we need to get stateID
1827  MatGetRow(stateBoundaryConID_, bRowGlobal, &nColsID, &colsID, &valsID);
1828  }
1829 
1830  // now loop over the row and set any column that match this level
1831  // in conMat
1832  for (label i = 0; i < nCols; i++)
1833  {
1834  idxJ = cols[i];
1835  label val = round(vals[i]); // val is the connectivity level extracted from stateBoundaryCon_ at this col
1836  // selectively add some states into conMat
1837  label addState;
1838  label stateID = -9999;
1839  // check if we need to get stateID
1840  if (connectedStates.size() != 0)
1841  {
1842  stateID = round(valsID[i]);
1843  }
1844 
1845  if (connectedStates.size() == 0)
1846  {
1847  addState = 1;
1848  }
1849  else if (connectedStateIDs[lv - 1].found(stateID))
1850  {
1851  addState = 1;
1852  }
1853  else
1854  {
1855  addState = 0;
1856  }
1857  // if the level match and the state is what you want
1858  if (val == lv && addState)
1859  {
1860  // need to do v[lv-1] here since v is an array with starting index 0
1861  PetscScalar valIn = v[lv - 1];
1862  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &valIn, INSERT_VALUES);
1863  }
1864  if (val == 10 && addFaces)
1865  {
1866  // this is a necessary connection
1867  PetscScalar valIn = v[lv - 1];
1868  MatSetValues(conMat, 1, &idxI, 1, &idxJ, &valIn, INSERT_VALUES);
1869  }
1870  }
1871 
1872  // restore the row of the matrix
1873  MatRestoreRow(stateBoundaryCon_, bRowGlobal, &nCols, &cols, &vals);
1874  if (connectedStates.size() != 0)
1875  {
1876  // check if we need to get stateID
1877  MatRestoreRow(stateBoundaryConID_, bRowGlobal, &nColsID, &colsID, &valsID);
1878  }
1879  }
1880  }
1881  }
1882 
1883  return;
1884 }
1885 
1886 label DAJacCon::coloringExists(const word postFix) const
1887 {
1888  /*
1889  Description:
1890  Check whether the coloring file exists
1891 
1892  Input:
1893  postFix: the post fix of the file name, e.g., the original
1894  name is dFdWColoring_1.bin, then the new name is
1895  dFdWColoring_drag_1.bin with postFix = _drag
1896 
1897  Output:
1898  return 1 if coloring files exist, otherwise, return 0
1899  */
1900 
1901  Info << "Checking if Coloring file exists.." << endl;
1902  label nProcs = Pstream::nProcs();
1903  word fileName = modelType_ + "Coloring" + postFix + "_" + Foam::name(nProcs);
1904  word checkFile = fileName + ".bin";
1905  std::ifstream fIn(checkFile);
1906  if (!fIn.fail())
1907  {
1908  Info << checkFile << " exists." << endl;
1909  return 1;
1910  }
1911  else
1912  {
1913  return 0;
1914  }
1915 }
1916 
1917 void DAJacCon::calcJacConColoring(const word postFix)
1918 {
1919  /*
1920  Description:
1921  Calculate the coloring for jacCon.
1922 
1923  Input:
1924  postFix: the post fix of the file name, e.g., the original
1925  name is dFdWColoring_1.bin, then the new name is
1926  dFdWColoring_drag_1.bin with postFix = _drag
1927 
1928  Output:
1929  jacConColors_: jacCon coloring and save to files.
1930  The naming convention for coloring vector is
1931  coloringVecName_nProcs.bin. This is necessary because
1932  using different CPU cores result in different jacCon
1933  and therefore different coloring
1934 
1935  nJacColors: number of jacCon colors
1936 
1937  */
1938 
1939  // first check if the file name exists, if yes, return and
1940  // don't compute the coloring
1941  Info << "Calculating " << modelType_ << " Coloring.." << endl;
1942  label nProcs = Pstream::nProcs();
1943  word fileName = modelType_ + "Coloring" + postFix + "_" + Foam::name(nProcs);
1944 
1945  VecZeroEntries(jacConColors_);
1946  if (daOption_.getOption<label>("adjUseColoring"))
1947  {
1948  // use parallelD2 coloring to compute colors
1949  daColoring_.parallelD2Coloring(jacCon_, jacConColors_, nJacConColors_);
1950  }
1951  else
1952  {
1953  // use brute force to compute colors, basically, we assign
1954  // color to its global index
1955  PetscScalar* jacConColorsArray;
1956  VecGetArray(jacConColors_, &jacConColorsArray);
1957  label Istart, Iend;
1958  VecGetOwnershipRange(jacConColors_, &Istart, &Iend);
1959  for (label i = Istart; i < Iend; i++)
1960  {
1961  label relIdx = i - Istart;
1962  jacConColorsArray[relIdx] = i * 1.0;
1963  }
1964  VecRestoreArray(jacConColors_, &jacConColorsArray);
1965 
1966  PetscReal maxVal;
1967  VecMax(jacConColors_, NULL, &maxVal);
1968  nJacConColors_ = maxVal + 1;
1969  }
1970 
1971  daColoring_.validateColoring(jacCon_, jacConColors_);
1972  Info << " nJacConColors: " << nJacConColors_ << endl;
1973  // write jacCon colors
1974  Info << "Writing Colors to " << fileName << endl;
1975  DAUtility::writeVectorBinary(jacConColors_, fileName);
1976 
1977  return;
1978 }
1979 
1980 void DAJacCon::readJacConColoring(const word postFix)
1981 {
1982  /*
1983  Description:
1984  Read the jacCon coloring from files and
1985  compute nJacConColors. The naming convention for
1986  coloring vector is coloringVecName_nProcs.bin
1987  This is necessary because using different CPU
1988  cores result in different jacCon and therefore
1989  different coloring
1990 
1991  Input:
1992  postFix: the post fix of the file name, e.g., the original
1993  name is dFdWColoring_1.bin, then the new name is
1994  dFdWColoring_drag_1.bin with postFix = _drag
1995 
1996  Output:
1997  jacConColors_: read from file
1998 
1999  nJacConColors: number of jacCon colors
2000  */
2001 
2002  label nProcs = Pstream::nProcs();
2003  word fileName = modelType_ + "Coloring" + postFix + "_" + Foam::name(nProcs);
2004  Info << "Reading Coloring " << fileName << endl;
2005 
2006  VecZeroEntries(jacConColors_);
2007  DAUtility::readVectorBinary(jacConColors_, fileName);
2008 
2009  daColoring_.validateColoring(jacCon_, jacConColors_);
2010 
2011  PetscReal maxVal;
2012  VecMax(jacConColors_, NULL, &maxVal);
2013  nJacConColors_ = maxVal + 1;
2014 }
2015 
2016 void DAJacCon::setupJacConPreallocation(const dictionary& options)
2017 {
2018  /*
2019  Description:
2020  Setup the connectivity mat preallocation vectors:
2021 
2022  dRdWTPreallocOn_
2023  dRdWTPreallocOff_
2024  dRdWPreallocOn_
2025  dRdWPreallocOff_
2026 
2027  Input:
2028  options.stateResConInfo: a hashtable that contains the connectivity
2029  information for dRdW, usually obtained from Foam::DAStateInfo
2030  */
2031 
2032  HashTable<List<List<word>>> stateResConInfo;
2033  options.readEntry<HashTable<List<List<word>>>>("stateResConInfo", stateResConInfo);
2034 
2035  label isPrealloc = 1;
2036  this->setupdRdWCon(stateResConInfo, isPrealloc);
2037 }
2038 
2039 void DAJacCon::setupdRdWCon(
2040  const HashTable<List<List<word>>>& stateResConInfo,
2041  const label isPrealloc)
2042 {
2043  /*
2044  Description:
2045  Calculates the state Jacobian connectivity mat DAJacCon::jacCon_ or
2046  computing the preallocation vectors DAJacCondRdW::dRdWPreallocOn_ and
2047  DAJacCondRdW::dRdWPreallocOff_
2048 
2049  Input:
2050  stateResConInfo: a hashtable that contains the connectivity
2051  information for dRdW, usually obtained from Foam::DAStateInfo
2052 
2053  isPrealloc == 1: calculate the preallocation vectors, else, calculate jacCon_
2054 
2055  Output:
2056  DAJacCondRdW::dRdWPreallocOn_: preallocation vector that stores the number of
2057  on-diagonal conectivity for each row
2058 
2059  DAJacCon::jacCon_: state Jacobian connectivity mat with dimension
2060  sizeAdjStates by sizeAdjStates. jacCon_ has the same non-zero pattern as Jacobian mat.
2061  The difference is that jacCon_ has values one for all non-zero values, so jacCon_
2062  may look like this
2063 
2064  1 1 0 0 1 0
2065  1 1 1 0 0 1
2066  0 1 1 1 0 0
2067  jacCon_ = 1 0 1 1 1 0
2068  0 1 0 1 1 1
2069  0 0 1 0 0 1
2070 
2071  Example:
2072  The way setupJacCon works is that we call the DAJacCon::addStateConnections function
2073  to add connectivity for each row of DAJacCon::jacCon_.
2074 
2075  Here we need to loop over all cellI and add a certain number levels of connected states.
2076  If the connectivity list reads:
2077 
2078  stateResConInfo
2079  {
2080  "URes"
2081  {
2082  {"U", "p", "phi"}, // level 0 connectivity
2083  {"U", "p", "phi"}, // level 1 connectivity
2084  {"U"}, // level 2 connectivity
2085  }
2086  }
2087 
2088  and the cell topology with a inter-proc boundary cen be either of the following:
2089  CASE 1:
2090  ---------
2091  | cellQ |
2092  -----------------------
2093  | cellP | cellJ | cellO | <------ proc1
2094  ------------------------------------------------ <----- inter-processor boundary
2095  | cellT | cellK | cellI | cellL | cellU | <------ proc0
2096  -----------------------------------------
2097  | cellN | cellM | cellR |
2098  ------------------------
2099  | cellS |
2100  ---------
2101 
2102  CASE 2:
2103  ---------
2104  | cellQ | <------ proc1
2105  -------------------------------------------------- <----- inter-processor boundary
2106  | cellP | cellJ | cellO | <------ proc0
2107  -----------------------------------------
2108  | cellT | cellK | cellI | cellL | cellU |
2109  -----------------------------------------
2110  | cellN | cellM | cellR |
2111  ------------------------
2112  | cellS |
2113  ---------
2114 
2115  Then, to add the connectivity correctly, we need to add all levels of connected
2116  states for cellI.
2117  Level 0 connectivity is straightforward becasue we don't need
2118  to provide connectedStatesInterProc
2119 
2120  To add level 1 connectivity, we need to:
2121  set connectedLevelLocal = 1
2122  set connectedStatesLocal = {U, p}
2123  set connectedStatesInterProc = {{U,p}, {U}}
2124  set addFace = 1
2125  NOTE: we need set level 1 and level 2 con in connectedStatesInterProc because the
2126  north face of cellI is a inter-proc boundary and there are two levels of connected
2127  state on the other side of the inter-proc boundary for CASE 1. This is the only chance we
2128  can add all two levels of connected state across the boundary for CASE 1. For CASE 2, we won't
2129  add any level 1 inter-proc states because non of the faces for cellI are inter-proc
2130  faces so calling DAJacCon::addBoundaryFaceConnections for cellI won't add anything
2131 
2132  To add level 2 connectivity, we need to
2133  set connectedLevelLocal = 2
2134  set connectedStatesLocal = {U}
2135  set connectedStatesInterProc = {{U}}
2136  set addFace = 0
2137  NOTE 1: we need only level 2 con (U) for connectedStatesInterProc because if we are in CASE 1,
2138  the level 2 of inter-proc states have been added. For CASE 2, we only need to add cellQ
2139  by calling DAJacCon::addBoundaryFaceConnections with cellJ
2140  NOTE 2: If we didn't call two levels of connectedStatesInterProc in the previous call for
2141  level 1 con, we can not add it for connectedLevelLocal = 2 becasue for CASE 2 there is no
2142  inter-proc boundary for cellI
2143 
2144  NOTE: how to provide connectedLevelLocal, connectedStatesLocal, and connectedStatesInterProc
2145  are done in this function
2146 
2147  */
2148 
2149  label globalIdx;
2150  // connectedStatesP: one row matrix that stores the actual connectivity,
2151  // element value 1 denotes a connected state. connectedStatesP is then used to
2152  // assign dRdWPreallocOn or dRdWCon
2153  Mat connectedStatesP;
2154 
2155  PetscInt nCols;
2156  const PetscInt* cols;
2157  const PetscScalar* vals;
2158 
2159  PetscInt nColsID;
2160  const PetscInt* colsID;
2161  const PetscScalar* valsID;
2162 
2163  if (isPrealloc)
2164  {
2165  VecZeroEntries(dRdWPreallocOn_);
2166  VecZeroEntries(dRdWPreallocOff_);
2167  VecZeroEntries(dRdWTPreallocOn_);
2168  VecZeroEntries(dRdWTPreallocOff_);
2169  }
2170 
2171  if (daOption_.getOption<label>("debug"))
2172  {
2173  if (isPrealloc)
2174  {
2175 
2176  Info << "Computing preallocating vectors for Jacobian connectivity mat" << endl;
2177  }
2178  else
2179  {
2180  Info << "Setup Jacobian connectivity mat" << endl;
2181  }
2182  }
2183 
2184  // loop over all cell residuals, we bascially need to compute all
2185  // the input parameters for the DAJacCondRdW::addStateConnections function, then call
2186  // the function to get connectedStatesP (matrix that contains one row of connectivity
2187  // in DAJacCondRdW::jacCon_). Check DAJacCondRdW::addStateConnections for detail usages
2188  forAll(daIndex_.adjStateNames, idxI)
2189  {
2190  // get stateName and residual names
2191  word stateName = daIndex_.adjStateNames[idxI];
2192  word resName = stateName + "Res";
2193 
2194  // check if this state is a cell state, we do surfaceScalarState residuals separately
2195  if (daIndex_.adjStateType[stateName] == "surfaceScalarState")
2196  {
2197  continue;
2198  }
2199 
2200  // maximal connectivity level information
2201  // Note that stateResConInfo starts with level zero,
2202  // so the maxConLeve is its size minus one
2203  label maxConLevel = stateResConInfo[resName].size() - 1;
2204 
2205  // if it is a vectorState, set compMax=3
2206  label compMax = 1;
2207  if (daIndex_.adjStateType[stateName] == "volVectorState")
2208  {
2209  compMax = 3;
2210  }
2211 
2212  forAll(mesh_.cells(), cellI)
2213  {
2214  for (label comp = 0; comp < compMax; comp++)
2215  {
2216 
2217  // zero the connections
2218  this->createConnectionMat(&connectedStatesP);
2219 
2220  // now add the con. We loop over all the connectivity levels
2221  forAll(stateResConInfo[resName], idxJ) // idxJ: con level
2222  {
2223 
2224  // set connectedStatesLocal: the locally connected state variables for this level
2225  wordList connectedStatesLocal(0);
2226  forAll(stateResConInfo[resName][idxJ], idxK)
2227  {
2228  word conName = stateResConInfo[resName][idxJ][idxK];
2229  // Exclude surfaceScalarState when appending connectedStatesLocal
2230  // whether to add it depends on addFace parameter
2231  if (daIndex_.adjStateType[conName] != "surfaceScalarState")
2232  {
2233  connectedStatesLocal.append(conName);
2234  }
2235  }
2236 
2237  // set connectedStatesInterProc: the globally connected state variables for this level
2238  List<List<word>> connectedStatesInterProc;
2239  if (idxJ == 0)
2240  {
2241  // pass a zero list, no need to add interProc connecitivity for level 0
2242  connectedStatesInterProc.setSize(0);
2243  }
2244  else if (idxJ != maxConLevel)
2245  {
2246  connectedStatesInterProc.setSize(maxConLevel - idxJ + 1);
2247  for (label k = 0; k < maxConLevel - idxJ + 1; k++)
2248  {
2249  label conSize = stateResConInfo[resName][k + idxJ].size();
2250  for (label l = 0; l < conSize; l++)
2251  {
2252  word conName = stateResConInfo[resName][k + idxJ][l];
2253  // Exclude surfaceScalarState when appending connectedStatesLocal
2254  // whether to add it depends on addFace parameter
2255  if (daIndex_.adjStateType[conName] != "surfaceScalarState")
2256  {
2257  connectedStatesInterProc[k].append(conName);
2258  }
2259  }
2260  }
2261  }
2262  else
2263  {
2264  connectedStatesInterProc.setSize(1);
2265  label conSize = stateResConInfo[resName][maxConLevel].size();
2266  for (label l = 0; l < conSize; l++)
2267  {
2268  word conName = stateResConInfo[resName][maxConLevel][l];
2269  // Exclude surfaceScalarState when appending connectedStatesLocal
2270  // whether to add it depends on addFace parameter
2271  if (daIndex_.adjStateType[conName] != "surfaceScalarState")
2272  {
2273  connectedStatesInterProc[0].append(conName);
2274  }
2275  }
2276  }
2277 
2278  // check if we need to addFace for this level
2279  label addFace = 0;
2280  forAll(stateInfo_["surfaceScalarStates"], idxK)
2281  {
2282  word conName = stateInfo_["surfaceScalarStates"][idxK];
2283  if (stateResConInfo[resName][idxJ].found(conName))
2284  {
2285  addFace = 1;
2286  }
2287  }
2288 
2289  // special treatment for pressureInletVelocity. No such treatment is needed
2290  // for dFdW because we have added phi in their conInfo
2291  if (daField_.specialBCs.found("pressureInletVelocity")
2292  && this->addPhi4PIV(stateName, cellI, comp))
2293  {
2294  addFace = 1;
2295  }
2296 
2297  // Add connectivity
2298  this->addStateConnections(
2299  connectedStatesP,
2300  cellI,
2301  idxJ,
2302  connectedStatesLocal,
2303  connectedStatesInterProc,
2304  addFace);
2305 
2306  //Info<<"lv: "<<idxJ<<" locaStates: "<<connectedStatesLocal<<" interProcStates: "
2307  // <<connectedStatesInterProc<<" addFace: "<<addFace<<endl;
2308  }
2309 
2310  // get the global index of the current state for the row index
2311  globalIdx = daIndex_.getGlobalAdjointStateIndex(stateName, cellI, comp);
2312 
2313  if (isPrealloc)
2314  {
2315  this->allocateJacobianConnections(
2316  dRdWPreallocOn_,
2317  dRdWPreallocOff_,
2318  dRdWTPreallocOn_,
2319  dRdWTPreallocOff_,
2320  connectedStatesP,
2321  globalIdx);
2322  }
2323  else
2324  {
2325  this->setupJacobianConnections(
2326  jacCon_,
2327  connectedStatesP,
2328  globalIdx);
2329  }
2330  }
2331  }
2332  }
2333 
2334  // loop over all face residuals
2335  forAll(stateInfo_["surfaceScalarStates"], idxI)
2336  {
2337  // get stateName and residual names
2338  word stateName = stateInfo_["surfaceScalarStates"][idxI];
2339  word resName = stateName + "Res";
2340 
2341  // maximal connectivity level information
2342  label maxConLevel = stateResConInfo[resName].size() - 1;
2343 
2344  forAll(mesh_.faces(), faceI)
2345  {
2346 
2347  //zero the connections
2348  this->createConnectionMat(&connectedStatesP);
2349 
2350  // Get the owner and neighbour cells for this face
2351  label idxO = -1, idxN = -1;
2352  if (faceI < daIndex_.nLocalInternalFaces)
2353  {
2354  idxO = mesh_.owner()[faceI];
2355  idxN = mesh_.neighbour()[faceI];
2356  }
2357  else
2358  {
2359  label relIdx = faceI - daIndex_.nLocalInternalFaces;
2360  label patchIdx = daIndex_.bFacePatchI[relIdx];
2361  label faceIdx = daIndex_.bFaceFaceI[relIdx];
2362 
2363  const UList<label>& pFaceCells = mesh_.boundaryMesh()[patchIdx].faceCells();
2364  idxN = pFaceCells[faceIdx];
2365  }
2366 
2367  // now add the con. We loop over all the connectivity levels
2368  forAll(stateResConInfo[resName], idxJ) // idxJ: con level
2369  {
2370 
2371  // set connectedStatesLocal: the locally connected state variables for this level
2372  wordList connectedStatesLocal(0);
2373  forAll(stateResConInfo[resName][idxJ], idxK)
2374  {
2375  word conName = stateResConInfo[resName][idxJ][idxK];
2376  // Exclude surfaceScalarState when appending connectedStatesLocal
2377  // whether to add it depends on addFace parameter
2378  if (daIndex_.adjStateType[conName] != "surfaceScalarState")
2379  {
2380  connectedStatesLocal.append(conName);
2381  }
2382  }
2383 
2384  // set connectedStatesInterProc: the globally connected state variables for this level
2385  List<List<word>> connectedStatesInterProc;
2386  if (idxJ == 0)
2387  {
2388  // pass a zero list, no need to add interProc connecitivity for level 0
2389  connectedStatesInterProc.setSize(0);
2390  }
2391  else if (idxJ != maxConLevel)
2392  {
2393  connectedStatesInterProc.setSize(maxConLevel - idxJ + 1);
2394  for (label k = 0; k < maxConLevel - idxJ + 1; k++)
2395  {
2396  label conSize = stateResConInfo[resName][k + idxJ].size();
2397  for (label l = 0; l < conSize; l++)
2398  {
2399  word conName = stateResConInfo[resName][k + idxJ][l];
2400  // Exclude surfaceScalarState when appending connectedStatesLocal
2401  // whether to add it depends on addFace parameter
2402  if (daIndex_.adjStateType[conName] != "surfaceScalarState")
2403  {
2404  connectedStatesInterProc[k].append(conName);
2405  }
2406  }
2407  }
2408  }
2409  else
2410  {
2411  connectedStatesInterProc.setSize(1);
2412  label conSize = stateResConInfo[resName][maxConLevel].size();
2413  for (label l = 0; l < conSize; l++)
2414  {
2415  word conName = stateResConInfo[resName][maxConLevel][l];
2416  // Exclude surfaceScalarState when appending connectedStatesLocal
2417  // whether to add it depends on addFace parameter
2418  if (daIndex_.adjStateType[conName] != "surfaceScalarState")
2419  {
2420  connectedStatesInterProc[0].append(conName);
2421  }
2422  }
2423  }
2424 
2425  // check if we need to addFace for this level
2426  label addFace = 0;
2427  forAll(stateInfo_["surfaceScalarStates"], idxK)
2428  {
2429  word conName = stateInfo_["surfaceScalarStates"][idxK];
2430  // NOTE: we need special treatment for boundary faces for level>0
2431  // since addFace for boundary face should add one more extra level of faces
2432  // This is because we only have idxN for a boundary face while the idxO can
2433  // be on the other side of the inter-proc boundary
2434  // In this case, we need to use idxJ-1 instead of idxJ information to tell whether to addFace
2435  label levelCheck;
2436  if (faceI < daIndex_.nLocalInternalFaces or idxJ == 0)
2437  {
2438  levelCheck = idxJ;
2439  }
2440  else
2441  {
2442  levelCheck = idxJ - 1;
2443  }
2444 
2445  if (stateResConInfo[resName][levelCheck].found(conName))
2446  {
2447  addFace = 1;
2448  }
2449  }
2450 
2451  // special treatment for pressureInletVelocity. No such treatment is needed
2452  // for dFdW because we have added phi in their conInfo
2453  if (daField_.specialBCs.found("pressureInletVelocity")
2454  && this->addPhi4PIV(stateName, faceI))
2455  {
2456  addFace = 1;
2457  }
2458 
2459  // Add connectivity for idxN
2460  this->addStateConnections(
2461  connectedStatesP,
2462  idxN,
2463  idxJ,
2464  connectedStatesLocal,
2465  connectedStatesInterProc,
2466  addFace);
2467 
2468  if (faceI < daIndex_.nLocalInternalFaces)
2469  {
2470  // Add connectivity for idxO
2471  this->addStateConnections(
2472  connectedStatesP,
2473  idxO,
2474  idxJ,
2475  connectedStatesLocal,
2476  connectedStatesInterProc,
2477  addFace);
2478  }
2479 
2480  //Info<<"lv: "<<idxJ<<" locaStates: "<<connectedStatesLocal<<" interProcStates: "
2481  // <<connectedStatesInterProc<<" addFace: "<<addFace<<endl;
2482  }
2483 
2484  // NOTE: if this faceI is on a coupled patch, the above connectivity is not enough to
2485  // cover the points on the other side of proc domain, we need to add 3 lvs of cells here
2486  if (faceI >= daIndex_.nLocalInternalFaces)
2487  {
2488  label relIdx = faceI - daIndex_.nLocalInternalFaces;
2489  label patchIdx = daIndex_.bFacePatchI[relIdx];
2490 
2491  label maxLevel = stateResConInfo[resName].size();
2492 
2493  if (mesh_.boundaryMesh()[patchIdx].coupled())
2494  {
2495 
2496  label bRow = this->getLocalCoupledBFaceIndex(faceI);
2497  label bRowGlobal = daIndex_.globalCoupledBFaceNumbering.toGlobal(bRow);
2498  MatGetRow(stateBoundaryCon_, bRowGlobal, &nCols, &cols, &vals);
2499  MatGetRow(stateBoundaryConID_, bRowGlobal, &nColsID, &colsID, &valsID);
2500  for (label i = 0; i < nCols; i++)
2501  {
2502  PetscInt idxJ = cols[i];
2503  label val = round(vals[i]);
2504  // we are going to add some selective states with connectivity level <= 3
2505  // first check the state
2506  label stateID = round(valsID[i]);
2507  word conName = daIndex_.adjStateNames[stateID];
2508  label addState = 0;
2509  // NOTE: we use val-1 here since phi actually has 3 levels of connectivity
2510  // however, when we assign stateResConInfo, we ignore the level 0
2511  // connectivity since they are idxN and idxO
2512  if (val != 10 && val < maxLevel + 1)
2513  {
2514  if (stateResConInfo[resName][val - 1].found(conName))
2515  {
2516  addState = 1;
2517  }
2518  }
2519  if (addState == 1 && val < maxLevel + 1 && val > 0)
2520  {
2521  this->setConnections(connectedStatesP, idxJ);
2522  }
2523  }
2524  MatRestoreRow(stateBoundaryCon_, bRowGlobal, &nCols, &cols, &vals);
2525  MatRestoreRow(stateBoundaryConID_, bRowGlobal, &nColsID, &colsID, &valsID);
2526  }
2527  }
2528 
2529  // get the global index of the current state for the row index
2530  globalIdx = daIndex_.getGlobalAdjointStateIndex(stateName, faceI);
2531 
2532  if (isPrealloc)
2533  {
2534  this->allocateJacobianConnections(
2535  dRdWPreallocOn_,
2536  dRdWPreallocOff_,
2537  dRdWTPreallocOn_,
2538  dRdWTPreallocOff_,
2539  connectedStatesP,
2540  globalIdx);
2541  }
2542  else
2543  {
2544  this->setupJacobianConnections(
2545  jacCon_,
2546  connectedStatesP,
2547  globalIdx);
2548  }
2549  }
2550  }
2551 
2552  if (isPrealloc)
2553  {
2554  VecAssemblyBegin(dRdWPreallocOn_);
2555  VecAssemblyEnd(dRdWPreallocOn_);
2556  VecAssemblyBegin(dRdWPreallocOff_);
2557  VecAssemblyEnd(dRdWPreallocOff_);
2558  VecAssemblyBegin(dRdWTPreallocOn_);
2559  VecAssemblyEnd(dRdWTPreallocOn_);
2560  VecAssemblyBegin(dRdWTPreallocOff_);
2561  VecAssemblyEnd(dRdWTPreallocOff_);
2562 
2563  //output the matrix to a file
2564  wordList writeJacobians;
2565  daOption_.getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
2566  if (writeJacobians.found("dRdWTPrealloc"))
2567  {
2568  DAUtility::writeVectorASCII(dRdWTPreallocOn_, "dRdWTPreallocOn");
2569  DAUtility::writeVectorASCII(dRdWTPreallocOff_, "dRdWTPreallocOff");
2570  DAUtility::writeVectorASCII(dRdWPreallocOn_, "dRdWPreallocOn");
2571  DAUtility::writeVectorASCII(dRdWPreallocOff_, "dRdWPreallocOff");
2572  }
2573  }
2574  else
2575  {
2576  MatAssemblyBegin(jacCon_, MAT_FINAL_ASSEMBLY);
2577  MatAssemblyEnd(jacCon_, MAT_FINAL_ASSEMBLY);
2578 
2579  //output the matrix to a file
2580  wordList writeJacobians;
2581  daOption_.getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
2582  if (writeJacobians.found("dRdWCon"))
2583  {
2584  //DAUtility::writeMatRowSize(jacCon_, "dRdWCon");
2585  DAUtility::writeMatrixBinary(jacCon_, "dRdWCon");
2586  }
2587  }
2588 
2589  if (daOption_.getOption<label>("debug"))
2590  {
2591  if (isPrealloc)
2592  {
2593  Info << "Preallocating state Jacobian connectivity mat: finished!" << endl;
2594  }
2595  else
2596  {
2597  Info << "Setup state Jacobian connectivity mat: finished!" << endl;
2598  }
2599  }
2600 }
2601 
2602 void DAJacCon::allocateJacobianConnections(
2603  Vec preallocOnProc,
2604  Vec preallocOffProc,
2605  Vec preallocOnProcT,
2606  Vec preallocOffProcT,
2607  Mat connections,
2608  const label row)
2609 {
2610  /*
2611  Description:
2612  Compute the matrix allocation vector based on one row connection mat
2613 
2614  Input:
2615  connections: a one row matrix that contains all the nonzeros for one
2616  row of DAJacCon::jacCon_
2617 
2618  row: which row to add for the preallocation vector
2619 
2620  Output:
2621  preallocOnProc: the vector that contains the number of nonzeros for each row
2622  in dRdW (on-diagonal block elements)
2623 
2624  preallocOffProc: the vector that contains the number of nonzeros for each row
2625  in dRdW (off-diagonal block elements)
2626 
2627  preallocOnProcT: the vector that contains the number of nonzeros for each row
2628  in dRdWT (on-diagonal block elements)
2629 
2630  preallocOffProcT: the vector that contains the number of nonzeros for each row
2631  in dRdWT (off-diagonal block elements)
2632 
2633  */
2634  PetscScalar v = 1.0;
2635  PetscInt nCols;
2636  const PetscInt* cols;
2637  const PetscScalar* vals;
2638 
2639  MatAssemblyBegin(connections, MAT_FINAL_ASSEMBLY);
2640  MatAssemblyEnd(connections, MAT_FINAL_ASSEMBLY);
2641 
2642  // Compute the transposed case
2643  // in this case connections represents a single column, so we need to
2644  // increment the counter in each row with a non-zero entry.
2645 
2646  label colMin = daIndex_.globalAdjointStateNumbering.toGlobal(0);
2647  label colMax = colMin + daIndex_.nLocalAdjointStates;
2648  // by construction rows should be limited to local rows
2649  MatGetRow(connections, 0, &nCols, &cols, &vals);
2650 
2651  // for the non-transposed case just sum up the row.
2652  // count up the total number of non zeros in this row
2653  label totalCount = 0; //2
2654  label localCount = 0;
2655  //int idx;
2656  for (label j = 0; j < nCols; j++)
2657  {
2658  // int idx = cols[j];
2659  scalar val = vals[j];
2660  if (DAUtility::isValueCloseToRef(val, 1.0))
2661  {
2662  // We can compute the first part of the non-transposed row here.
2663  totalCount++;
2664  label idx = cols[j];
2665  // Set the transposed version as well
2666  if (colMin <= idx && idx < colMax)
2667  {
2668  //this entry is a local entry, increment the corresponding row
2669  VecSetValue(preallocOnProcT, idx, v, ADD_VALUES);
2670  localCount++;
2671  }
2672  else
2673  {
2674  // this is an off proc entry.
2675  VecSetValue(preallocOffProcT, idx, v, ADD_VALUES);
2676  }
2677  }
2678  }
2679 
2680  label offProcCount = totalCount - localCount;
2681  VecSetValue(preallocOnProc, row, localCount, INSERT_VALUES);
2682  VecSetValue(preallocOffProc, row, offProcCount, INSERT_VALUES);
2683 
2684  // restore the row of the matrix
2685  MatRestoreRow(connections, 0, &nCols, &cols, &vals);
2686  MatDestroy(&connections);
2687 
2688  return;
2689 }
2690 
2691 void DAJacCon::calcColoredColumns(
2692  const label colorI,
2693  Vec coloredColumn) const
2694 {
2695  /*
2696  Description:
2697  Compute the colored column vector: coloredColumn. This vector will then
2698  be used to assign resVec to dRdW in DAPartDeriv::calcPartDeriv
2699 
2700  Input:
2701  colorI: the ith color index
2702 
2703  Output:
2704  coloredColumn: For a given colorI, coloredColumn vector contains the column
2705  index for non-zero elements in the DAJacCon::jacCon_ matrix. If there is
2706  no non-zero element for this color, set the value to -1
2707 
2708  Example:
2709 
2710  If the DAJacCon::jacCon_ matrix reads,
2711 
2712  color0 color1
2713  | |
2714  1 0 0 0
2715  jacCon = 0 1 1 0
2716  0 0 1 0
2717  0 0 0 1
2718  | |
2719  color0 color0
2720 
2721  and the coloring vector DAJacCon::jacConColors_ = {0, 0, 1, 0}.
2722 
2723  **************************
2724  ***** If colorI = 0 ******
2725  **************************
2726  Calling calcColoredColumns(0, coloredColumn) will return
2727 
2728  coloredColumn = {0, 1, -1, 3}
2729 
2730  In this case, we have three columns (0, 1, and 3) for color=0, the nonzero pattern is:
2731 
2732  color0
2733  |
2734  1 0 0 0
2735  jacCon = 0 1 0 0
2736  0 0 0 0
2737  0 0 0 1
2738  | |
2739  color0 color0
2740 
2741  So for the 0th, 1th, and 3th rows, the non-zero elements in the jacCon matrix are at the
2742  0th, 1th, and 3th columns, respectively, which gives coloredColumn = {0, 1, -1, 3}
2743 
2744  **************************
2745  ***** If colorI = 1 ******
2746  **************************
2747  Calling calcColoredColumns(1, coloredColumn) will return
2748 
2749  coloredColumn = {-1, 2, 2, -1}
2750 
2751  In this case, we have one column (2) for color=1, , the nonzero pattern is:
2752 
2753  color1
2754  |
2755  0 0 0 0
2756  jacCon = 0 0 1 0
2757  0 0 1 0
2758  0 0 0 0
2759 
2760  So for the 1th and 2th rows, the non-zero elements in the jacCon matrix are at the
2761  2th and 2th columns, respectively, which gives coloredColumn = {-1, 2, 2, -1}
2762 
2763  */
2764 
2765  if (daOption_.getOption<label>("adjUseColoring"))
2766  {
2767 
2768  Vec colorIdx;
2769  label Istart, Iend;
2770 
2771  /* for the current color, determine which row/column pairs match up. */
2772 
2773  // create a vector to hold the column indices associated with this color
2774  VecDuplicate(jacConColors_, &colorIdx);
2775  VecZeroEntries(colorIdx);
2776 
2777  // Start by looping over the color vector. Set each column index associated
2778  // with the current color to its own value in the color idx vector
2779  // get the values on this proc
2780  VecGetOwnershipRange(jacConColors_, &Istart, &Iend);
2781 
2782  // create the arrays to access them directly
2783  const PetscScalar* jacConColorsArray;
2784  PetscScalar* colorIdxArray;
2785  VecGetArrayRead(jacConColors_, &jacConColorsArray);
2786  VecGetArray(colorIdx, &colorIdxArray);
2787 
2788  // loop over the entries to find the ones that match this color
2789  for (label j = Istart; j < Iend; j++)
2790  {
2791  label idx = j - Istart;
2792  if (DAUtility::isValueCloseToRef(jacConColorsArray[idx], colorI * 1.0))
2793  {
2794  // use 1 based indexing here and then subtract 1 from all values in
2795  // the mat mult. This will handle the zero index case in the first row
2796  colorIdxArray[idx] = j + 1;
2797  }
2798  }
2799  VecRestoreArrayRead(jacConColors_, &jacConColorsArray);
2800  VecRestoreArray(colorIdx, &colorIdxArray);
2801 
2802  //VecAssemblyBegin(colorIdx);
2803  //VecAssemblyEnd(colorIdx);
2804 
2805  //Set coloredColumn to -1 to account for the 1 based indexing in the above loop
2806  VecSet(coloredColumn, -1);
2807  // Now do a MatVec with the conMat to get the row coloredColumn pairs.
2808  MatMultAdd(jacCon_, colorIdx, coloredColumn, coloredColumn);
2809 
2810  // destroy the temporary vector
2811  VecDestroy(&colorIdx);
2812  }
2813  else
2814  {
2815  // uncolored case, we just set all elements to colorI
2816  PetscScalar val = colorI * 1.0;
2817  VecSet(coloredColumn, val);
2818 
2819  VecAssemblyBegin(coloredColumn);
2820  VecAssemblyEnd(coloredColumn);
2821  }
2822 }
2823 
2824 void DAJacCon::checkSpecialBCs()
2825 {
2826  /*
2827  Description:
2828  Check if there is special boundary conditions that need special treatment in jacCon_
2829  */
2830 
2831  // *******************************************************************
2832  // pressureInletVelocity
2833  // *******************************************************************
2834  if (daField_.specialBCs.found("pressureInletVelocity"))
2835  {
2836  // we need special treatment for connectivity
2837  Info << "pressureInletVelocity detected, applying special treatment for connectivity." << endl;
2838  // initialize and compute isPIVBCState_
2839  VecCreate(PETSC_COMM_WORLD, &isPIVBCState_);
2840  VecSetSizes(isPIVBCState_, daIndex_.nLocalAdjointStates, PETSC_DECIDE);
2841  VecSetFromOptions(isPIVBCState_);
2842  VecZeroEntries(isPIVBCState_);
2843 
2844  // compute isPIVBCState_
2845  // Note we need to read the U field, instead of getting it from db
2846  // this is because coloringSolver does not read U
2847  Info << "Calculating pressureInletVelocity state vec..." << endl;
2848  volVectorField U(
2849  IOobject(
2850  "U",
2851  mesh_.time().timeName(),
2852  mesh_,
2853  IOobject::READ_IF_PRESENT,
2854  IOobject::NO_WRITE,
2855  false),
2856  mesh_,
2857  dimensionedVector("U", dimensionSet(0, 1, -1, 0, 0, 0, 0), vector::zero),
2858  zeroGradientFvPatchField<vector>::typeName);
2859 
2860  forAll(U.boundaryField(), patchI)
2861  {
2862  if (U.boundaryFieldRef()[patchI].type() == "pressureInletVelocity")
2863  {
2864  const UList<label>& pFaceCells = mesh_.boundaryMesh()[patchI].faceCells();
2865  forAll(U.boundaryFieldRef()[patchI], faceI)
2866  {
2867 
2868  label faceCellI = pFaceCells[faceI]; // lv 1 face neibouring cells
2869  this->setPIVVec(isPIVBCState_, faceCellI);
2870 
2871  forAll(mesh_.cellCells()[faceCellI], cellJ)
2872  {
2873 
2874  label faceCellJ = mesh_.cellCells()[faceCellI][cellJ]; // lv 2 face neibouring cells
2875  this->setPIVVec(isPIVBCState_, faceCellJ);
2876 
2877  forAll(mesh_.cellCells()[faceCellJ], cellK)
2878  {
2879  label faceCellK = mesh_.cellCells()[faceCellI][cellJ]; // lv 3 face neibouring cells
2880  this->setPIVVec(isPIVBCState_, faceCellK);
2881  }
2882  }
2883  }
2884  }
2885  }
2886 
2887  VecAssemblyBegin(isPIVBCState_);
2888  VecAssemblyEnd(isPIVBCState_);
2889 
2890  wordList writeJacobians;
2891  daOption_.getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
2892  if (writeJacobians.found("isPIVVec"))
2893  {
2894  DAUtility::writeVectorASCII(isPIVBCState_, "isPIVVec");
2895  }
2896  }
2897 
2898  // *******************************************************************
2899  // append more special BCs
2900  // *******************************************************************
2901 
2902  return;
2903 }
2904 
2905 void DAJacCon::setPIVVec(
2906  Vec isPIV,
2907  const label cellI)
2908 {
2909  forAll(daIndex_.adjStateNames, idxI)
2910  {
2911  // get stateName and residual names
2912  word stateName = daIndex_.adjStateNames[idxI];
2913 
2914  // check if this state is a cell state, we do surfaceScalarState separately
2915  if (daIndex_.adjStateType[stateName] == "surfaceScalarState")
2916  {
2917  forAll(mesh_.cells()[cellI], faceI)
2918  {
2919  label cellFaceI = mesh_.cells()[cellI][faceI];
2920  label rowI = daIndex_.getGlobalAdjointStateIndex(stateName, cellFaceI);
2921  VecSetValue(isPIV, rowI, 1.0, INSERT_VALUES);
2922  }
2923  }
2924  else
2925  {
2926  // if it is a vectorState, set compMax=3
2927  label compMax = 1;
2928  if (daIndex_.adjStateType[stateName] == "volVectorState")
2929  {
2930  compMax = 3;
2931  }
2932 
2933  for (label comp = 0; comp < compMax; comp++)
2934  {
2935  label rowI = daIndex_.getGlobalAdjointStateIndex(stateName, cellI, comp);
2936  VecSetValue(isPIV, rowI, 1.0, INSERT_VALUES);
2937  }
2938  }
2939  }
2940 }
2941 
2942 label DAJacCon::addPhi4PIV(
2943  const word stateName,
2944  const label idxI,
2945  const label comp)
2946 {
2947 
2948  label localI = daIndex_.getLocalAdjointStateIndex(stateName, idxI, comp);
2949  PetscScalar* isPIVArray;
2950  VecGetArray(isPIVBCState_, &isPIVArray);
2951  if (DAUtility::isValueCloseToRef(isPIVArray[localI], 1.0))
2952  {
2953  VecRestoreArray(isPIVBCState_, &isPIVArray);
2954  return 1;
2955  }
2956  else
2957  {
2958  VecRestoreArray(isPIVBCState_, &isPIVArray);
2959  return 0;
2960  }
2961 
2962  VecRestoreArray(isPIVBCState_, &isPIVArray);
2963  return 0;
2964 }
2965 
2966 void DAJacCon::clear()
2967 {
2968  /*
2969  Description:
2970  Clear all members to avoid memory leak because we will initalize
2971  multiple objects of DAJacCon. Here we need to delete all members
2972  in the parent and child classes
2973  */
2974 
2975  VecDestroy(&dRdWTPreallocOn_);
2976  VecDestroy(&dRdWTPreallocOff_);
2977  VecDestroy(&dRdWPreallocOn_);
2978  VecDestroy(&dRdWPreallocOff_);
2979 
2980  this->clearDAJacConMembers();
2981 }
2982 
2983 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2984 
2985 } // End namespace Foam
2986 
2987 // ************************************************************************* //
Foam::DAJacCon::modelType_
const word modelType_
the name of the jacCon matrix
Definition: DAJacCon.H:47
Foam::DAJacCon::initializeJacCon
void initializeJacCon(const dictionary &options)
initialize the state Jacobian connectivity matrix
Definition: DAJacCon.C:254
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Foam::DAJacCon::daColoring_
DAColoring daColoring_
DAColoring object.
Definition: DAJacCon.H:62
Foam::DAJacCon::setPIVVec
void setPIVVec(Vec iSPIV, const label idxI)
function used to add connectivity for pressureInletVelocity
Definition: DAJacCon.C:2905
Foam::DAUtility::readVectorBinary
static void readVectorBinary(Vec vecIn, const word prefix)
read petsc vector in binary format
Definition: DAUtility.C:282
Foam::DAStateInfo::New
static autoPtr< DAStateInfo > New(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel)
Definition: DAStateInfo.C:43
Foam::DAIndex::globalAdjointStateNumbering
globalIndex globalAdjointStateNumbering
Definition: DAIndex.H:167
Foam::DAJacCon::checkSpecialBCs
void checkSpecialBCs()
check if there is special boundary conditions that need special treatment in jacCon_
Definition: DAJacCon.C:2824
Foam::DAJacCon::calcColoredColumns
void calcColoredColumns(const label colorI, Vec coloredColumn) const
calculate the colored column vector
Definition: DAJacCon.C:2691
Foam::DAJacCon::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAJacCon.H:50
Foam::DAJacCon::jacCon_
Mat jacCon_
Jacobian connectivity mat.
Definition: DAJacCon.H:80
Foam::DAIndex::bFacePatchI
labelList bFacePatchI
Definition: DAIndex.H:110
Foam::DAJacCon::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAJacCon.H:53
Foam::DAJacCon::clearDAJacConMembers
void clearDAJacConMembers()
clear members in DAJacCon
Definition: DAJacCon.H:173
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:2942
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::setupStateBoundaryConID
void setupStateBoundaryConID(Mat *stateBoundaryConID)
calculate DAJacCon::stateBoundaryConID_
Definition: DAJacCon.C:1443
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:108
Foam::DAJacCon::setupJacCon
void setupJacCon(const dictionary &options)
assign 1 to all non-zero elements for the Jacobian connecitivyt matrix
Definition: DAJacCon.C:286
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:1207
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:753
Foam::DAOption::getOption
classType getOption(const word key) const
get an option from subDict and key
Definition: DAOption.H:92
Foam::DAJacCon::dRdWPreallocOff_
Vec dRdWPreallocOff_
Definition: DAJacCon.H:93
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:45
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:111
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:707
Foam::DAJacCon::calcNeiBFaceGlobalCompact
void calcNeiBFaceGlobalCompact(labelList &neiBFaceGlobalCompact)
calculate DAJacCon::neiBFaceGlobalCompact_
Definition: DAJacCon.C:686
solverName
word solverName
Definition: createAdjoint.H:14
Foam::DAJacCon::clear
void clear()
clear members in parent and child objects
Definition: DAJacCon.C:2966
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::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:1509
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:1568
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:315
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
Foam::DAJacCon::allocateJacobianConnections
void allocateJacobianConnections(Vec preallocOnProc, Vec preallocOffProc, Vec preallocOnProcT, Vec preallocOffProcT, Mat connections, const label row)
compute preallocation vectors
Definition: DAJacCon.C:2602
dafoam_plot3dtransform.vals
vals
Definition: dafoam_plot3dtransform.py:39
Foam::DAJacCon::preallocatedRdW
void preallocatedRdW(Mat dRMat, const label transposed) const
preallocate dRdW matrix using the preallocVec
Definition: DAJacCon.C:228
Foam::DAModel
Definition: DAModel.H:57
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:188
k
dimensionedScalar & k
Definition: createRefsHeatTransfer.H:6
Foam
Definition: checkGeometry.C:32
Foam::DAJacCon::dRdWTPreallocOn_
Vec dRdWTPreallocOn_
Definition: DAJacCon.H:90
Foam::DAJacCon::stateBoundaryCon_
Mat stateBoundaryCon_
matrix to store boundary connectivity levels for state Jacobians
Definition: DAJacCon.H:71
Foam::DAJacCon::daField_
DAField daField_
DAField object.
Definition: DAJacCon.H:65
Foam::DAJacCon::daIndex_
const DAIndex & daIndex_
DAIndex object.
Definition: DAJacCon.H:59
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
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:1886
Foam::DAUtility::writeMatrixBinary
static void writeMatrixBinary(const Mat matIn, const word prefix)
write petsc matrix in binary format
Definition: DAUtility.C:411
Foam::DAJacCon::setupStateBoundaryCon
void setupStateBoundaryCon(Mat *stateBoundaryCon)
calculate DAJacCon::stateBoundaryCon_
Definition: DAJacCon.C:800
Foam::DAJacCon::dRdWTPreallocOff_
Vec dRdWTPreallocOff_
Definition: DAJacCon.H:91
Foam::DAJacCon::dRdWPreallocOn_
Vec dRdWPreallocOn_
Definition: DAJacCon.H:92
Foam::DAUtility::writeVectorASCII
static void writeVectorASCII(const Vec vecIn, const word prefix)
write petsc vector in ascii format
Definition: DAUtility.C:347
Foam::DAJacCon::nJacConColors_
label nJacConColors_
number of jacCon colors
Definition: DAJacCon.H:86
Foam::DAJacCon::neiBFaceGlobalCompact_
labelList neiBFaceGlobalCompact_
neibough face global index for a given local boundary face
Definition: DAJacCon.H:77
Foam::DAJacCon::createConnectionMat
void createConnectionMat(Mat *connectedStates)
allocate connectedState matrix
Definition: DAJacCon.C:146
Foam::DAJacCon::setupJacConPreallocation
void setupJacConPreallocation(const dictionary &options)
calculate the preallocation vector for initializing the JacCon mat, if necessary
Definition: DAJacCon.C:2016
Foam::DAJacCon::setConnections
void setConnections(Mat conMat, const label idx) const
add value 1 for the colume idx to conMat
Definition: DAJacCon.C:669
Foam::DAJacCon::calcJacConColoring
void calcJacConColoring(const word postFix="")
compute graph coloring for Jacobian connectivity matrix
Definition: DAJacCon.C:1917
Foam::DAJacCon::stateBoundaryConID_
Mat stateBoundaryConID_
matrix to store boundary connectivity ID for state Jacobians
Definition: DAJacCon.H:74
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
Foam::DAJacCon::preallocateJacobianMatrix
void preallocateJacobianMatrix(Mat dRMat, const Vec preallocOnProc, const Vec preallocOffProc) const
compute preallocation vectors
Definition: DAJacCon.C:173
Foam::DAJacCon::jacConColors_
Vec jacConColors_
jacCon matrix colors
Definition: DAJacCon.H:83
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:1637
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:1693
Foam::DAJacCon::setupdRdWCon
void setupdRdWCon(const HashTable< List< List< word >>> &stateResConInfo, const label isPrealloc)
Definition: DAJacCon.C:2039
Foam::DAJacCon::readJacConColoring
void readJacConColoring(const word postFix="")
read colors for JacCon
Definition: DAJacCon.C:1980
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:304
Foam::DAIndex::bFaceFaceI
labelList bFaceFaceI
Definition: DAIndex.H:116
Foam::DAJacCon::stateInfo_
HashTable< wordList > stateInfo_
the regState_ list from DAStateInfo object
Definition: DAJacCon.H:68
Foam::DAIndex::nLocalAdjointStates
label nLocalAdjointStates
number of local adjoint states (including all cells and faces)
Definition: DAIndex.H:145
Foam::DAJacCon::initializePetscVecs
void initializePetscVecs()
initialize petsc vectors
Definition: DAJacCon.C:83

Generated by   doxygen 1.8.17