Actual source code: pcpatch.c
1: #include <petsc/private/pcpatchimpl.h>
2: #include <petsc/private/kspimpl.h>
3: #include <petsc/private/vecimpl.h>
4: #include <petsc/private/dmpleximpl.h>
5: #include <petscsf.h>
6: #include <petscbt.h>
7: #include <petscds.h>
8: #include <../src/mat/impls/dense/seq/dense.h>
10: PetscBool PCPatchcite = PETSC_FALSE;
11: const char PCPatchCitation[] = "@article{FarrellKnepleyWechsungMitchell2020,\n"
12: "title = {{PCPATCH}: software for the topological construction of multigrid relaxation methods},\n"
13: "author = {Patrick E Farrell and Matthew G Knepley and Lawrence Mitchell and Florian Wechsung},\n"
14: "journal = {ACM Transaction on Mathematical Software},\n"
15: "eprint = {http://arxiv.org/abs/1912.08516},\n"
16: "volume = {47},\n"
17: "number = {3},\n"
18: "pages = {1--22},\n"
19: "year = {2021},\n"
20: "petsc_uses={KSP,DMPlex}\n}\n";
22: PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Apply, PC_Patch_Prealloc;
24: static inline PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
25: {
26: PetscCall(PetscViewerPushFormat(viewer, format));
27: PetscCall(PetscObjectView(obj, viewer));
28: PetscCall(PetscViewerPopFormat(viewer));
29: return PETSC_SUCCESS;
30: }
32: static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
33: {
34: PetscInt starSize;
35: PetscInt *star = NULL, si;
37: PetscFunctionBegin;
38: PetscCall(PetscHSetIClear(ht));
39: /* To start with, add the point we care about */
40: PetscCall(PetscHSetIAdd(ht, point));
41: /* Loop over all the points that this point connects to */
42: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
43: for (si = 0; si < starSize * 2; si += 2) PetscCall(PetscHSetIAdd(ht, star[si]));
44: PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
48: static PetscErrorCode PCPatchConstruct_Vanka(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
49: {
50: PC_PATCH *patch = (PC_PATCH *)vpatch;
51: PetscInt starSize;
52: PetscInt *star = NULL;
53: PetscBool shouldIgnore = PETSC_FALSE;
54: PetscInt cStart, cEnd, iStart, iEnd, si;
56: PetscFunctionBegin;
57: PetscCall(PetscHSetIClear(ht));
58: /* To start with, add the point we care about */
59: PetscCall(PetscHSetIAdd(ht, point));
60: /* Should we ignore any points of a certain dimension? */
61: if (patch->vankadim >= 0) {
62: shouldIgnore = PETSC_TRUE;
63: PetscCall(DMPlexGetDepthStratum(dm, patch->vankadim, &iStart, &iEnd));
64: }
65: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
66: /* Loop over all the cells that this point connects to */
67: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
68: for (si = 0; si < starSize * 2; si += 2) {
69: const PetscInt cell = star[si];
70: PetscInt closureSize;
71: PetscInt *closure = NULL, ci;
73: if (cell < cStart || cell >= cEnd) continue;
74: /* now loop over all entities in the closure of that cell */
75: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
76: for (ci = 0; ci < closureSize * 2; ci += 2) {
77: const PetscInt newpoint = closure[ci];
79: /* We've been told to ignore entities of this type.*/
80: if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue;
81: PetscCall(PetscHSetIAdd(ht, newpoint));
82: }
83: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
84: }
85: PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: static PetscErrorCode PCPatchConstruct_Pardecomp(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
90: {
91: PC_PATCH *patch = (PC_PATCH *)vpatch;
92: DMLabel ghost = NULL;
93: const PetscInt *leaves = NULL;
94: PetscInt nleaves = 0, pStart, pEnd, loc;
95: PetscBool isFiredrake;
96: PetscBool flg;
97: PetscInt starSize;
98: PetscInt *star = NULL;
99: PetscInt opoint, overlapi;
101: PetscFunctionBegin;
102: PetscCall(PetscHSetIClear(ht));
104: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
106: PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
107: if (isFiredrake) {
108: PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
109: PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
110: } else {
111: PetscSF sf;
112: PetscCall(DMGetPointSF(dm, &sf));
113: PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
114: nleaves = PetscMax(nleaves, 0);
115: }
117: for (opoint = pStart; opoint < pEnd; ++opoint) {
118: if (ghost) PetscCall(DMLabelHasPoint(ghost, opoint, &flg));
119: else {
120: PetscCall(PetscFindInt(opoint, nleaves, leaves, &loc));
121: flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
122: }
123: /* Not an owned entity, don't make a cell patch. */
124: if (flg) continue;
125: PetscCall(PetscHSetIAdd(ht, opoint));
126: }
128: /* Now build the overlap for the patch */
129: for (overlapi = 0; overlapi < patch->pardecomp_overlap; ++overlapi) {
130: PetscInt index = 0;
131: PetscInt *htpoints = NULL;
132: PetscInt htsize;
133: PetscInt i;
135: PetscCall(PetscHSetIGetSize(ht, &htsize));
136: PetscCall(PetscMalloc1(htsize, &htpoints));
137: PetscCall(PetscHSetIGetElems(ht, &index, htpoints));
139: for (i = 0; i < htsize; ++i) {
140: PetscInt hpoint = htpoints[i];
141: PetscInt si;
143: PetscCall(DMPlexGetTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star));
144: for (si = 0; si < starSize * 2; si += 2) {
145: const PetscInt starp = star[si];
146: PetscInt closureSize;
147: PetscInt *closure = NULL, ci;
149: /* now loop over all entities in the closure of starp */
150: PetscCall(DMPlexGetTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure));
151: for (ci = 0; ci < closureSize * 2; ci += 2) {
152: const PetscInt closstarp = closure[ci];
153: PetscCall(PetscHSetIAdd(ht, closstarp));
154: }
155: PetscCall(DMPlexRestoreTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure));
156: }
157: PetscCall(DMPlexRestoreTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star));
158: }
159: PetscCall(PetscFree(htpoints));
160: }
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /* The user's already set the patches in patch->userIS. Build the hash tables */
165: static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
166: {
167: PC_PATCH *patch = (PC_PATCH *)vpatch;
168: IS patchis = patch->userIS[point];
169: PetscInt n;
170: const PetscInt *patchdata;
171: PetscInt pStart, pEnd, i;
173: PetscFunctionBegin;
174: PetscCall(PetscHSetIClear(ht));
175: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
176: PetscCall(ISGetLocalSize(patchis, &n));
177: PetscCall(ISGetIndices(patchis, &patchdata));
178: for (i = 0; i < n; ++i) {
179: const PetscInt ownedpoint = patchdata[i];
181: PetscCheck(ownedpoint >= pStart && ownedpoint < pEnd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %" PetscInt_FMT " was not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", ownedpoint, pStart, pEnd);
182: PetscCall(PetscHSetIAdd(ht, ownedpoint));
183: }
184: PetscCall(ISRestoreIndices(patchis, &patchdata));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs)
189: {
190: PC_PATCH *patch = (PC_PATCH *)pc->data;
192: PetscFunctionBegin;
193: if (n == 1 && bs[0] == 1) {
194: patch->sectionSF = sf[0];
195: PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
196: } else {
197: PetscInt allRoots = 0, allLeaves = 0;
198: PetscInt leafOffset = 0;
199: PetscInt *ilocal = NULL;
200: PetscSFNode *iremote = NULL;
201: PetscInt *remoteOffsets = NULL;
202: PetscInt index = 0;
203: PetscHMapI rankToIndex;
204: PetscInt numRanks = 0;
205: PetscSFNode *remote = NULL;
206: PetscSF rankSF;
207: PetscInt *ranks = NULL;
208: PetscInt *offsets = NULL;
209: MPI_Datatype contig;
210: PetscHSetI ranksUniq;
211: PetscMPIInt in;
213: /* First figure out how many dofs there are in the concatenated numbering.
214: allRoots: number of owned global dofs;
215: allLeaves: number of visible dofs (global + ghosted).
216: */
217: for (PetscInt i = 0; i < n; ++i) {
218: PetscInt nroots, nleaves;
220: PetscCall(PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL));
221: allRoots += nroots * bs[i];
222: allLeaves += nleaves * bs[i];
223: }
224: PetscCall(PetscMalloc1(allLeaves, &ilocal));
225: PetscCall(PetscMalloc1(allLeaves, &iremote));
226: /* Now build an SF that just contains process connectivity. */
227: PetscCall(PetscHSetICreate(&ranksUniq));
228: for (PetscInt i = 0; i < n; ++i) {
229: const PetscMPIInt *ranks = NULL;
230: PetscMPIInt nranks;
232: PetscCall(PetscSFSetUp(sf[i]));
233: PetscCall(PetscSFGetRootRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL));
234: /* These are all the ranks who communicate with me. */
235: for (PetscMPIInt j = 0; j < nranks; ++j) PetscCall(PetscHSetIAdd(ranksUniq, (PetscInt)ranks[j]));
236: }
237: PetscCall(PetscHSetIGetSize(ranksUniq, &numRanks));
238: PetscCall(PetscMalloc1(numRanks, &remote));
239: PetscCall(PetscMalloc1(numRanks, &ranks));
240: PetscCall(PetscHSetIGetElems(ranksUniq, &index, ranks));
242: PetscCall(PetscHMapICreate(&rankToIndex));
243: for (PetscInt i = 0; i < numRanks; ++i) {
244: remote[i].rank = ranks[i];
245: remote[i].index = 0;
246: PetscCall(PetscHMapISet(rankToIndex, ranks[i], i));
247: }
248: PetscCall(PetscFree(ranks));
249: PetscCall(PetscHSetIDestroy(&ranksUniq));
250: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &rankSF));
251: PetscCall(PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
252: PetscCall(PetscSFSetUp(rankSF));
253: /* OK, use it to communicate the root offset on the remote processes for each subspace. */
254: PetscCall(PetscMalloc1(n, &offsets));
255: PetscCall(PetscMalloc1(n * numRanks, &remoteOffsets));
257: offsets[0] = 0;
258: for (PetscInt i = 1; i < n; ++i) {
259: PetscInt nroots;
261: PetscCall(PetscSFGetGraph(sf[i - 1], &nroots, NULL, NULL, NULL));
262: offsets[i] = offsets[i - 1] + nroots * bs[i - 1];
263: }
264: /* Offsets are the offsets on the current process of the global dof numbering for the subspaces. */
265: PetscCall(PetscMPIIntCast(n, &in));
266: PetscCallMPI(MPI_Type_contiguous(in, MPIU_INT, &contig));
267: PetscCallMPI(MPI_Type_commit(&contig));
269: PetscCall(PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets, MPI_REPLACE));
270: PetscCall(PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets, MPI_REPLACE));
271: PetscCallMPI(MPI_Type_free(&contig));
272: PetscCall(PetscFree(offsets));
273: PetscCall(PetscSFDestroy(&rankSF));
274: /* Now remoteOffsets contains the offsets on the remote
275: processes who communicate with me. So now we can
276: concatenate the list of SFs into a single one. */
277: index = 0;
278: for (PetscInt i = 0; i < n; ++i) {
279: const PetscSFNode *remote = NULL;
280: const PetscInt *local = NULL;
281: PetscInt nroots, nleaves, j;
283: PetscCall(PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote));
284: for (j = 0; j < nleaves; ++j) {
285: PetscInt rank = remote[j].rank;
286: PetscInt idx, rootOffset, k;
288: PetscCall(PetscHMapIGet(rankToIndex, rank, &idx));
289: PetscCheck(idx != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?");
290: /* Offset on given rank for ith subspace */
291: rootOffset = remoteOffsets[n * idx + i];
292: for (k = 0; k < bs[i]; ++k) {
293: ilocal[index] = (local ? local[j] : j) * bs[i] + k + leafOffset;
294: iremote[index].rank = remote[j].rank;
295: iremote[index].index = remote[j].index * bs[i] + k + rootOffset;
296: ++index;
297: }
298: }
299: leafOffset += nleaves * bs[i];
300: }
301: PetscCall(PetscHMapIDestroy(&rankToIndex));
302: PetscCall(PetscFree(remoteOffsets));
303: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->sectionSF));
304: PetscCall(PetscSFSetGraph(patch->sectionSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
305: }
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /* TODO: Docs */
310: static PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim)
311: {
312: PC_PATCH *patch = (PC_PATCH *)pc->data;
314: PetscFunctionBegin;
315: *dim = patch->ignoredim;
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@
320: PCPatchSetSaveOperators - Set whether the per-patch sub-matrices should be built and kept, instead of being reassembled at each application
322: Logically Collective
324: Input Parameters:
325: + pc - the `PCPATCH` preconditioner
326: - flg - `PETSC_TRUE` to store the assembled sub-matrices for each patch, `PETSC_FALSE` to rebuild them on demand
328: Level: intermediate
330: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchGetSaveOperators()`, `PCPatchSetPrecomputeElementTensors()`
331: @*/
332: PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg)
333: {
334: PC_PATCH *patch = (PC_PATCH *)pc->data;
336: PetscFunctionBegin;
337: patch->save_operators = flg;
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: PCPatchGetSaveOperators - Get whether the per-patch sub-matrices are built and kept between applications of the `PCPATCH` preconditioner
344: Not Collective
346: Input Parameter:
347: . pc - the `PCPATCH` preconditioner
349: Output Parameter:
350: . flg - `PETSC_TRUE` if the assembled sub-matrices are stored, `PETSC_FALSE` if they are rebuilt on demand
352: Level: intermediate
354: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchSetSaveOperators()`, `PCPatchSetPrecomputeElementTensors()`
355: @*/
356: PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg)
357: {
358: PC_PATCH *patch = (PC_PATCH *)pc->data;
360: PetscFunctionBegin;
361: *flg = patch->save_operators;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: PCPatchSetPrecomputeElementTensors - Set whether element tensors should be precomputed once and reused when assembling each patch matrix
368: Logically Collective
370: Input Parameters:
371: + pc - the `PCPATCH` preconditioner
372: - flg - `PETSC_TRUE` to precompute the element tensors, `PETSC_FALSE` to recompute them for each patch
374: Level: intermediate
376: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchGetPrecomputeElementTensors()`, `PCPatchSetSaveOperators()`
377: @*/
378: PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg)
379: {
380: PC_PATCH *patch = (PC_PATCH *)pc->data;
382: PetscFunctionBegin;
383: patch->precomputeElementTensors = flg;
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*@
388: PCPatchGetPrecomputeElementTensors - Get whether element tensors are precomputed once and reused when assembling each patch matrix
390: Not Collective
392: Input Parameter:
393: . pc - the `PCPATCH` preconditioner
395: Output Parameter:
396: . flg - `PETSC_TRUE` if the element tensors are precomputed, `PETSC_FALSE` if they are recomputed for each patch
398: Level: intermediate
400: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchSetPrecomputeElementTensors()`, `PCPatchSetSaveOperators()`
401: @*/
402: PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg)
403: {
404: PC_PATCH *patch = (PC_PATCH *)pc->data;
406: PetscFunctionBegin;
407: *flg = patch->precomputeElementTensors;
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: /*@
412: PCPatchSetPartitionOfUnity - Set whether the patch contributions should be weighted by a partition of unity when combining local solves
414: Logically Collective
416: Input Parameters:
417: + pc - the `PCPATCH` preconditioner
418: - flg - `PETSC_TRUE` to weight local patch updates by a partition of unity, `PETSC_FALSE` to sum them directly
420: Level: intermediate
422: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchGetPartitionOfUnity()`
423: @*/
424: PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg)
425: {
426: PC_PATCH *patch = (PC_PATCH *)pc->data;
428: PetscFunctionBegin;
429: patch->partition_of_unity = flg;
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: /*@
434: PCPatchGetPartitionOfUnity - Get whether the patch contributions are weighted by a partition of unity when combining local solves
436: Not Collective
438: Input Parameter:
439: . pc - the `PCPATCH` preconditioner
441: Output Parameter:
442: . flg - `PETSC_TRUE` if local patch updates are weighted by a partition of unity, `PETSC_FALSE` if they are summed directly
444: Level: intermediate
446: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchSetPartitionOfUnity()`
447: @*/
448: PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg)
449: {
450: PC_PATCH *patch = (PC_PATCH *)pc->data;
452: PetscFunctionBegin;
453: *flg = patch->partition_of_unity;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /* TODO: Docs */
458: static PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
459: {
460: PC_PATCH *patch = (PC_PATCH *)pc->data;
462: PetscFunctionBegin;
463: PetscCheck(type == PC_COMPOSITE_ADDITIVE || type == PC_COMPOSITE_MULTIPLICATIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Only supports additive or multiplicative as the local type");
464: patch->local_composition_type = type;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: PCPatchGetSubKSP - Get the per-patch `KSP` objects used to solve each local patch problem in a `PCPATCH` preconditioner
471: Not Collective
473: Input Parameter:
474: . pc - the `PCPATCH` preconditioner
476: Output Parameters:
477: + npatch - number of local patches (may be `NULL`)
478: - ksp - newly allocated array of length `npatch` holding the per-patch `KSP` objects; the caller must free the array with `PetscFree()`
480: Level: advanced
482: Note:
483: `PCSetUp()` must have been called on the `PCPATCH` (typically through `KSPSetUp()` on the outer `KSP`) before calling this routine.
485: .seealso: [](ch_ksp), `PCPATCH`, `KSP`, `PCASMGetSubKSP()`
486: @*/
487: PetscErrorCode PCPatchGetSubKSP(PC pc, PetscInt *npatch, KSP *ksp[])
488: {
489: PC_PATCH *patch = (PC_PATCH *)pc->data;
491: PetscFunctionBegin;
492: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here");
493: PetscCall(PetscMalloc1(patch->npatch, ksp));
494: for (PetscInt i = 0; i < patch->npatch; ++i) (*ksp)[i] = (KSP)patch->solver[i];
495: if (npatch) *npatch = patch->npatch;
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /*@
500: PCPatchSetSubMatType - Set the `MatType` used to store the per-patch sub-matrices in a `PCPATCH` preconditioner
502: Logically Collective
504: Input Parameters:
505: + pc - the `PCPATCH` preconditioner
506: - sub_mat_type - the `MatType` to use for the per-patch sub-matrices (e.g. `MATDENSE`, `MATSEQAIJ`)
508: Level: advanced
510: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchGetSubMatType()`, `MatType`
511: @*/
512: PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type)
513: {
514: PC_PATCH *patch = (PC_PATCH *)pc->data;
516: PetscFunctionBegin;
517: PetscCall(PetscFree(patch->sub_mat_type));
518: PetscCall(PetscStrallocpy(sub_mat_type, (char **)&patch->sub_mat_type));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*@
523: PCPatchGetSubMatType - Get the `MatType` used to store the per-patch sub-matrices in a `PCPATCH` preconditioner
525: Not Collective
527: Input Parameter:
528: . pc - the `PCPATCH` preconditioner
530: Output Parameter:
531: . sub_mat_type - the `MatType` used for the per-patch sub-matrices
533: Level: advanced
535: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchSetSubMatType()`, `MatType`
536: @*/
537: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
538: {
539: PC_PATCH *patch = (PC_PATCH *)pc->data;
541: PetscFunctionBegin;
542: *sub_mat_type = patch->sub_mat_type;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*@
547: PCPatchSetCellNumbering - Set the `PetscSection` that provides a numbering of the cells used to define patches in a `PCPATCH` preconditioner
549: Logically Collective
551: Input Parameters:
552: + pc - the `PCPATCH` preconditioner
553: - cellNumbering - the `PetscSection` giving the cell numbering; its reference count is incremented
555: Level: advanced
557: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchGetCellNumbering()`, `PetscSection`
558: @*/
559: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
560: {
561: PC_PATCH *patch = (PC_PATCH *)pc->data;
563: PetscFunctionBegin;
564: patch->cellNumbering = cellNumbering;
565: PetscCall(PetscObjectReference((PetscObject)cellNumbering));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: /*@
570: PCPatchGetCellNumbering - Get the `PetscSection` that provides the numbering of the cells used to define patches in a `PCPATCH` preconditioner
572: Not Collective
574: Input Parameter:
575: . pc - the `PCPATCH` preconditioner
577: Output Parameter:
578: . cellNumbering - the `PetscSection` giving the cell numbering
580: Level: advanced
582: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchSetCellNumbering()`, `PetscSection`
583: @*/
584: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
585: {
586: PC_PATCH *patch = (PC_PATCH *)pc->data;
588: PetscFunctionBegin;
589: *cellNumbering = patch->cellNumbering;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*@C
594: PCPatchSetConstructType - Set the way patches are constructed for a `PCPATCH` preconditioner
596: Logically Collective
598: Input Parameters:
599: + pc - the `PCPATCH` preconditioner
600: . ctype - the `PCPatchConstructType` selecting the patch construction strategy (e.g. `PC_PATCH_STAR`, `PC_PATCH_VANKA`, `PC_PATCH_PARDECOMP`, `PC_PATCH_USER`, `PC_PATCH_PYTHON`)
601: . func - user callback that builds the patches, used only when `ctype` is `PC_PATCH_USER` or `PC_PATCH_PYTHON`; may be `NULL` otherwise
602: - ctx - optional application context passed to `func`
604: Calling sequence of `func`:
605: + pc - the `PCPATCH` preconditioner
606: . npatch - number of patches
607: . patches - the `IS` that define each patch
608: . patchIterationSet - how the patches are iterated over
609: - ctx - optional application context
611: Level: advanced
613: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchGetConstructType()`, `PCPatchConstructType`
614: @*/
615: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC pc, PetscInt *npatch, IS *patches[], IS *patchIterationSet, PetscCtx ctx), PetscCtx ctx)
616: {
617: PC_PATCH *patch = (PC_PATCH *)pc->data;
619: PetscFunctionBegin;
620: patch->ctype = ctype;
621: switch (ctype) {
622: case PC_PATCH_STAR:
623: patch->user_patches = PETSC_FALSE;
624: patch->patchconstructop = PCPatchConstruct_Star;
625: break;
626: case PC_PATCH_VANKA:
627: patch->user_patches = PETSC_FALSE;
628: patch->patchconstructop = PCPatchConstruct_Vanka;
629: break;
630: case PC_PATCH_PARDECOMP:
631: patch->user_patches = PETSC_FALSE;
632: patch->patchconstructop = PCPatchConstruct_Pardecomp;
633: break;
634: case PC_PATCH_USER:
635: case PC_PATCH_PYTHON:
636: patch->user_patches = PETSC_TRUE;
637: patch->patchconstructop = PCPatchConstruct_User;
638: if (func) {
639: patch->userpatchconstructionop = func;
640: patch->userpatchconstructctx = ctx;
641: }
642: break;
643: default:
644: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
645: }
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*@C
650: PCPatchGetConstructType - Get the strategy currently used to construct patches for a `PCPATCH` preconditioner
652: Not Collective
654: Input Parameter:
655: . pc - the `PCPATCH` preconditioner
657: Output Parameters:
658: + ctype - the `PCPatchConstructType`
659: . func - the callback that builds the patches when `ctype` is `PC_PATCH_USER` or `PC_PATCH_PYTHON`; otherwise unchanged
660: - ctx - the application context associated with `func`; otherwise unchanged
662: Calling sequence of `func`:
663: + pc - the `PCPATCH` preconditioner
664: . npatch - number of patches
665: . patches - the `IS` that define each patch
666: . patchIterationSet - how the patches are iterated over
667: - ctx - optional application context
669: Level: advanced
671: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchSetConstructType()`, `PCPatchConstructType`
672: @*/
673: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC pc, PetscInt *npatch, IS *patches[], IS *patchIterationSet, PetscCtx ctx), PetscCtxRt ctx)
674: {
675: PC_PATCH *patch = (PC_PATCH *)pc->data;
677: PetscFunctionBegin;
678: *ctype = patch->ctype;
679: switch (patch->ctype) {
680: case PC_PATCH_STAR:
681: case PC_PATCH_VANKA:
682: case PC_PATCH_PARDECOMP:
683: break;
684: case PC_PATCH_USER:
685: case PC_PATCH_PYTHON:
686: *func = patch->userpatchconstructionop;
687: *(void **)ctx = patch->userpatchconstructctx;
688: break;
689: default:
690: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
691: }
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@C
696: PCPatchSetDiscretisationInfo - Provide the per-subspace discretisation information required by a `PCPATCH` preconditioner to build patch problems
698: Logically Collective
700: Input Parameters:
701: + pc - the `PCPATCH` preconditioner
702: . nsubspaces - the number of discretisation subspaces (e.g. fields)
703: . dms - array of length `nsubspaces` of `DM`s, one per subspace, from which the local sections and section `PetscSF`s are obtained
704: . bs - array of length `nsubspaces` giving the block size of each subspace
705: . nodesPerCell - array of length `nsubspaces` giving the number of nodes per cell for each subspace
706: . cellNodeMap - array of length `nsubspaces`; entry `i` is a cell-to-node map (array) of length `(cEnd - cStart) * nodesPerCell[i]`
707: . subspaceOffsets - array of length `nsubspaces + 1` giving the starting global dof offset of each subspace
708: . numGhostBcs - number of ghost (off-process) boundary-condition dofs
709: . ghostBcNodes - array of length `numGhostBcs` of the ghost boundary-condition dof indices
710: . numGlobalBcs - number of global boundary-condition dofs
711: - globalBcNodes - array of length `numGlobalBcs` of the global boundary-condition dof indices
713: Level: advanced
715: .seealso: [](ch_ksp), `PCPATCH`, `PCPatchSetComputeOperator()`, `PCPatchSetComputeFunction()`
716: @*/
717: PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM dms[], PetscInt bs[], PetscInt nodesPerCell[], const PetscInt **cellNodeMap, const PetscInt subspaceOffsets[], PetscInt numGhostBcs, const PetscInt ghostBcNodes[], PetscInt numGlobalBcs, const PetscInt globalBcNodes[])
718: {
719: PC_PATCH *patch = (PC_PATCH *)pc->data;
720: DM dm, plex;
721: PetscSF *sfs;
722: PetscInt cStart, cEnd, i, j;
724: PetscFunctionBegin;
725: PetscCall(PCGetDM(pc, &dm));
726: PetscCall(DMConvert(dm, DMPLEX, &plex));
727: dm = plex;
728: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
729: PetscCall(PetscMalloc1(nsubspaces, &sfs));
730: PetscCall(PetscMalloc1(nsubspaces, &patch->dofSection));
731: PetscCall(PetscMalloc1(nsubspaces, &patch->bs));
732: PetscCall(PetscMalloc1(nsubspaces, &patch->nodesPerCell));
733: PetscCall(PetscMalloc1(nsubspaces, &patch->cellNodeMap));
734: PetscCall(PetscMalloc1(nsubspaces + 1, &patch->subspaceOffsets));
736: patch->nsubspaces = nsubspaces;
737: patch->totalDofsPerCell = 0;
738: for (i = 0; i < nsubspaces; ++i) {
739: PetscCall(DMGetLocalSection(dms[i], &patch->dofSection[i]));
740: PetscCall(PetscObjectReference((PetscObject)patch->dofSection[i]));
741: PetscCall(DMGetSectionSF(dms[i], &sfs[i]));
742: patch->bs[i] = bs[i];
743: patch->nodesPerCell[i] = nodesPerCell[i];
744: patch->totalDofsPerCell += nodesPerCell[i] * bs[i];
745: PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
746: for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
747: patch->subspaceOffsets[i] = subspaceOffsets[i];
748: }
749: PetscCall(PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs));
750: PetscCall(PetscFree(sfs));
752: patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
753: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
754: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
755: PetscCall(DMDestroy(&dm));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: /* TODO: Docs */
760: static PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
761: {
762: PC_PATCH *patch = (PC_PATCH *)pc->data;
763: PetscInt cStart, cEnd, i, j;
765: PetscFunctionBegin;
766: patch->combined = PETSC_TRUE;
767: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
768: PetscCall(DMGetNumFields(dm, &patch->nsubspaces));
769: PetscCall(PetscCalloc1(patch->nsubspaces, &patch->dofSection));
770: PetscCall(PetscMalloc1(patch->nsubspaces, &patch->bs));
771: PetscCall(PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell));
772: PetscCall(PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap));
773: PetscCall(PetscCalloc1(patch->nsubspaces + 1, &patch->subspaceOffsets));
774: PetscCall(DMGetLocalSection(dm, &patch->dofSection[0]));
775: PetscCall(PetscObjectReference((PetscObject)patch->dofSection[0]));
776: PetscCall(PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]));
777: patch->totalDofsPerCell = 0;
778: for (i = 0; i < patch->nsubspaces; ++i) {
779: patch->bs[i] = 1;
780: patch->nodesPerCell[i] = nodesPerCell[i];
781: patch->totalDofsPerCell += nodesPerCell[i];
782: PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
783: for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
784: }
785: PetscCall(DMGetSectionSF(dm, &patch->sectionSF));
786: PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
787: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
788: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
789: PetscFunctionReturn(PETSC_SUCCESS);
790: }
792: /*@C
793: PCPatchSetComputeFunction - Set the callback function used to compute patch residuals
795: Logically Collective
797: Input Parameters:
798: + pc - The `PC`
799: . func - The callback function
800: - ctx - The application context
802: Calling sequence of `func`:
803: + pc - The `PC`
804: . point - The point
805: . x - The input solution (not used in linear problems)
806: . f - The patch residual vector
807: . cellIS - An array of the cell numbers
808: . n - The size of `dofsArray`
809: . dofsArray - The dofmap for the dofs to be solved for
810: . dofsArrayWithAll - The dofmap for all dofs on the patch
811: - ctx - The application context
813: Level: advanced
815: Note:
816: The entries of `f` (the output residual vector) have been set to zero before the call.
818: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunctionInteriorFacets()`
819: @*/
820: PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Vec f, IS cellIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, PetscCtx ctx), PetscCtx ctx)
821: {
822: PC_PATCH *patch = (PC_PATCH *)pc->data;
824: PetscFunctionBegin;
825: patch->usercomputef = func;
826: patch->usercomputefctx = ctx;
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: /*@C
831: PCPatchSetComputeFunctionInteriorFacets - Set the callback function used to compute facet integrals for patch residuals
833: Logically Collective
835: Input Parameters:
836: + pc - The `PC`
837: . func - The callback function
838: - ctx - The application context
840: Calling sequence of `func`:
841: + pc - The `PC`
842: . point - The point
843: . x - The input solution (not used in linear problems)
844: . f - The patch residual vector
845: . facetIS - An array of the facet numbers
846: . n - The size of `dofsArray`
847: . dofsArray - The dofmap for the dofs to be solved for
848: . dofsArrayWithAll - The dofmap for all dofs on the patch
849: - ctx - The application context
851: Level: advanced
853: Note:
854: The entries of `f` (the output residual vector) have been set to zero before the call.
856: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunction()`
857: @*/
858: PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Vec f, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, PetscCtx ctx), PetscCtx ctx)
859: {
860: PC_PATCH *patch = (PC_PATCH *)pc->data;
862: PetscFunctionBegin;
863: patch->usercomputefintfacet = func;
864: patch->usercomputefintfacetctx = ctx;
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*@C
869: PCPatchSetComputeOperator - Set the callback function used to compute patch matrices
871: Logically Collective
873: Input Parameters:
874: + pc - The `PC`
875: . func - The callback function
876: - ctx - The application context
878: Calling sequence of `func`:
879: + pc - The `PC`
880: . point - The point
881: . x - The input solution (not used in linear problems)
882: . mat - The patch matrix
883: . facetIS - An array of the cell numbers
884: . n - The size of `dofsArray`
885: . dofsArray - The dofmap for the dofs to be solved for
886: . dofsArrayWithAll - The dofmap for all dofs on the patch
887: - ctx - The application context
889: Level: advanced
891: Note:
892: The matrix entries have been set to zero before the call.
894: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
895: @*/
896: PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, PetscCtx ctx), PetscCtx ctx)
897: {
898: PC_PATCH *patch = (PC_PATCH *)pc->data;
900: PetscFunctionBegin;
901: patch->usercomputeop = func;
902: patch->usercomputeopctx = ctx;
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: /*@C
907: PCPatchSetComputeOperatorInteriorFacets - Set the callback function used to compute facet integrals for patch matrices
909: Logically Collective
911: Input Parameters:
912: + pc - The `PC`
913: . func - The callback function
914: - ctx - The application context
916: Calling sequence of `func`:
917: + pc - The `PC`
918: . point - The point
919: . x - The input solution (not used in linear problems)
920: . mat - The patch matrix
921: . facetIS - An array of the facet numbers
922: . n - The size of `dofsArray`
923: . dofsArray - The dofmap for the dofs to be solved for
924: . dofsArrayWithAll - The dofmap for all dofs on the patch
925: - ctx - The application context
927: Level: advanced
929: Note:
930: The matrix entries have been set to zero before the call.
932: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
933: @*/
934: PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, PetscCtx ctx), PetscCtx ctx)
935: {
936: PC_PATCH *patch = (PC_PATCH *)pc->data;
938: PetscFunctionBegin;
939: patch->usercomputeopintfacet = func;
940: patch->usercomputeopintfacetctx = ctx;
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: /*@C
945: PCPatchSetComputeOperatorExteriorFacets - Set the callback function used to compute exterior facet integrals for patch matrices
947: Logically Collective
949: Input Parameters:
950: + pc - The `PC`
951: . func - The callback function
952: - ctx - The application context
954: Calling sequence of `func`:
955: + pc - The `PC`
956: . point - The point
957: . x - The input solution (not used in linear problems)
958: . mat - The patch matrix
959: . facetIS - An array of the facet numbers
960: . n - The size of `dofsArray`
961: . dofsArray - The dofmap for the dofs to be solved for
962: . dofsArrayWithAll - The dofmap for all dofs on the patch
963: - ctx - The application context
965: Level: advanced
967: Note:
968: The matrix entries have been set to zero before the call.
970: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchSetComputeOperatorInteriorFacets()`, `PCPatchSetComputeFunctionExteriorFacets()`, `PCPatchSetDiscretisationInfo()`
971: @*/
972: PetscErrorCode PCPatchSetComputeOperatorExteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, PetscCtx ctx), PetscCtx ctx)
973: {
974: PC_PATCH *patch = (PC_PATCH *)pc->data;
976: PetscFunctionBegin;
977: patch->usercomputeopextfacet = func;
978: patch->usercomputeopextfacetctx = ctx;
979: PetscFunctionReturn(PETSC_SUCCESS);
980: }
982: /*@C
983: PCPatchSetComputeFunctionExteriorFacets - Set the callback function used to compute exterior facet integrals for patch residuals
985: Logically Collective
987: Input Parameters:
988: + pc - The `PC`
989: . func - The callback function
990: - ctx - The application context
992: Calling sequence of `func`:
993: + pc - The `PC`
994: . point - The point
995: . x - The input solution (not used in linear problems)
996: . f - The patch residual vector
997: . facetIS - An array of the facet numbers
998: . n - The size of `dofsArray`
999: . dofsArray - The dofmap for the dofs to be solved for
1000: . dofsArrayWithAll - The dofmap for all dofs on the patch
1001: - ctx - The application context
1003: Level: advanced
1005: Note:
1006: The entries of `f` (the output residual vector) have been set to zero before the call.
1008: .seealso: [](ch_ksp), `PCPatchSetComputeFunction()`, `PCPatchSetComputeFunctionInteriorFacets()`, `PCPatchSetComputeOperatorExteriorFacets()`, `PCPatchSetDiscretisationInfo()`
1009: @*/
1010: PetscErrorCode PCPatchSetComputeFunctionExteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Vec f, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, PetscCtx ctx), PetscCtx ctx)
1011: {
1012: PC_PATCH *patch = (PC_PATCH *)pc->data;
1014: PetscFunctionBegin;
1015: patch->usercomputefextfacet = func;
1016: patch->usercomputefextfacetctx = ctx;
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1020: /* On entry, ht contains the topological entities whose dofs we are responsible for solving for;
1021: on exit, cht contains all the topological entities we need to compute their residuals.
1022: In full generality this should incorporate knowledge of the sparsity pattern of the matrix;
1023: here we assume a standard FE sparsity pattern.*/
1024: /* TODO: Use DMPlexGetAdjacency() */
1025: static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht)
1026: {
1027: DM dm, plex;
1028: PC_PATCH *patch = (PC_PATCH *)pc->data;
1029: PetscHashIter hi;
1030: PetscInt point;
1031: PetscInt *star = NULL, *closure = NULL;
1032: PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci;
1033: PetscInt *fStar = NULL, *fClosure = NULL;
1034: PetscInt fBegin, fEnd, fsi, fci, fStarSize, fClosureSize;
1036: PetscFunctionBegin;
1037: PetscCall(PCGetDM(pc, &dm));
1038: PetscCall(DMConvert(dm, DMPLEX, &plex));
1039: dm = plex;
1040: PetscCall(DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd));
1041: PetscCall(PCPatchGetIgnoreDim(pc, &ignoredim));
1042: if (ignoredim >= 0) PetscCall(DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd));
1043: PetscCall(PetscHSetIClear(cht));
1044: PetscHashIterBegin(ht, hi);
1045: while (!PetscHashIterAtEnd(ht, hi)) {
1046: PetscHashIterGetKey(ht, hi, point);
1047: PetscHashIterNext(ht, hi);
1049: /* Loop over all the cells that this point connects to */
1050: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
1051: for (si = 0; si < starSize * 2; si += 2) {
1052: const PetscInt ownedpoint = star[si];
1053: /* TODO Check for point in cht before running through closure again */
1054: /* now loop over all entities in the closure of that cell */
1055: PetscCall(DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure));
1056: for (ci = 0; ci < closureSize * 2; ci += 2) {
1057: const PetscInt seenpoint = closure[ci];
1058: if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue;
1059: PetscCall(PetscHSetIAdd(cht, seenpoint));
1060: /* Facet integrals couple dofs across facets, so in that case for each of
1061: the facets we need to add all dofs on the other side of the facet to
1062: the seen dofs. */
1063: if (patch->usercomputeopintfacet) {
1064: if (fBegin <= seenpoint && seenpoint < fEnd) {
1065: PetscCall(DMPlexGetTransitiveClosure(dm, seenpoint, PETSC_FALSE, &fStarSize, &fStar));
1066: for (fsi = 0; fsi < fStarSize * 2; fsi += 2) {
1067: PetscCall(DMPlexGetTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, &fClosureSize, &fClosure));
1068: for (fci = 0; fci < fClosureSize * 2; fci += 2) PetscCall(PetscHSetIAdd(cht, fClosure[fci]));
1069: PetscCall(DMPlexRestoreTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, NULL, &fClosure));
1070: }
1071: PetscCall(DMPlexRestoreTransitiveClosure(dm, seenpoint, PETSC_FALSE, NULL, &fStar));
1072: }
1073: }
1074: }
1075: PetscCall(DMPlexRestoreTransitiveClosure(dm, ownedpoint, PETSC_TRUE, NULL, &closure));
1076: }
1077: PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, NULL, &star));
1078: }
1079: PetscCall(DMDestroy(&dm));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
1084: {
1085: PetscFunctionBegin;
1086: if (combined) {
1087: if (f < 0) {
1088: if (dof) PetscCall(PetscSectionGetDof(dofSection[0], p, dof));
1089: if (off) PetscCall(PetscSectionGetOffset(dofSection[0], p, off));
1090: } else {
1091: if (dof) PetscCall(PetscSectionGetFieldDof(dofSection[0], p, f, dof));
1092: if (off) PetscCall(PetscSectionGetFieldOffset(dofSection[0], p, f, off));
1093: }
1094: } else {
1095: if (f < 0) {
1096: PC_PATCH *patch = (PC_PATCH *)pc->data;
1097: PetscInt fdof, g;
1099: if (dof) {
1100: *dof = 0;
1101: for (g = 0; g < patch->nsubspaces; ++g) {
1102: PetscCall(PetscSectionGetDof(dofSection[g], p, &fdof));
1103: *dof += fdof;
1104: }
1105: }
1106: if (off) {
1107: *off = 0;
1108: for (g = 0; g < patch->nsubspaces; ++g) {
1109: PetscCall(PetscSectionGetOffset(dofSection[g], p, &fdof));
1110: *off += fdof;
1111: }
1112: }
1113: } else {
1114: if (dof) PetscCall(PetscSectionGetDof(dofSection[f], p, dof));
1115: if (off) PetscCall(PetscSectionGetOffset(dofSection[f], p, off));
1116: }
1117: }
1118: PetscFunctionReturn(PETSC_SUCCESS);
1119: }
1121: /* Given a hash table with a set of topological entities (pts), compute the degrees of
1122: freedom in global concatenated numbering on those entities.
1123: For Vanka smoothing, this needs to do something special: ignore dofs of the
1124: constraint subspace on entities that aren't the base entity we're building the patch
1125: around. */
1126: static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI *subspaces_to_exclude)
1127: {
1128: PC_PATCH *patch = (PC_PATCH *)pc->data;
1129: PetscHashIter hi;
1130: PetscInt ldof, loff;
1131: PetscInt k, p;
1133: PetscFunctionBegin;
1134: PetscCall(PetscHSetIClear(dofs));
1135: for (k = 0; k < patch->nsubspaces; ++k) {
1136: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1137: PetscInt bs = patch->bs[k];
1138: PetscInt j, l;
1140: if (subspaces_to_exclude != NULL) {
1141: PetscBool should_exclude_k = PETSC_FALSE;
1142: PetscCall(PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k));
1143: if (should_exclude_k) {
1144: /* only get this subspace dofs at the base entity, not any others */
1145: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff));
1146: if (0 == ldof) continue;
1147: for (j = loff; j < ldof + loff; ++j) {
1148: for (l = 0; l < bs; ++l) {
1149: PetscInt dof = bs * j + l + subspaceOffset;
1150: PetscCall(PetscHSetIAdd(dofs, dof));
1151: }
1152: }
1153: continue; /* skip the other dofs of this subspace */
1154: }
1155: }
1157: PetscHashIterBegin(pts, hi);
1158: while (!PetscHashIterAtEnd(pts, hi)) {
1159: PetscHashIterGetKey(pts, hi, p);
1160: PetscHashIterNext(pts, hi);
1161: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff));
1162: if (0 == ldof) continue;
1163: for (j = loff; j < ldof + loff; ++j) {
1164: for (l = 0; l < bs; ++l) {
1165: PetscInt dof = bs * j + l + subspaceOffset;
1166: PetscCall(PetscHSetIAdd(dofs, dof));
1167: }
1168: }
1169: }
1170: }
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
1175: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
1176: {
1177: PetscHashIter hi;
1178: PetscInt key;
1179: PetscBool flg;
1181: PetscFunctionBegin;
1182: PetscCall(PetscHSetIClear(C));
1183: PetscHashIterBegin(B, hi);
1184: while (!PetscHashIterAtEnd(B, hi)) {
1185: PetscHashIterGetKey(B, hi, key);
1186: PetscHashIterNext(B, hi);
1187: PetscCall(PetscHSetIHas(A, key, &flg));
1188: if (!flg) PetscCall(PetscHSetIAdd(C, key));
1189: }
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: // PetscClangLinter pragma disable: -fdoc-sowing-chars
1194: /*
1195: PCPatchCreateCellPatches - create patches.
1197: Input Parameter:
1198: . dm - The DMPlex object defining the mesh
1200: Output Parameters:
1201: + cellCounts - Section with counts of cells around each vertex
1202: . cells - IS of the cell point indices of cells in each patch
1203: . pointCounts - Section with counts of cells around each vertex
1204: - point - IS of the cell point indices of cells in each patch
1205: */
1206: static PetscErrorCode PCPatchCreateCellPatches(PC pc)
1207: {
1208: PC_PATCH *patch = (PC_PATCH *)pc->data;
1209: DMLabel ghost = NULL;
1210: DM dm, plex;
1211: PetscHSetI ht = NULL, cht = NULL;
1212: PetscSection cellCounts, pointCounts, intFacetCounts, extFacetCounts;
1213: PetscInt *cellsArray, *pointsArray, *intFacetsArray, *extFacetsArray, *intFacetsToPatchCell, *extFacetsToPatchCell;
1214: PetscInt numCells, numPoints, numIntFacets, numExtFacets;
1215: const PetscInt *leaves;
1216: PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, v;
1217: PetscBool isFiredrake;
1219: PetscFunctionBegin;
1220: /* Used to keep track of the cells in the patch. */
1221: PetscCall(PetscHSetICreate(&ht));
1222: PetscCall(PetscHSetICreate(&cht));
1224: PetscCall(PCGetDM(pc, &dm));
1225: PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC");
1226: PetscCall(DMConvert(dm, DMPLEX, &plex));
1227: dm = plex;
1228: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1229: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1231: if (patch->user_patches) {
1232: PetscCall(patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx));
1233: vStart = 0;
1234: vEnd = patch->npatch;
1235: } else if (patch->ctype == PC_PATCH_PARDECOMP) {
1236: vStart = 0;
1237: vEnd = 1;
1238: } else if (patch->codim < 0) {
1239: if (patch->dim < 0) PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1240: else PetscCall(DMPlexGetDepthStratum(dm, patch->dim, &vStart, &vEnd));
1241: } else PetscCall(DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd));
1242: patch->npatch = vEnd - vStart;
1244: /* These labels mark the owned points. We only create patches around points that this process owns. */
1245: PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
1246: if (isFiredrake) {
1247: PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
1248: PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
1249: } else {
1250: PetscSF sf;
1252: PetscCall(DMGetPointSF(dm, &sf));
1253: PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
1254: nleaves = PetscMax(nleaves, 0);
1255: }
1257: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts));
1258: PetscCall(PetscObjectSetName((PetscObject)patch->cellCounts, "Patch Cell Layout"));
1259: cellCounts = patch->cellCounts;
1260: PetscCall(PetscSectionSetChart(cellCounts, vStart, vEnd));
1261: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts));
1262: PetscCall(PetscObjectSetName((PetscObject)patch->pointCounts, "Patch Point Layout"));
1263: pointCounts = patch->pointCounts;
1264: PetscCall(PetscSectionSetChart(pointCounts, vStart, vEnd));
1265: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts));
1266: PetscCall(PetscObjectSetName((PetscObject)patch->extFacetCounts, "Patch Exterior Facet Layout"));
1267: extFacetCounts = patch->extFacetCounts;
1268: PetscCall(PetscSectionSetChart(extFacetCounts, vStart, vEnd));
1269: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts));
1270: PetscCall(PetscObjectSetName((PetscObject)patch->intFacetCounts, "Patch Interior Facet Layout"));
1271: intFacetCounts = patch->intFacetCounts;
1272: PetscCall(PetscSectionSetChart(intFacetCounts, vStart, vEnd));
1273: /* Count cells and points in the patch surrounding each entity */
1274: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
1275: for (v = vStart; v < vEnd; ++v) {
1276: PetscHashIter hi;
1277: PetscInt chtSize, loc = -1;
1278: PetscBool flg;
1280: if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
1281: if (ghost) PetscCall(DMLabelHasPoint(ghost, v, &flg));
1282: else {
1283: PetscCall(PetscFindInt(v, nleaves, leaves, &loc));
1284: flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
1285: }
1286: /* Not an owned entity, don't make a cell patch. */
1287: if (flg) continue;
1288: }
1290: PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
1291: PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
1292: PetscCall(PetscHSetIGetSize(cht, &chtSize));
1293: /* empty patch, continue */
1294: if (chtSize == 0) continue;
1296: /* safe because size(cht) > 0 from above */
1297: PetscHashIterBegin(cht, hi);
1298: while (!PetscHashIterAtEnd(cht, hi)) {
1299: PetscInt point, pdof;
1301: PetscHashIterGetKey(cht, hi, point);
1302: if (fStart <= point && point < fEnd) {
1303: const PetscInt *support;
1304: PetscInt supportSize, p;
1305: PetscBool interior = PETSC_TRUE;
1306: PetscCall(DMPlexGetSupport(dm, point, &support));
1307: PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
1308: if (supportSize == 1) {
1309: interior = PETSC_FALSE;
1310: } else {
1311: for (p = 0; p < supportSize; p++) {
1312: PetscBool found;
1313: /* FIXME: can I do this while iterating over cht? */
1314: PetscCall(PetscHSetIHas(cht, support[p], &found));
1315: if (!found) {
1316: interior = PETSC_FALSE;
1317: break;
1318: }
1319: }
1320: }
1321: if (interior) {
1322: PetscCall(PetscSectionAddDof(intFacetCounts, v, 1));
1323: } else {
1324: PetscCall(PetscSectionAddDof(extFacetCounts, v, 1));
1325: }
1326: }
1327: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
1328: if (pdof) PetscCall(PetscSectionAddDof(pointCounts, v, 1));
1329: if (point >= cStart && point < cEnd) PetscCall(PetscSectionAddDof(cellCounts, v, 1));
1330: PetscHashIterNext(cht, hi);
1331: }
1332: }
1333: if (isFiredrake) PetscCall(DMLabelDestroyIndex(ghost));
1335: PetscCall(PetscSectionSetUp(cellCounts));
1336: PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1337: PetscCall(PetscMalloc1(numCells, &cellsArray));
1338: PetscCall(PetscSectionSetUp(pointCounts));
1339: PetscCall(PetscSectionGetStorageSize(pointCounts, &numPoints));
1340: PetscCall(PetscMalloc1(numPoints, &pointsArray));
1342: PetscCall(PetscSectionSetUp(intFacetCounts));
1343: PetscCall(PetscSectionSetUp(extFacetCounts));
1344: PetscCall(PetscSectionGetStorageSize(intFacetCounts, &numIntFacets));
1345: PetscCall(PetscSectionGetStorageSize(extFacetCounts, &numExtFacets));
1346: PetscCall(PetscMalloc1(numIntFacets, &intFacetsArray));
1347: PetscCall(PetscMalloc1(numIntFacets * 2, &intFacetsToPatchCell));
1348: PetscCall(PetscMalloc1(numExtFacets, &extFacetsArray));
1349: PetscCall(PetscMalloc1(numExtFacets, &extFacetsToPatchCell));
1351: /* Now that we know how much space we need, run through again and actually remember the cells. */
1352: for (v = vStart; v < vEnd; v++) {
1353: PetscHashIter hi;
1354: PetscInt dof, off, cdof, coff, efdof, efoff, ifdof, ifoff, pdof, n = 0, cn = 0, ifn = 0, efn = 0;
1356: PetscCall(PetscSectionGetDof(pointCounts, v, &dof));
1357: PetscCall(PetscSectionGetOffset(pointCounts, v, &off));
1358: PetscCall(PetscSectionGetDof(cellCounts, v, &cdof));
1359: PetscCall(PetscSectionGetOffset(cellCounts, v, &coff));
1360: PetscCall(PetscSectionGetDof(intFacetCounts, v, &ifdof));
1361: PetscCall(PetscSectionGetOffset(intFacetCounts, v, &ifoff));
1362: PetscCall(PetscSectionGetDof(extFacetCounts, v, &efdof));
1363: PetscCall(PetscSectionGetOffset(extFacetCounts, v, &efoff));
1364: if (dof <= 0) continue;
1365: PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
1366: PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
1367: PetscHashIterBegin(cht, hi);
1368: while (!PetscHashIterAtEnd(cht, hi)) {
1369: PetscInt point;
1371: PetscHashIterGetKey(cht, hi, point);
1372: if (fStart <= point && point < fEnd) {
1373: const PetscInt *support;
1374: PetscInt supportSize, p;
1375: PetscBool interior = PETSC_TRUE;
1376: PetscCall(DMPlexGetSupport(dm, point, &support));
1377: PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
1378: if (supportSize == 1) {
1379: interior = PETSC_FALSE;
1380: } else {
1381: for (p = 0; p < supportSize; p++) {
1382: PetscBool found;
1383: /* FIXME: can I do this while iterating over cht? */
1384: PetscCall(PetscHSetIHas(cht, support[p], &found));
1385: if (!found) {
1386: interior = PETSC_FALSE;
1387: break;
1388: }
1389: }
1390: }
1391: if (interior) {
1392: intFacetsToPatchCell[2 * (ifoff + ifn)] = support[0];
1393: intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = support[1];
1394: intFacetsArray[ifoff + ifn++] = point;
1395: } else {
1396: /* Find the support cell that is in the patch */
1397: PetscInt supportCell = -1;
1398: for (p = 0; p < supportSize; p++) {
1399: PetscBool found;
1400: PetscCall(PetscHSetIHas(cht, support[p], &found));
1401: if (found && support[p] >= cStart && support[p] < cEnd) {
1402: supportCell = support[p];
1403: break;
1404: }
1405: }
1406: extFacetsToPatchCell[efoff + efn] = supportCell;
1407: extFacetsArray[efoff + efn++] = point;
1408: }
1409: }
1410: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
1411: if (pdof) pointsArray[off + n++] = point;
1412: if (point >= cStart && point < cEnd) cellsArray[coff + cn++] = point;
1413: PetscHashIterNext(cht, hi);
1414: }
1415: PetscCheck(ifn == ifdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of interior facets in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, ifn, ifdof);
1416: PetscCheck(efn == efdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of exterior facets in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, efn, efdof);
1417: PetscCheck(cn == cdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of cells in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, cn, cdof);
1418: PetscCheck(n == dof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of points in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, n, dof);
1420: for (ifn = 0; ifn < ifdof; ifn++) {
1421: PetscInt cell0 = intFacetsToPatchCell[2 * (ifoff + ifn)];
1422: PetscInt cell1 = intFacetsToPatchCell[2 * (ifoff + ifn) + 1];
1423: PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
1424: for (n = 0; n < cdof; n++) {
1425: if (!found0 && cell0 == cellsArray[coff + n]) {
1426: intFacetsToPatchCell[2 * (ifoff + ifn)] = n;
1427: found0 = PETSC_TRUE;
1428: }
1429: if (!found1 && cell1 == cellsArray[coff + n]) {
1430: intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = n;
1431: found1 = PETSC_TRUE;
1432: }
1433: if (found0 && found1) break;
1434: }
1435: PetscCheck(found0 && found1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support");
1436: }
1437: for (efn = 0; efn < efdof; efn++) {
1438: PetscInt cell0 = extFacetsToPatchCell[efoff + efn];
1439: PetscBool found0 = PETSC_FALSE;
1440: for (n = 0; n < cdof; n++) {
1441: if (cell0 == cellsArray[coff + n]) {
1442: extFacetsToPatchCell[efoff + efn] = n;
1443: found0 = PETSC_TRUE;
1444: break;
1445: }
1446: }
1447: PetscCheck(found0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point number for exterior facet support");
1448: }
1449: }
1450: PetscCall(PetscHSetIDestroy(&ht));
1451: PetscCall(PetscHSetIDestroy(&cht));
1453: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells));
1454: PetscCall(PetscObjectSetName((PetscObject)patch->cells, "Patch Cells"));
1455: if (patch->viewCells) {
1456: PetscCall(ObjectView((PetscObject)patch->cellCounts, patch->viewerCells, patch->formatCells));
1457: PetscCall(ObjectView((PetscObject)patch->cells, patch->viewerCells, patch->formatCells));
1458: }
1459: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets));
1460: PetscCall(PetscObjectSetName((PetscObject)patch->intFacets, "Patch Interior Facets"));
1461: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2 * numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell));
1462: PetscCall(PetscObjectSetName((PetscObject)patch->intFacetsToPatchCell, "Patch Interior Facets local support"));
1463: if (patch->viewIntFacets) {
1464: PetscCall(ObjectView((PetscObject)patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets));
1465: PetscCall(ObjectView((PetscObject)patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets));
1466: PetscCall(ObjectView((PetscObject)patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets));
1467: }
1468: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets));
1469: PetscCall(PetscObjectSetName((PetscObject)patch->extFacets, "Patch Exterior Facets"));
1470: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsToPatchCell, PETSC_OWN_POINTER, &patch->extFacetsToPatchCell));
1471: PetscCall(PetscObjectSetName((PetscObject)patch->extFacetsToPatchCell, "Patch Exterior Facets local support"));
1472: if (patch->viewExtFacets) {
1473: PetscCall(ObjectView((PetscObject)patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets));
1474: PetscCall(ObjectView((PetscObject)patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets));
1475: }
1476: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points));
1477: PetscCall(PetscObjectSetName((PetscObject)patch->points, "Patch Points"));
1478: if (patch->viewPoints) {
1479: PetscCall(ObjectView((PetscObject)patch->pointCounts, patch->viewerPoints, patch->formatPoints));
1480: PetscCall(ObjectView((PetscObject)patch->points, patch->viewerPoints, patch->formatPoints));
1481: }
1482: PetscCall(DMDestroy(&dm));
1483: PetscFunctionReturn(PETSC_SUCCESS);
1484: }
1486: /*
1487: PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches
1489: Input Parameters:
1490: + dm - The DMPlex object defining the mesh
1491: . cellCounts - Section with counts of cells around each vertex
1492: . cells - IS of the cell point indices of cells in each patch
1493: . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1494: . nodesPerCell - number of nodes per cell.
1495: - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)
1497: Output Parameters:
1498: + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering
1499: . gtolCounts - Section with counts of dofs per cell patch
1500: - gtol - IS mapping from global dofs to local dofs for each patch.
1501: */
1502: static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc)
1503: {
1504: PC_PATCH *patch = (PC_PATCH *)pc->data;
1505: PetscSection cellCounts = patch->cellCounts;
1506: PetscSection pointCounts = patch->pointCounts;
1507: PetscSection gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL;
1508: IS cells = patch->cells;
1509: IS points = patch->points;
1510: PetscSection cellNumbering = patch->cellNumbering;
1511: PetscInt Nf = patch->nsubspaces;
1512: PetscInt numCells, numPoints;
1513: PetscInt numDofs;
1514: PetscInt numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll;
1515: PetscInt totalDofsPerCell = patch->totalDofsPerCell;
1516: PetscInt vStart, vEnd, v;
1517: const PetscInt *cellsArray, *pointsArray;
1518: PetscInt *newCellsArray = NULL;
1519: PetscInt *dofsArray = NULL;
1520: PetscInt *dofsArrayWithArtificial = NULL;
1521: PetscInt *dofsArrayWithAll = NULL;
1522: PetscInt *offsArray = NULL;
1523: PetscInt *offsArrayWithArtificial = NULL;
1524: PetscInt *offsArrayWithAll = NULL;
1525: PetscInt *asmArray = NULL;
1526: PetscInt *asmArrayWithArtificial = NULL;
1527: PetscInt *asmArrayWithAll = NULL;
1528: PetscInt *globalDofsArray = NULL;
1529: PetscInt *globalDofsArrayWithArtificial = NULL;
1530: PetscInt *globalDofsArrayWithAll = NULL;
1531: PetscInt globalIndex = 0;
1532: PetscInt key = 0;
1533: PetscInt asmKey = 0;
1534: DM dm = NULL, plex;
1535: const PetscInt *bcNodes = NULL;
1536: PetscHMapI ht;
1537: PetscHMapI htWithArtificial;
1538: PetscHMapI htWithAll;
1539: PetscHSetI globalBcs;
1540: PetscInt numBcs;
1541: PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs;
1542: PetscInt pStart, pEnd, p, i;
1543: char option[PETSC_MAX_PATH_LEN];
1544: PetscBool isNonlinear;
1546: PetscFunctionBegin;
1547: PetscCall(PCGetDM(pc, &dm));
1548: PetscCall(DMConvert(dm, DMPLEX, &plex));
1549: dm = plex;
1550: /* dofcounts section is cellcounts section * dofPerCell */
1551: PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1552: PetscCall(PetscSectionGetStorageSize(patch->pointCounts, &numPoints));
1553: numDofs = numCells * totalDofsPerCell;
1554: PetscCall(PetscMalloc1(numDofs, &dofsArray));
1555: PetscCall(PetscMalloc1(numPoints * Nf, &offsArray));
1556: PetscCall(PetscMalloc1(numDofs, &asmArray));
1557: PetscCall(PetscMalloc1(numCells, &newCellsArray));
1558: PetscCall(PetscSectionGetChart(cellCounts, &vStart, &vEnd));
1559: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts));
1560: gtolCounts = patch->gtolCounts;
1561: PetscCall(PetscSectionSetChart(gtolCounts, vStart, vEnd));
1562: PetscCall(PetscObjectSetName((PetscObject)patch->gtolCounts, "Patch Global Index Section"));
1564: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1565: PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithArtificial));
1566: PetscCall(PetscMalloc1(numDofs, &asmArrayWithArtificial));
1567: PetscCall(PetscMalloc1(numDofs, &dofsArrayWithArtificial));
1568: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial));
1569: gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1570: PetscCall(PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd));
1571: PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs"));
1572: }
1574: isNonlinear = patch->isNonlinear;
1575: if (isNonlinear) {
1576: PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithAll));
1577: PetscCall(PetscMalloc1(numDofs, &asmArrayWithAll));
1578: PetscCall(PetscMalloc1(numDofs, &dofsArrayWithAll));
1579: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll));
1580: gtolCountsWithAll = patch->gtolCountsWithAll;
1581: PetscCall(PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd));
1582: PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs"));
1583: }
1585: /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1586: conditions */
1587: PetscCall(PetscHSetICreate(&globalBcs));
1588: PetscCall(ISGetIndices(patch->ghostBcNodes, &bcNodes));
1589: PetscCall(ISGetSize(patch->ghostBcNodes, &numBcs));
1590: for (i = 0; i < numBcs; ++i) PetscCall(PetscHSetIAdd(globalBcs, bcNodes[i])); /* these are already in concatenated numbering */
1591: PetscCall(ISRestoreIndices(patch->ghostBcNodes, &bcNodes));
1592: PetscCall(ISDestroy(&patch->ghostBcNodes)); /* memory optimisation */
1594: /* Hash tables for artificial BC construction */
1595: PetscCall(PetscHSetICreate(&ownedpts));
1596: PetscCall(PetscHSetICreate(&seenpts));
1597: PetscCall(PetscHSetICreate(&owneddofs));
1598: PetscCall(PetscHSetICreate(&seendofs));
1599: PetscCall(PetscHSetICreate(&artificialbcs));
1601: PetscCall(ISGetIndices(cells, &cellsArray));
1602: PetscCall(ISGetIndices(points, &pointsArray));
1603: PetscCall(PetscHMapICreate(&ht));
1604: PetscCall(PetscHMapICreate(&htWithArtificial));
1605: PetscCall(PetscHMapICreate(&htWithAll));
1606: for (v = vStart; v < vEnd; ++v) {
1607: PetscInt localIndex = 0;
1608: PetscInt localIndexWithArtificial = 0;
1609: PetscInt localIndexWithAll = 0;
1610: PetscInt dof, off, i, j, k, l;
1612: PetscCall(PetscHMapIClear(ht));
1613: PetscCall(PetscHMapIClear(htWithArtificial));
1614: PetscCall(PetscHMapIClear(htWithAll));
1615: PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1616: PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1617: if (dof <= 0) continue;
1619: /* Calculate the global numbers of the artificial BC dofs here first */
1620: PetscCall(patch->patchconstructop((void *)patch, dm, v, ownedpts));
1621: PetscCall(PCPatchCompleteCellPatch(pc, ownedpts, seenpts));
1622: PetscCall(PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude));
1623: PetscCall(PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL));
1624: PetscCall(PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs));
1625: if (patch->viewPatches) {
1626: PetscHSetI globalbcdofs;
1627: PetscHashIter hi;
1628: MPI_Comm comm = PetscObjectComm((PetscObject)pc);
1630: PetscCall(PetscHSetICreate(&globalbcdofs));
1631: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": owned dofs:\n", v));
1632: PetscHashIterBegin(owneddofs, hi);
1633: while (!PetscHashIterAtEnd(owneddofs, hi)) {
1634: PetscInt globalDof;
1636: PetscHashIterGetKey(owneddofs, hi, globalDof);
1637: PetscHashIterNext(owneddofs, hi);
1638: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1639: }
1640: PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1641: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": seen dofs:\n", v));
1642: PetscHashIterBegin(seendofs, hi);
1643: while (!PetscHashIterAtEnd(seendofs, hi)) {
1644: PetscInt globalDof;
1645: PetscBool flg;
1647: PetscHashIterGetKey(seendofs, hi, globalDof);
1648: PetscHashIterNext(seendofs, hi);
1649: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1651: PetscCall(PetscHSetIHas(globalBcs, globalDof, &flg));
1652: if (flg) PetscCall(PetscHSetIAdd(globalbcdofs, globalDof));
1653: }
1654: PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1655: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": global BCs:\n", v));
1656: PetscCall(PetscHSetIGetSize(globalbcdofs, &numBcs));
1657: if (numBcs > 0) {
1658: PetscHashIterBegin(globalbcdofs, hi);
1659: while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1660: PetscInt globalDof;
1661: PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1662: PetscHashIterNext(globalbcdofs, hi);
1663: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1664: }
1665: }
1666: PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1667: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": artificial BCs:\n", v));
1668: PetscCall(PetscHSetIGetSize(artificialbcs, &numBcs));
1669: if (numBcs > 0) {
1670: PetscHashIterBegin(artificialbcs, hi);
1671: while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1672: PetscInt globalDof;
1673: PetscHashIterGetKey(artificialbcs, hi, globalDof);
1674: PetscHashIterNext(artificialbcs, hi);
1675: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1676: }
1677: }
1678: PetscCall(PetscSynchronizedPrintf(comm, "\n\n"));
1679: PetscCall(PetscHSetIDestroy(&globalbcdofs));
1680: }
1681: for (k = 0; k < patch->nsubspaces; ++k) {
1682: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1683: PetscInt nodesPerCell = patch->nodesPerCell[k];
1684: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1685: PetscInt bs = patch->bs[k];
1687: for (i = off; i < off + dof; ++i) {
1688: /* Walk over the cells in this patch. */
1689: const PetscInt c = cellsArray[i];
1690: PetscInt cell = c;
1692: /* TODO Change this to an IS */
1693: if (cellNumbering) {
1694: PetscCall(PetscSectionGetDof(cellNumbering, c, &cell));
1695: PetscCheck(cell > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " doesn't appear in cell numbering map", c);
1696: PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1697: }
1698: newCellsArray[i] = cell;
1699: for (j = 0; j < nodesPerCell; ++j) {
1700: /* For each global dof, map it into contiguous local storage. */
1701: const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + subspaceOffset;
1702: /* finally, loop over block size */
1703: for (l = 0; l < bs; ++l) {
1704: PetscInt localDof;
1705: PetscBool isGlobalBcDof, isArtificialBcDof;
1707: /* first, check if this is either a globally enforced or locally enforced BC dof */
1708: PetscCall(PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof));
1709: PetscCall(PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof));
1711: /* if it's either, don't ever give it a local dof number */
1712: if (isGlobalBcDof || isArtificialBcDof) {
1713: dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1714: } else {
1715: PetscCall(PetscHMapIGet(ht, globalDof + l, &localDof));
1716: if (localDof == -1) {
1717: localDof = localIndex++;
1718: PetscCall(PetscHMapISet(ht, globalDof + l, localDof));
1719: }
1720: PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1721: /* And store. */
1722: dofsArray[globalIndex] = localDof;
1723: }
1725: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1726: if (isGlobalBcDof) {
1727: dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1728: } else {
1729: PetscCall(PetscHMapIGet(htWithArtificial, globalDof + l, &localDof));
1730: if (localDof == -1) {
1731: localDof = localIndexWithArtificial++;
1732: PetscCall(PetscHMapISet(htWithArtificial, globalDof + l, localDof));
1733: }
1734: PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1735: /* And store.*/
1736: dofsArrayWithArtificial[globalIndex] = localDof;
1737: }
1738: }
1740: if (isNonlinear) {
1741: /* Build the dofmap for the function space with _all_ dofs,
1742: including those in any kind of boundary condition */
1743: PetscCall(PetscHMapIGet(htWithAll, globalDof + l, &localDof));
1744: if (localDof == -1) {
1745: localDof = localIndexWithAll++;
1746: PetscCall(PetscHMapISet(htWithAll, globalDof + l, localDof));
1747: }
1748: PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1749: /* And store.*/
1750: dofsArrayWithAll[globalIndex] = localDof;
1751: }
1752: globalIndex++;
1753: }
1754: }
1755: }
1756: }
1757: /* How many local dofs in this patch? */
1758: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1759: PetscCall(PetscHMapIGetSize(htWithArtificial, &dof));
1760: PetscCall(PetscSectionSetDof(gtolCountsWithArtificial, v, dof));
1761: }
1762: if (isNonlinear) {
1763: PetscCall(PetscHMapIGetSize(htWithAll, &dof));
1764: PetscCall(PetscSectionSetDof(gtolCountsWithAll, v, dof));
1765: }
1766: PetscCall(PetscHMapIGetSize(ht, &dof));
1767: PetscCall(PetscSectionSetDof(gtolCounts, v, dof));
1768: }
1770: PetscCall(DMDestroy(&dm));
1771: PetscCheck(globalIndex == numDofs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%" PetscInt_FMT ") doesn't match found number (%" PetscInt_FMT ")", numDofs, globalIndex);
1772: PetscCall(PetscSectionSetUp(gtolCounts));
1773: PetscCall(PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs));
1774: PetscCall(PetscMalloc1(numGlobalDofs, &globalDofsArray));
1776: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1777: PetscCall(PetscSectionSetUp(gtolCountsWithArtificial));
1778: PetscCall(PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial));
1779: PetscCall(PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial));
1780: }
1781: if (isNonlinear) {
1782: PetscCall(PetscSectionSetUp(gtolCountsWithAll));
1783: PetscCall(PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll));
1784: PetscCall(PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll));
1785: }
1786: /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */
1787: for (v = vStart; v < vEnd; ++v) {
1788: PetscHashIter hi;
1789: PetscInt dof, off, Np, ooff, i, j, k, l;
1791: PetscCall(PetscHMapIClear(ht));
1792: PetscCall(PetscHMapIClear(htWithArtificial));
1793: PetscCall(PetscHMapIClear(htWithAll));
1794: PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1795: PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1796: PetscCall(PetscSectionGetDof(pointCounts, v, &Np));
1797: PetscCall(PetscSectionGetOffset(pointCounts, v, &ooff));
1798: if (dof <= 0) continue;
1800: for (k = 0; k < patch->nsubspaces; ++k) {
1801: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1802: PetscInt nodesPerCell = patch->nodesPerCell[k];
1803: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1804: PetscInt bs = patch->bs[k];
1805: PetscInt goff;
1807: for (i = off; i < off + dof; ++i) {
1808: /* Reconstruct mapping of global-to-local on this patch. */
1809: const PetscInt c = cellsArray[i];
1810: PetscInt cell = c;
1812: if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1813: for (j = 0; j < nodesPerCell; ++j) {
1814: for (l = 0; l < bs; ++l) {
1815: const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1816: const PetscInt localDof = dofsArray[key];
1817: if (localDof >= 0) PetscCall(PetscHMapISet(ht, globalDof, localDof));
1818: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1819: const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1820: if (localDofWithArtificial >= 0) PetscCall(PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial));
1821: }
1822: if (isNonlinear) {
1823: const PetscInt localDofWithAll = dofsArrayWithAll[key];
1824: if (localDofWithAll >= 0) PetscCall(PetscHMapISet(htWithAll, globalDof, localDofWithAll));
1825: }
1826: key++;
1827: }
1828: }
1829: }
1831: /* Shove it in the output data structure. */
1832: PetscCall(PetscSectionGetOffset(gtolCounts, v, &goff));
1833: PetscHashIterBegin(ht, hi);
1834: while (!PetscHashIterAtEnd(ht, hi)) {
1835: PetscInt globalDof, localDof;
1837: PetscHashIterGetKey(ht, hi, globalDof);
1838: PetscHashIterGetVal(ht, hi, localDof);
1839: if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1840: PetscHashIterNext(ht, hi);
1841: }
1843: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1844: PetscCall(PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff));
1845: PetscHashIterBegin(htWithArtificial, hi);
1846: while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1847: PetscInt globalDof, localDof;
1848: PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1849: PetscHashIterGetVal(htWithArtificial, hi, localDof);
1850: if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1851: PetscHashIterNext(htWithArtificial, hi);
1852: }
1853: }
1854: if (isNonlinear) {
1855: PetscCall(PetscSectionGetOffset(gtolCountsWithAll, v, &goff));
1856: PetscHashIterBegin(htWithAll, hi);
1857: while (!PetscHashIterAtEnd(htWithAll, hi)) {
1858: PetscInt globalDof, localDof;
1859: PetscHashIterGetKey(htWithAll, hi, globalDof);
1860: PetscHashIterGetVal(htWithAll, hi, localDof);
1861: if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1862: PetscHashIterNext(htWithAll, hi);
1863: }
1864: }
1866: for (p = 0; p < Np; ++p) {
1867: const PetscInt point = pointsArray[ooff + p];
1868: PetscInt globalDof, localDof;
1870: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof));
1871: PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
1872: offsArray[(ooff + p) * Nf + k] = localDof;
1873: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1874: PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
1875: offsArrayWithArtificial[(ooff + p) * Nf + k] = localDof;
1876: }
1877: if (isNonlinear) {
1878: PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
1879: offsArrayWithAll[(ooff + p) * Nf + k] = localDof;
1880: }
1881: }
1882: }
1884: PetscCall(PetscHSetIDestroy(&globalBcs));
1885: PetscCall(PetscHSetIDestroy(&ownedpts));
1886: PetscCall(PetscHSetIDestroy(&seenpts));
1887: PetscCall(PetscHSetIDestroy(&owneddofs));
1888: PetscCall(PetscHSetIDestroy(&seendofs));
1889: PetscCall(PetscHSetIDestroy(&artificialbcs));
1891: /* At this point, we have a hash table ht built that maps globalDof -> localDof.
1892: We need to create the dof table laid out cellwise first, then by subspace,
1893: as the assembler assembles cell-wise and we need to stuff the different
1894: contributions of the different function spaces to the right places. So we loop
1895: over cells, then over subspaces. */
1896: if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */
1897: for (i = off; i < off + dof; ++i) {
1898: const PetscInt c = cellsArray[i];
1899: PetscInt cell = c;
1901: if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1902: for (k = 0; k < patch->nsubspaces; ++k) {
1903: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1904: PetscInt nodesPerCell = patch->nodesPerCell[k];
1905: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1906: PetscInt bs = patch->bs[k];
1908: for (j = 0; j < nodesPerCell; ++j) {
1909: for (l = 0; l < bs; ++l) {
1910: const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1911: PetscInt localDof;
1913: PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
1914: /* If it's not in the hash table, i.e. is a BC dof,
1915: then the PetscHSetIMap above gives -1, which matches
1916: exactly the convention for PETSc's matrix assembly to
1917: ignore the dof. So we don't need to do anything here */
1918: asmArray[asmKey] = localDof;
1919: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1920: PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
1921: asmArrayWithArtificial[asmKey] = localDof;
1922: }
1923: if (isNonlinear) {
1924: PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
1925: asmArrayWithAll[asmKey] = localDof;
1926: }
1927: asmKey++;
1928: }
1929: }
1930: }
1931: }
1932: }
1933: }
1934: if (1 == patch->nsubspaces) {
1935: PetscCall(PetscArraycpy(asmArray, dofsArray, numDofs));
1936: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs));
1937: if (isNonlinear) PetscCall(PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs));
1938: }
1940: PetscCall(PetscHMapIDestroy(&ht));
1941: PetscCall(PetscHMapIDestroy(&htWithArtificial));
1942: PetscCall(PetscHMapIDestroy(&htWithAll));
1943: PetscCall(ISRestoreIndices(cells, &cellsArray));
1944: PetscCall(ISRestoreIndices(points, &pointsArray));
1945: PetscCall(PetscFree(dofsArray));
1946: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscFree(dofsArrayWithArtificial));
1947: if (isNonlinear) PetscCall(PetscFree(dofsArrayWithAll));
1948: /* Create placeholder section for map from points to patch dofs */
1949: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection));
1950: PetscCall(PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces));
1951: if (patch->combined) {
1952: PetscInt numFields;
1953: PetscCall(PetscSectionGetNumFields(patch->dofSection[0], &numFields));
1954: PetscCheck(numFields == patch->nsubspaces, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Mismatch between number of section fields %" PetscInt_FMT " and number of subspaces %" PetscInt_FMT, numFields, patch->nsubspaces);
1955: PetscCall(PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd));
1956: PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
1957: for (p = pStart; p < pEnd; ++p) {
1958: PetscInt dof, fdof, f;
1960: PetscCall(PetscSectionGetDof(patch->dofSection[0], p, &dof));
1961: PetscCall(PetscSectionSetDof(patch->patchSection, p, dof));
1962: for (f = 0; f < patch->nsubspaces; ++f) {
1963: PetscCall(PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof));
1964: PetscCall(PetscSectionSetFieldDof(patch->patchSection, p, f, fdof));
1965: }
1966: }
1967: } else {
1968: PetscInt pStartf, pEndf, f;
1969: pStart = PETSC_INT_MAX;
1970: pEnd = PETSC_INT_MIN;
1971: for (f = 0; f < patch->nsubspaces; ++f) {
1972: PetscCall(PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf));
1973: pStart = PetscMin(pStart, pStartf);
1974: pEnd = PetscMax(pEnd, pEndf);
1975: }
1976: PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
1977: for (f = 0; f < patch->nsubspaces; ++f) {
1978: PetscCall(PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf));
1979: for (p = pStartf; p < pEndf; ++p) {
1980: PetscInt fdof;
1981: PetscCall(PetscSectionGetDof(patch->dofSection[f], p, &fdof));
1982: PetscCall(PetscSectionAddDof(patch->patchSection, p, fdof));
1983: PetscCall(PetscSectionSetFieldDof(patch->patchSection, p, f, fdof));
1984: }
1985: }
1986: }
1987: PetscCall(PetscSectionSetUp(patch->patchSection));
1988: PetscCall(PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE));
1989: /* Replace cell indices with firedrake-numbered ones. */
1990: PetscCall(ISGeneralSetIndices(cells, numCells, (const PetscInt *)newCellsArray, PETSC_OWN_POINTER));
1991: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol));
1992: PetscCall(PetscObjectSetName((PetscObject)patch->gtol, "Global Indices"));
1993: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname));
1994: PetscCall(PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject)pc, option));
1995: PetscCall(ISViewFromOptions(patch->gtol, (PetscObject)pc, option));
1996: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs));
1997: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArray, PETSC_OWN_POINTER, &patch->offs));
1998: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1999: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial));
2000: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial));
2001: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial));
2002: }
2003: if (isNonlinear) {
2004: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll));
2005: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll));
2006: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll));
2007: }
2008: PetscFunctionReturn(PETSC_SUCCESS);
2009: }
2011: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
2012: {
2013: PC_PATCH *patch = (PC_PATCH *)pc->data;
2014: PetscBool flg;
2015: PetscInt csize, rsize;
2016: const char *prefix = NULL;
2018: PetscFunctionBegin;
2019: if (withArtificial) {
2020: /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
2021: PetscInt pStart;
2022: PetscCall(PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL));
2023: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize));
2024: csize = rsize;
2025: } else {
2026: PetscInt pStart;
2027: PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
2028: PetscCall(PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize));
2029: csize = rsize;
2030: }
2032: PetscCall(MatCreate(PETSC_COMM_SELF, mat));
2033: PetscCall(PCGetOptionsPrefix(pc, &prefix));
2034: PetscCall(MatSetOptionsPrefix(*mat, prefix));
2035: PetscCall(MatAppendOptionsPrefix(*mat, "pc_patch_sub_"));
2036: if (patch->sub_mat_type) PetscCall(MatSetType(*mat, patch->sub_mat_type));
2037: else if (!patch->sub_mat_type) PetscCall(MatSetType(*mat, MATDENSE));
2038: PetscCall(MatSetSizes(*mat, rsize, csize, rsize, csize));
2039: PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATDENSE, &flg));
2040: if (!flg) PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg));
2041: /* Sparse patch matrices */
2042: if (!flg) {
2043: PetscBT bt;
2044: PetscInt *dnnz = NULL;
2045: const PetscInt *dofsArray = NULL;
2046: PetscInt pStart, pEnd, ncell, offset, c, i, j;
2048: if (withArtificial) {
2049: PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
2050: } else {
2051: PetscCall(ISGetIndices(patch->dofs, &dofsArray));
2052: }
2053: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
2054: point += pStart;
2055: PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
2056: PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
2057: PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
2058: PetscCall(PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0));
2059: /* A PetscBT uses N^2 bits to store the sparsity pattern on a
2060: * patch. This is probably OK if the patches are not too big,
2061: * but uses too much memory. We therefore switch based on rsize. */
2062: if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
2063: PetscScalar *zeroes;
2064: PetscInt rows;
2066: PetscCall(PetscCalloc1(rsize, &dnnz));
2067: PetscCall(PetscBTCreate(rsize * rsize, &bt));
2068: for (c = 0; c < ncell; ++c) {
2069: const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
2070: for (i = 0; i < patch->totalDofsPerCell; ++i) {
2071: const PetscInt row = idx[i];
2072: if (row < 0) continue;
2073: for (j = 0; j < patch->totalDofsPerCell; ++j) {
2074: const PetscInt col = idx[j];
2075: const PetscInt key = row * rsize + col;
2076: if (col < 0) continue;
2077: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
2078: }
2079: }
2080: }
2082: if (patch->usercomputeopintfacet) {
2083: const PetscInt *intFacetsArray = NULL;
2084: PetscInt i, numIntFacets, intFacetOffset;
2085: const PetscInt *facetCells = NULL;
2087: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
2088: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
2089: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
2090: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2091: for (i = 0; i < numIntFacets; i++) {
2092: const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
2093: const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
2094: PetscInt celli, cellj;
2096: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
2097: const PetscInt row = dofsArray[(offset + cell0) * patch->totalDofsPerCell + celli];
2098: if (row < 0) continue;
2099: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
2100: const PetscInt col = dofsArray[(offset + cell1) * patch->totalDofsPerCell + cellj];
2101: const PetscInt key = row * rsize + col;
2102: if (col < 0) continue;
2103: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
2104: }
2105: }
2107: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
2108: const PetscInt row = dofsArray[(offset + cell1) * patch->totalDofsPerCell + celli];
2109: if (row < 0) continue;
2110: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
2111: const PetscInt col = dofsArray[(offset + cell0) * patch->totalDofsPerCell + cellj];
2112: const PetscInt key = row * rsize + col;
2113: if (col < 0) continue;
2114: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
2115: }
2116: }
2117: }
2118: }
2119: PetscCall(PetscBTDestroy(&bt));
2120: PetscCall(MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL));
2121: PetscCall(PetscFree(dnnz));
2123: PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &zeroes));
2124: for (c = 0; c < ncell; ++c) {
2125: const PetscInt *idx = &dofsArray[(offset + c) * patch->totalDofsPerCell];
2126: PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES));
2127: }
2128: PetscCall(MatGetLocalSize(*mat, &rows, NULL));
2129: for (i = 0; i < rows; ++i) PetscCall(MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES));
2131: if (patch->usercomputeopintfacet) {
2132: const PetscInt *intFacetsArray = NULL;
2133: PetscInt i, numIntFacets, intFacetOffset;
2134: const PetscInt *facetCells = NULL;
2136: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
2137: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
2138: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
2139: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2140: for (i = 0; i < numIntFacets; i++) {
2141: const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
2142: const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
2143: const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
2144: const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
2145: PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES));
2146: PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES));
2147: }
2148: }
2150: /* Exterior facet preallocation: each exterior facet touches one cell */
2151: if (patch->usercomputeopextfacet) {
2152: PetscInt i, numExtFacets, extFacetOffset;
2153: const PetscInt *extFacetCells = NULL;
2155: PetscCall(PetscSectionGetDof(patch->extFacetCounts, point, &numExtFacets));
2156: PetscCall(PetscSectionGetOffset(patch->extFacetCounts, point, &extFacetOffset));
2157: PetscCall(ISGetIndices(patch->extFacetsToPatchCell, &extFacetCells));
2158: for (i = 0; i < numExtFacets; i++) {
2159: const PetscInt cell0 = extFacetCells[extFacetOffset + i];
2160: const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
2161: PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES));
2162: }
2163: PetscCall(ISRestoreIndices(patch->extFacetsToPatchCell, &extFacetCells));
2164: }
2166: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
2167: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
2169: PetscCall(PetscFree(zeroes));
2171: } else { /* rsize too big, use MATPREALLOCATOR */
2172: Mat preallocator;
2173: PetscScalar *vals;
2175: PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &vals));
2176: PetscCall(MatCreate(PETSC_COMM_SELF, &preallocator));
2177: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
2178: PetscCall(MatSetSizes(preallocator, rsize, rsize, rsize, rsize));
2179: PetscCall(MatSetUp(preallocator));
2181: for (c = 0; c < ncell; ++c) {
2182: const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
2183: PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES));
2184: }
2186: if (patch->usercomputeopintfacet) {
2187: const PetscInt *intFacetsArray = NULL;
2188: PetscInt i, numIntFacets, intFacetOffset;
2189: const PetscInt *facetCells = NULL;
2191: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
2192: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
2193: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
2194: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2195: for (i = 0; i < numIntFacets; i++) {
2196: const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
2197: const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
2198: const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
2199: const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
2200: PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES));
2201: PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES));
2202: }
2203: }
2205: /* Exterior facet preallocation: each exterior facet touches one cell */
2206: if (patch->usercomputeopextfacet) {
2207: PetscInt i, numExtFacets, extFacetOffset;
2208: const PetscInt *extFacetCells = NULL;
2210: PetscCall(PetscSectionGetDof(patch->extFacetCounts, point, &numExtFacets));
2211: PetscCall(PetscSectionGetOffset(patch->extFacetCounts, point, &extFacetOffset));
2212: PetscCall(ISGetIndices(patch->extFacetsToPatchCell, &extFacetCells));
2213: for (i = 0; i < numExtFacets; i++) {
2214: const PetscInt cell0 = extFacetCells[extFacetOffset + i];
2215: const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
2216: PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES));
2217: }
2218: PetscCall(ISRestoreIndices(patch->extFacetsToPatchCell, &extFacetCells));
2219: }
2221: PetscCall(PetscFree(vals));
2222: PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
2223: PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
2224: PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat));
2225: PetscCall(MatDestroy(&preallocator));
2226: }
2227: PetscCall(PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0));
2228: if (withArtificial) {
2229: PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
2230: } else {
2231: PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
2232: }
2233: }
2234: PetscCall(MatSetUp(*mat));
2235: PetscFunctionReturn(PETSC_SUCCESS);
2236: }
2238: static PetscErrorCode PCPatchComputeFunction_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Vec F, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, PetscCtx ctx)
2239: {
2240: PC_PATCH *patch = (PC_PATCH *)pc->data;
2241: DM dm, plex;
2242: PetscSection s;
2243: const PetscInt *parray, *oarray;
2244: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
2246: PetscFunctionBegin;
2247: PetscCheck(!patch->precomputeElementTensors, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute function");
2248: PetscCall(PCGetDM(pc, &dm));
2249: PetscCall(DMConvert(dm, DMPLEX, &plex));
2250: dm = plex;
2251: PetscCall(DMGetLocalSection(dm, &s));
2252: /* Set offset into patch */
2253: PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
2254: PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
2255: PetscCall(ISGetIndices(patch->points, &parray));
2256: PetscCall(ISGetIndices(patch->offs, &oarray));
2257: for (f = 0; f < Nf; ++f) {
2258: for (p = 0; p < Np; ++p) {
2259: const PetscInt point = parray[poff + p];
2260: PetscInt dof;
2262: PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
2263: PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
2264: if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
2265: else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
2266: }
2267: }
2268: PetscCall(ISRestoreIndices(patch->points, &parray));
2269: PetscCall(ISRestoreIndices(patch->offs, &oarray));
2270: if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
2271: PetscCall(DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx));
2272: PetscCall(DMDestroy(&dm));
2273: PetscFunctionReturn(PETSC_SUCCESS);
2274: }
2276: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
2277: {
2278: PC_PATCH *patch = (PC_PATCH *)pc->data;
2279: const PetscInt *dofsArray;
2280: const PetscInt *dofsArrayWithAll;
2281: const PetscInt *cellsArray;
2282: PetscInt ncell, offset, pStart, pEnd;
2284: PetscFunctionBegin;
2285: PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
2286: PetscCheck(patch->usercomputef || patch->usercomputefintfacet || patch->usercomputefextfacet, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeFunction(), PCPatchSetComputeFunctionInteriorFacets(), or PCPatchSetComputeFunctionExteriorFacets() to set callback");
2287: PetscCall(ISGetIndices(patch->dofs, &dofsArray));
2288: PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
2289: PetscCall(ISGetIndices(patch->cells, &cellsArray));
2290: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
2292: point += pStart;
2293: PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
2295: PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
2296: PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
2297: if (ncell <= 0) {
2298: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2299: PetscFunctionReturn(PETSC_SUCCESS);
2300: }
2301: PetscCall(VecSet(F, 0.0));
2302: if (patch->usercomputef) {
2303: /* Cannot reuse the same IS because the geometry info is being cached in it */
2304: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS));
2305: PetscCallBack("PCPatch callback", patch->usercomputef(pc, point, x, F, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, dofsArrayWithAll + offset * patch->totalDofsPerCell, patch->usercomputefctx));
2306: PetscCall(ISDestroy(&patch->cellIS));
2307: }
2308: if (patch->usercomputefextfacet) {
2309: PetscInt numExtFacets, extFacetOffset;
2310: PetscCall(PetscSectionGetDof(patch->extFacetCounts, point, &numExtFacets));
2311: PetscCall(PetscSectionGetOffset(patch->extFacetCounts, point, &extFacetOffset));
2312: if (numExtFacets > 0) {
2313: PetscInt *facetDofs = NULL;
2314: const PetscInt *extFacetsArray = NULL, *extFacetCells = NULL;
2315: PetscInt idx = 0;
2316: IS facetIS = NULL;
2318: PetscCall(ISGetIndices(patch->extFacetsToPatchCell, &extFacetCells));
2319: PetscCall(ISGetIndices(patch->extFacets, &extFacetsArray));
2320: PetscCall(PetscMalloc1(patch->totalDofsPerCell * numExtFacets, &facetDofs));
2321: for (PetscInt i = 0; i < numExtFacets; i++) {
2322: const PetscInt cell = extFacetCells[extFacetOffset + i];
2323: for (PetscInt d = 0; d < patch->totalDofsPerCell; d++) {
2324: facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2325: idx++;
2326: }
2327: }
2328: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray + extFacetOffset, PETSC_USE_POINTER, &facetIS));
2329: PetscCall(patch->usercomputefextfacet(pc, point, x, F, facetIS, numExtFacets * patch->totalDofsPerCell, facetDofs, dofsArrayWithAll + offset * patch->totalDofsPerCell, patch->usercomputefextfacetctx));
2330: PetscCall(ISDestroy(&facetIS));
2331: PetscCall(ISRestoreIndices(patch->extFacetsToPatchCell, &extFacetCells));
2332: PetscCall(ISRestoreIndices(patch->extFacets, &extFacetsArray));
2333: PetscCall(PetscFree(facetDofs));
2334: }
2335: }
2336: PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
2337: PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
2338: PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
2339: if (patch->viewMatrix) {
2340: char name[PETSC_MAX_PATH_LEN];
2342: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch vector for Point %" PetscInt_FMT, point));
2343: PetscCall(PetscObjectSetName((PetscObject)F, name));
2344: PetscCall(ObjectView((PetscObject)F, patch->viewerMatrix, patch->formatMatrix));
2345: }
2346: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2347: PetscFunctionReturn(PETSC_SUCCESS);
2348: }
2350: static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, PetscCtx ctx)
2351: {
2352: PC_PATCH *patch = (PC_PATCH *)pc->data;
2353: DM dm, plex;
2354: PetscSection s;
2355: const PetscInt *parray, *oarray;
2356: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
2358: PetscFunctionBegin;
2359: PetscCall(PCGetDM(pc, &dm));
2360: PetscCall(DMConvert(dm, DMPLEX, &plex));
2361: dm = plex;
2362: PetscCall(DMGetLocalSection(dm, &s));
2363: /* Set offset into patch */
2364: PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
2365: PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
2366: PetscCall(ISGetIndices(patch->points, &parray));
2367: PetscCall(ISGetIndices(patch->offs, &oarray));
2368: for (f = 0; f < Nf; ++f) {
2369: for (p = 0; p < Np; ++p) {
2370: const PetscInt point = parray[poff + p];
2371: PetscInt dof;
2373: PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
2374: PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
2375: if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
2376: else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
2377: }
2378: }
2379: PetscCall(ISRestoreIndices(patch->points, &parray));
2380: PetscCall(ISRestoreIndices(patch->offs, &oarray));
2381: if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
2382: /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
2383: PetscCall(DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx));
2384: PetscCall(DMDestroy(&dm));
2385: PetscFunctionReturn(PETSC_SUCCESS);
2386: }
2388: /* This function zeros mat on entry */
2389: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2390: {
2391: PC_PATCH *patch = (PC_PATCH *)pc->data;
2392: const PetscInt *dofsArray;
2393: const PetscInt *dofsArrayWithAll = NULL;
2394: const PetscInt *cellsArray;
2395: PetscInt ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2396: PetscBool isNonlinear;
2398: PetscFunctionBegin;
2399: PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
2400: isNonlinear = patch->isNonlinear;
2401: PetscCheck(patch->usercomputeop || patch->usercomputeopintfacet || patch->usercomputeopextfacet, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator(), PCPatchSetComputeOperatorInteriorFacets(), or PCPatchSetComputeOperatorExteriorFacets() to set callback");
2402: if (withArtificial) {
2403: PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
2404: } else {
2405: PetscCall(ISGetIndices(patch->dofs, &dofsArray));
2406: }
2407: if (isNonlinear) PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
2408: PetscCall(ISGetIndices(patch->cells, &cellsArray));
2409: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
2411: point += pStart;
2412: PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
2414: PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
2415: PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
2416: if (ncell <= 0) {
2417: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2418: PetscFunctionReturn(PETSC_SUCCESS);
2419: }
2420: PetscCall(MatZeroEntries(mat));
2421: if (patch->usercomputeop) {
2422: if (patch->precomputeElementTensors) {
2423: PetscInt i;
2424: PetscInt ndof = patch->totalDofsPerCell;
2425: const PetscScalar *elementTensors;
2427: PetscCall(VecGetArrayRead(patch->cellMats, &elementTensors));
2428: for (i = 0; i < ncell; i++) {
2429: const PetscInt cell = cellsArray[i + offset];
2430: const PetscInt *idx = dofsArray + (offset + i) * ndof;
2431: const PetscScalar *v = elementTensors + patch->precomputedTensorLocations[cell] * ndof * ndof;
2432: PetscCall(MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES));
2433: }
2434: PetscCall(VecRestoreArrayRead(patch->cellMats, &elementTensors));
2435: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2436: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2437: } else {
2438: /* Cannot reuse the same IS because the geometry info is being cached in it */
2439: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS));
2440: PetscCallBack("PCPatch callback",
2441: patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, PetscSafePointerPlusOffset(dofsArrayWithAll, offset * patch->totalDofsPerCell), patch->usercomputeopctx));
2442: }
2443: }
2444: if (patch->usercomputeopintfacet) {
2445: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
2446: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
2447: if (numIntFacets > 0) {
2448: /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2449: PetscInt *facetDofs = NULL, *facetDofsWithAll = NULL;
2450: const PetscInt *intFacetsArray = NULL;
2451: PetscInt idx = 0;
2452: PetscInt i, c, d;
2453: PetscInt fStart;
2454: DM dm, plex;
2455: IS facetIS = NULL;
2456: const PetscInt *facetCells = NULL;
2458: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
2459: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2460: PetscCall(PCGetDM(pc, &dm));
2461: PetscCall(DMConvert(dm, DMPLEX, &plex));
2462: dm = plex;
2463: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, NULL));
2464: /* FIXME: Pull this malloc out. */
2465: PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs));
2466: if (dofsArrayWithAll) PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll));
2467: if (patch->precomputeElementTensors) {
2468: PetscInt nFacetDof = 2 * patch->totalDofsPerCell;
2469: const PetscScalar *elementTensors;
2471: PetscCall(VecGetArrayRead(patch->intFacetMats, &elementTensors));
2473: for (i = 0; i < numIntFacets; i++) {
2474: const PetscInt facet = intFacetsArray[i + intFacetOffset];
2475: const PetscScalar *v = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart] * nFacetDof * nFacetDof;
2476: idx = 0;
2477: /*
2478: 0--1
2479: |\-|
2480: |+\|
2481: 2--3
2482: [0, 2, 3, 0, 1, 3]
2483: */
2484: for (c = 0; c < 2; c++) {
2485: const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
2486: for (d = 0; d < patch->totalDofsPerCell; d++) {
2487: facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2488: idx++;
2489: }
2490: }
2491: PetscCall(MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES));
2492: }
2493: PetscCall(VecRestoreArrayRead(patch->intFacetMats, &elementTensors));
2494: } else {
2495: /*
2496: 0--1
2497: |\-|
2498: |+\|
2499: 2--3
2500: [0, 2, 3, 0, 1, 3]
2501: */
2502: for (i = 0; i < numIntFacets; i++) {
2503: for (c = 0; c < 2; c++) {
2504: const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
2505: for (d = 0; d < patch->totalDofsPerCell; d++) {
2506: facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2507: if (dofsArrayWithAll) facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell) * patch->totalDofsPerCell + d];
2508: idx++;
2509: }
2510: }
2511: }
2512: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS));
2513: PetscCall(patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2 * numIntFacets * patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx));
2514: PetscCall(ISDestroy(&facetIS));
2515: }
2516: PetscCall(ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells));
2517: PetscCall(ISRestoreIndices(patch->intFacets, &intFacetsArray));
2518: PetscCall(PetscFree(facetDofs));
2519: PetscCall(PetscFree(facetDofsWithAll));
2520: PetscCall(DMDestroy(&dm));
2521: }
2522: }
2523: if (patch->usercomputeopextfacet) {
2524: PetscInt numExtFacets, extFacetOffset;
2525: PetscCall(PetscSectionGetDof(patch->extFacetCounts, point, &numExtFacets));
2526: PetscCall(PetscSectionGetOffset(patch->extFacetCounts, point, &extFacetOffset));
2527: if (numExtFacets > 0) {
2528: /* For each exterior facet, grab the one cell (in local numbering, and build dof numbering for that cell) */
2529: PetscInt *facetDofs = NULL, *facetDofsWithAll = NULL;
2530: const PetscInt *extFacetsArray = NULL, *extFacetCells = NULL;
2531: PetscInt idx = 0;
2532: IS facetIS = NULL;
2534: PetscCall(ISGetIndices(patch->extFacetsToPatchCell, &extFacetCells));
2535: PetscCall(ISGetIndices(patch->extFacets, &extFacetsArray));
2536: /* FIXME: Pull this malloc out. */
2537: PetscCall(PetscMalloc1(patch->totalDofsPerCell * numExtFacets, &facetDofs));
2538: if (dofsArrayWithAll) PetscCall(PetscMalloc1(patch->totalDofsPerCell * numExtFacets, &facetDofsWithAll));
2539: for (PetscInt i = 0; i < numExtFacets; i++) {
2540: const PetscInt cell = extFacetCells[extFacetOffset + i];
2541: for (PetscInt d = 0; d < patch->totalDofsPerCell; d++) {
2542: facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2543: if (dofsArrayWithAll) facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell) * patch->totalDofsPerCell + d];
2544: idx++;
2545: }
2546: }
2547: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray + extFacetOffset, PETSC_USE_POINTER, &facetIS));
2548: PetscCall(patch->usercomputeopextfacet(pc, point, x, mat, facetIS, numExtFacets * patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopextfacetctx));
2549: PetscCall(ISDestroy(&facetIS));
2550: PetscCall(ISRestoreIndices(patch->extFacetsToPatchCell, &extFacetCells));
2551: PetscCall(ISRestoreIndices(patch->extFacets, &extFacetsArray));
2552: PetscCall(PetscFree(facetDofs));
2553: PetscCall(PetscFree(facetDofsWithAll));
2554: }
2555: }
2557: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2558: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2560: if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2561: MatFactorInfo info;
2562: PetscBool flg;
2563: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg));
2564: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Invalid Mat type for dense inverse");
2565: PetscCall(MatFactorInfoInitialize(&info));
2566: PetscCall(MatLUFactor(mat, NULL, NULL, &info));
2567: PetscCall(MatSeqDenseInvertFactors_Private(mat));
2568: }
2569: PetscCall(ISDestroy(&patch->cellIS));
2570: if (withArtificial) {
2571: PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
2572: } else {
2573: PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
2574: }
2575: if (isNonlinear) PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
2576: PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
2577: if (patch->viewMatrix) {
2578: char name[PETSC_MAX_PATH_LEN];
2580: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch matrix for Point %" PetscInt_FMT, point));
2581: PetscCall(PetscObjectSetName((PetscObject)mat, name));
2582: PetscCall(ObjectView((PetscObject)mat, patch->viewerMatrix, patch->formatMatrix));
2583: }
2584: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2585: PetscFunctionReturn(PETSC_SUCCESS);
2586: }
2588: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2589: {
2590: Vec data;
2591: PetscScalar *array;
2592: PetscInt bs, nz, i, j, cell;
2594: PetscFunctionBegin;
2595: PetscCall(MatShellGetContext(mat, &data));
2596: PetscCall(VecGetBlockSize(data, &bs));
2597: PetscCall(VecGetSize(data, &nz));
2598: PetscCall(VecGetArray(data, &array));
2599: PetscCheck(m == n, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2600: cell = idxm[0] / bs; /* use the fact that this is called once per cell */
2601: for (i = 0; i < m; i++) {
2602: PetscCheck(idxm[i] == idxn[i], PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2603: for (j = 0; j < n; j++) {
2604: const PetscScalar v_ = v[i * bs + j];
2605: /* Indexing is special to the data structure we have! */
2606: if (addv == INSERT_VALUES) {
2607: array[cell * bs * bs + i * bs + j] = v_;
2608: } else {
2609: array[cell * bs * bs + i * bs + j] += v_;
2610: }
2611: }
2612: }
2613: PetscCall(VecRestoreArray(data, &array));
2614: PetscFunctionReturn(PETSC_SUCCESS);
2615: }
2617: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2618: {
2619: PC_PATCH *patch = (PC_PATCH *)pc->data;
2620: const PetscInt *cellsArray;
2621: PetscInt ncell, offset;
2622: const PetscInt *dofMapArray;
2623: PetscInt i, j;
2624: IS dofMap;
2625: IS cellIS;
2626: const PetscInt ndof = patch->totalDofsPerCell;
2627: Mat vecMat;
2628: PetscInt cStart, cEnd;
2629: DM dm, plex;
2631: PetscFunctionBegin;
2632: PetscCall(ISGetSize(patch->cells, &ncell));
2633: if (!ncell) { /* No cells to assemble over -> skip */
2634: PetscFunctionReturn(PETSC_SUCCESS);
2635: }
2637: PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
2639: PetscCall(PCGetDM(pc, &dm));
2640: PetscCall(DMConvert(dm, DMPLEX, &plex));
2641: dm = plex;
2642: if (!patch->allCells) {
2643: PetscHSetI cells;
2644: PetscHashIter hi;
2645: PetscInt pStart, pEnd;
2646: PetscInt *allCells = NULL;
2647: PetscCall(PetscHSetICreate(&cells));
2648: PetscCall(ISGetIndices(patch->cells, &cellsArray));
2649: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
2650: for (i = pStart; i < pEnd; i++) {
2651: PetscCall(PetscSectionGetDof(patch->cellCounts, i, &ncell));
2652: PetscCall(PetscSectionGetOffset(patch->cellCounts, i, &offset));
2653: if (ncell <= 0) continue;
2654: for (j = 0; j < ncell; j++) PetscCall(PetscHSetIAdd(cells, cellsArray[offset + j]));
2655: }
2656: PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
2657: PetscCall(PetscHSetIGetSize(cells, &ncell));
2658: PetscCall(PetscMalloc1(ncell, &allCells));
2659: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2660: PetscCall(PetscMalloc1(cEnd - cStart, &patch->precomputedTensorLocations));
2661: i = 0;
2662: PetscHashIterBegin(cells, hi);
2663: while (!PetscHashIterAtEnd(cells, hi)) {
2664: PetscHashIterGetKey(cells, hi, allCells[i]);
2665: patch->precomputedTensorLocations[allCells[i]] = i;
2666: PetscHashIterNext(cells, hi);
2667: i++;
2668: }
2669: PetscCall(PetscHSetIDestroy(&cells));
2670: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells));
2671: }
2672: PetscCall(ISGetSize(patch->allCells, &ncell));
2673: if (!patch->cellMats) {
2674: PetscCall(VecCreateSeq(PETSC_COMM_SELF, ncell * ndof * ndof, &patch->cellMats));
2675: PetscCall(VecSetBlockSize(patch->cellMats, ndof));
2676: }
2677: PetscCall(VecSet(patch->cellMats, 0));
2679: PetscCall(MatCreateShell(PETSC_COMM_SELF, ncell * ndof, ncell * ndof, ncell * ndof, ncell * ndof, (void *)patch->cellMats, &vecMat));
2680: PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (PetscErrorCodeFn *)MatSetValues_PCPatch_Private));
2681: PetscCall(ISGetSize(patch->allCells, &ncell));
2682: PetscCall(ISCreateStride(PETSC_COMM_SELF, ndof * ncell, 0, 1, &dofMap));
2683: PetscCall(ISGetIndices(dofMap, &dofMapArray));
2684: PetscCall(ISGetIndices(patch->allCells, &cellsArray));
2685: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS));
2686: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2687: PetscCallBack("PCPatch callback", patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof * ncell, dofMapArray, NULL, patch->usercomputeopctx));
2688: PetscCall(ISDestroy(&cellIS));
2689: PetscCall(MatDestroy(&vecMat));
2690: PetscCall(ISRestoreIndices(patch->allCells, &cellsArray));
2691: PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
2692: PetscCall(ISDestroy(&dofMap));
2694: if (patch->usercomputeopintfacet) {
2695: PetscInt nIntFacets;
2696: IS intFacetsIS;
2697: const PetscInt *intFacetsArray = NULL;
2698: if (!patch->allIntFacets) {
2699: PetscHSetI facets;
2700: PetscHashIter hi;
2701: PetscInt pStart, pEnd, fStart, fEnd;
2702: PetscInt *allIntFacets = NULL;
2703: PetscCall(PetscHSetICreate(&facets));
2704: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2705: PetscCall(PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd));
2706: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2707: for (i = pStart; i < pEnd; i++) {
2708: PetscCall(PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets));
2709: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, i, &offset));
2710: if (nIntFacets <= 0) continue;
2711: for (j = 0; j < nIntFacets; j++) PetscCall(PetscHSetIAdd(facets, intFacetsArray[offset + j]));
2712: }
2713: PetscCall(ISRestoreIndices(patch->intFacets, &intFacetsArray));
2714: PetscCall(PetscHSetIGetSize(facets, &nIntFacets));
2715: PetscCall(PetscMalloc1(nIntFacets, &allIntFacets));
2716: PetscCall(PetscMalloc1(fEnd - fStart, &patch->precomputedIntFacetTensorLocations));
2717: i = 0;
2718: PetscHashIterBegin(facets, hi);
2719: while (!PetscHashIterAtEnd(facets, hi)) {
2720: PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2721: patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2722: PetscHashIterNext(facets, hi);
2723: i++;
2724: }
2725: PetscCall(PetscHSetIDestroy(&facets));
2726: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets));
2727: }
2728: PetscCall(ISGetSize(patch->allIntFacets, &nIntFacets));
2729: if (!patch->intFacetMats) {
2730: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nIntFacets * ndof * ndof * 4, &patch->intFacetMats));
2731: PetscCall(VecSetBlockSize(patch->intFacetMats, ndof * 2));
2732: }
2733: PetscCall(VecSet(patch->intFacetMats, 0));
2735: PetscCall(MatCreateShell(PETSC_COMM_SELF, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, (void *)patch->intFacetMats, &vecMat));
2736: PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (PetscErrorCodeFn *)MatSetValues_PCPatch_Private));
2737: PetscCall(ISCreateStride(PETSC_COMM_SELF, 2 * ndof * nIntFacets, 0, 1, &dofMap));
2738: PetscCall(ISGetIndices(dofMap, &dofMapArray));
2739: PetscCall(ISGetIndices(patch->allIntFacets, &intFacetsArray));
2740: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS));
2741: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2742: PetscCallBack("PCPatch callback (interior facets)", patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2 * ndof * nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx));
2743: PetscCall(ISDestroy(&intFacetsIS));
2744: PetscCall(MatDestroy(&vecMat));
2745: PetscCall(ISRestoreIndices(patch->allIntFacets, &intFacetsArray));
2746: PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
2747: PetscCall(ISDestroy(&dofMap));
2748: }
2749: PetscCall(DMDestroy(&dm));
2750: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2751: PetscFunctionReturn(PETSC_SUCCESS);
2752: }
2754: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2755: {
2756: PC_PATCH *patch = (PC_PATCH *)pc->data;
2757: const PetscScalar *xArray = NULL;
2758: PetscScalar *yArray = NULL;
2759: const PetscInt *gtolArray = NULL;
2760: PetscInt dof, offset, lidx;
2762: PetscFunctionBeginHot;
2763: PetscCall(VecGetArrayRead(x, &xArray));
2764: PetscCall(VecGetArray(y, &yArray));
2765: if (scattertype == SCATTER_WITHARTIFICIAL) {
2766: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof));
2767: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset));
2768: PetscCall(ISGetIndices(patch->gtolWithArtificial, >olArray));
2769: } else if (scattertype == SCATTER_WITHALL) {
2770: PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof));
2771: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset));
2772: PetscCall(ISGetIndices(patch->gtolWithAll, >olArray));
2773: } else {
2774: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
2775: PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2776: PetscCall(ISGetIndices(patch->gtol, >olArray));
2777: }
2778: PetscCheck(mode != INSERT_VALUES || scat == SCATTER_FORWARD, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward");
2779: PetscCheck(mode != ADD_VALUES || scat == SCATTER_REVERSE, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse");
2780: for (lidx = 0; lidx < dof; ++lidx) {
2781: const PetscInt gidx = gtolArray[offset + lidx];
2783: if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
2784: else yArray[gidx] += xArray[lidx]; /* Reverse */
2785: }
2786: if (scattertype == SCATTER_WITHARTIFICIAL) {
2787: PetscCall(ISRestoreIndices(patch->gtolWithArtificial, >olArray));
2788: } else if (scattertype == SCATTER_WITHALL) {
2789: PetscCall(ISRestoreIndices(patch->gtolWithAll, >olArray));
2790: } else {
2791: PetscCall(ISRestoreIndices(patch->gtol, >olArray));
2792: }
2793: PetscCall(VecRestoreArrayRead(x, &xArray));
2794: PetscCall(VecRestoreArray(y, &yArray));
2795: PetscFunctionReturn(PETSC_SUCCESS);
2796: }
2798: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2799: {
2800: PC_PATCH *patch = (PC_PATCH *)pc->data;
2801: const char *prefix;
2802: PetscInt i;
2804: PetscFunctionBegin;
2805: if (!pc->setupcalled) {
2806: PetscCheck(patch->save_operators || !patch->denseinverse, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
2807: if (!patch->denseinverse) {
2808: PetscCall(PetscMalloc1(patch->npatch, &patch->solver));
2809: PetscCall(PCGetOptionsPrefix(pc, &prefix));
2810: for (i = 0; i < patch->npatch; ++i) {
2811: KSP ksp;
2812: PC subpc;
2814: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
2815: PetscCall(KSPSetNestLevel(ksp, pc->kspnestlevel));
2816: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
2817: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
2818: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
2819: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
2820: PetscCall(KSPGetPC(ksp, &subpc));
2821: PetscCall(PetscObjectIncrementTabLevel((PetscObject)subpc, (PetscObject)pc, 1));
2822: patch->solver[i] = (PetscObject)ksp;
2823: }
2824: }
2825: }
2826: if (patch->save_operators) {
2827: if (patch->precomputeElementTensors) PetscCall(PCPatchPrecomputePatchTensors_Private(pc));
2828: for (i = 0; i < patch->npatch; ++i) {
2829: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE));
2830: if (!patch->denseinverse) {
2831: PetscCall(KSPSetOperators((KSP)patch->solver[i], patch->mat[i], patch->mat[i]));
2832: } else if (patch->mat[i] && !patch->densesolve) {
2833: /* Setup matmult callback */
2834: PetscCall(MatGetOperation(patch->mat[i], MATOP_MULT, (PetscErrorCodeFn **)&patch->densesolve));
2835: }
2836: }
2837: }
2838: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2839: for (i = 0; i < patch->npatch; ++i) {
2840: /* Instead of padding patch->patchUpdate with zeros to get */
2841: /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2842: /* just get rid of the columns that correspond to the dofs with */
2843: /* artificial bcs. That's of course fairly inefficient, hopefully we */
2844: /* can just assemble the rectangular matrix in the first place. */
2845: Mat matSquare;
2846: IS rowis;
2847: PetscInt dof;
2849: PetscCall(MatGetSize(patch->mat[i], &dof, NULL));
2850: if (dof == 0) {
2851: patch->matWithArtificial[i] = NULL;
2852: continue;
2853: }
2855: PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
2856: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));
2858: PetscCall(MatGetSize(matSquare, &dof, NULL));
2859: PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
2860: if (pc->setupcalled) {
2861: PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]));
2862: } else {
2863: PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]));
2864: }
2865: PetscCall(ISDestroy(&rowis));
2866: PetscCall(MatDestroy(&matSquare));
2867: }
2868: }
2869: PetscFunctionReturn(PETSC_SUCCESS);
2870: }
2872: static PetscErrorCode PCSetUp_PATCH(PC pc)
2873: {
2874: PC_PATCH *patch = (PC_PATCH *)pc->data;
2875: PetscInt i;
2876: PetscBool isNonlinear;
2877: PetscInt maxDof = -1, maxDofWithArtificial = -1;
2879: PetscFunctionBegin;
2880: if (!pc->setupcalled) {
2881: PetscInt pStart, pEnd, p;
2882: PetscInt localSize;
2884: PetscCall(PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0));
2886: isNonlinear = patch->isNonlinear;
2887: if (!patch->nsubspaces) {
2888: DM dm, plex;
2889: PetscSection s;
2890: PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;
2892: PetscCall(PCGetDM(pc, &dm));
2893: PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2894: PetscCall(DMConvert(dm, DMPLEX, &plex));
2895: dm = plex;
2896: PetscCall(DMGetLocalSection(dm, &s));
2897: PetscCall(PetscSectionGetNumFields(s, &Nf));
2898: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2899: for (p = pStart; p < pEnd; ++p) {
2900: PetscInt cdof;
2901: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2902: numGlobalBcs += cdof;
2903: }
2904: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2905: PetscCall(PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs));
2906: for (f = 0; f < Nf; ++f) {
2907: PetscFE fe;
2908: PetscDualSpace sp;
2909: PetscInt cdoff = 0;
2911: PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
2912: /* PetscCall(PetscFEGetNumComponents(fe, &Nc[f])); */
2913: PetscCall(PetscFEGetDualSpace(fe, &sp));
2914: PetscCall(PetscDualSpaceGetDimension(sp, &Nb[f]));
2916: PetscCall(PetscMalloc1((cEnd - cStart) * Nb[f], &cellDofs[f]));
2917: for (c = cStart; c < cEnd; ++c) {
2918: PetscInt *closure = NULL;
2919: PetscInt clSize = 0, cl;
2921: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2922: for (cl = 0; cl < clSize * 2; cl += 2) {
2923: const PetscInt p = closure[cl];
2924: PetscInt fdof, d, foff;
2926: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2927: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2928: for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2929: }
2930: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2931: }
2932: PetscCheck(cdoff == (cEnd - cStart) * Nb[f], PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "Total number of cellDofs %" PetscInt_FMT " for field %" PetscInt_FMT " should be Nc (%" PetscInt_FMT ") * cellDof (%" PetscInt_FMT ")", cdoff, f, cEnd - cStart, Nb[f]);
2933: }
2934: numGlobalBcs = 0;
2935: for (p = pStart; p < pEnd; ++p) {
2936: const PetscInt *ind;
2937: PetscInt off, cdof, d;
2939: PetscCall(PetscSectionGetOffset(s, p, &off));
2940: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2941: PetscCall(PetscSectionGetConstraintIndices(s, p, &ind));
2942: for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2943: }
2945: PetscCall(PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **)cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs));
2946: for (f = 0; f < Nf; ++f) PetscCall(PetscFree(cellDofs[f]));
2947: PetscCall(PetscFree3(Nb, cellDofs, globalBcs));
2948: PetscCall(PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL));
2949: PetscCall(PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL));
2950: PetscCall(DMDestroy(&dm));
2951: }
2953: localSize = patch->subspaceOffsets[patch->nsubspaces];
2954: PetscCall(VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS));
2955: PetscCall(VecSetUp(patch->localRHS));
2956: PetscCall(VecDuplicate(patch->localRHS, &patch->localUpdate));
2957: PetscCall(PCPatchCreateCellPatches(pc));
2958: PetscCall(PCPatchCreateCellPatchDiscretisationInfo(pc));
2960: /* OK, now build the work vectors */
2961: PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd));
2963: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial));
2964: if (isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll));
2965: for (p = pStart; p < pEnd; ++p) {
2966: PetscInt dof;
2968: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
2969: maxDof = PetscMax(maxDof, dof);
2970: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2971: const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
2972: PetscInt numPatchDofs, offset;
2973: PetscInt numPatchDofsWithArtificial, offsetWithArtificial;
2974: PetscInt dofWithoutArtificialCounter = 0;
2975: PetscInt *patchWithoutArtificialToWithArtificialArray;
2977: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof));
2978: maxDofWithArtificial = PetscMax(maxDofWithArtificial, dof);
2980: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
2981: /* the index in the patch with all dofs */
2982: PetscCall(ISGetIndices(patch->gtol, >olArray));
2984: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
2985: if (numPatchDofs == 0) {
2986: patch->dofMappingWithoutToWithArtificial[p - pStart] = NULL;
2987: continue;
2988: }
2990: PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2991: PetscCall(ISGetIndices(patch->gtolWithArtificial, >olArrayWithArtificial));
2992: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial));
2993: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial));
2995: PetscCall(PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray));
2996: for (i = 0; i < numPatchDofsWithArtificial; i++) {
2997: if (gtolArrayWithArtificial[i + offsetWithArtificial] == gtolArray[offset + dofWithoutArtificialCounter]) {
2998: patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2999: dofWithoutArtificialCounter++;
3000: if (dofWithoutArtificialCounter == numPatchDofs) break;
3001: }
3002: }
3003: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p - pStart]));
3004: PetscCall(ISRestoreIndices(patch->gtol, >olArray));
3005: PetscCall(ISRestoreIndices(patch->gtolWithArtificial, >olArrayWithArtificial));
3006: }
3007: }
3008: for (p = pStart; p < pEnd; ++p) {
3009: if (isNonlinear) {
3010: const PetscInt *gtolArray, *gtolArrayWithAll = NULL;
3011: PetscInt numPatchDofs, offset;
3012: PetscInt numPatchDofsWithAll, offsetWithAll;
3013: PetscInt dofWithoutAllCounter = 0;
3014: PetscInt *patchWithoutAllToWithAllArray;
3016: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
3017: /* the index in the patch with all dofs */
3018: PetscCall(ISGetIndices(patch->gtol, >olArray));
3020: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
3021: if (numPatchDofs == 0) {
3022: patch->dofMappingWithoutToWithAll[p - pStart] = NULL;
3023: continue;
3024: }
3026: PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
3027: PetscCall(ISGetIndices(patch->gtolWithAll, >olArrayWithAll));
3028: PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll));
3029: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll));
3031: PetscCall(PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray));
3033: for (i = 0; i < numPatchDofsWithAll; i++) {
3034: if (gtolArrayWithAll[i + offsetWithAll] == gtolArray[offset + dofWithoutAllCounter]) {
3035: patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
3036: dofWithoutAllCounter++;
3037: if (dofWithoutAllCounter == numPatchDofs) break;
3038: }
3039: }
3040: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p - pStart]));
3041: PetscCall(ISRestoreIndices(patch->gtol, >olArray));
3042: PetscCall(ISRestoreIndices(patch->gtolWithAll, >olArrayWithAll));
3043: }
3044: }
3045: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3046: PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial));
3047: PetscCall(VecSetUp(patch->patchRHSWithArtificial));
3048: }
3049: PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS));
3050: PetscCall(VecSetUp(patch->patchRHS));
3051: PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate));
3052: PetscCall(VecSetUp(patch->patchUpdate));
3053: if (patch->save_operators) {
3054: PetscCall(PetscMalloc1(patch->npatch, &patch->mat));
3055: for (i = 0; i < patch->npatch; ++i) PetscCall(PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE));
3056: }
3057: PetscCall(PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0));
3059: /* If desired, calculate weights for dof multiplicity */
3060: if (patch->partition_of_unity) {
3061: PetscScalar *input = NULL;
3062: PetscScalar *output = NULL;
3063: Vec global;
3065: PetscCall(VecDuplicate(patch->localRHS, &patch->dof_weights));
3066: if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
3067: for (i = 0; i < patch->npatch; ++i) {
3068: PetscInt dof;
3070: PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &dof));
3071: if (dof <= 0) continue;
3072: PetscCall(VecSet(patch->patchRHS, 1.0));
3073: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
3074: }
3075: } else {
3076: /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
3077: PetscCall(VecSet(patch->dof_weights, 1.0));
3078: }
3080: PetscCall(VecDuplicate(patch->dof_weights, &global));
3081: PetscCall(VecSet(global, 0.));
3083: PetscCall(VecGetArray(patch->dof_weights, &input));
3084: PetscCall(VecGetArray(global, &output));
3085: PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
3086: PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
3087: PetscCall(VecRestoreArray(patch->dof_weights, &input));
3088: PetscCall(VecRestoreArray(global, &output));
3090: PetscCall(VecReciprocal(global));
3092: PetscCall(VecGetArray(patch->dof_weights, &output));
3093: PetscCall(VecGetArray(global, &input));
3094: PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
3095: PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
3096: PetscCall(VecRestoreArray(patch->dof_weights, &output));
3097: PetscCall(VecRestoreArray(global, &input));
3098: PetscCall(VecDestroy(&global));
3099: }
3100: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators && !patch->isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->matWithArtificial));
3101: }
3102: PetscCall((*patch->setupsolver)(pc));
3103: PetscFunctionReturn(PETSC_SUCCESS);
3104: }
3106: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
3107: {
3108: PC_PATCH *patch = (PC_PATCH *)pc->data;
3109: KSP ksp;
3110: Mat op;
3111: PetscInt m, n;
3113: PetscFunctionBegin;
3114: if (patch->denseinverse) {
3115: PetscCall((*patch->densesolve)(patch->mat[i], x, y));
3116: PetscFunctionReturn(PETSC_SUCCESS);
3117: }
3118: ksp = (KSP)patch->solver[i];
3119: if (!patch->save_operators) {
3120: Mat mat;
3122: PetscCall(PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE));
3123: /* Populate operator here. */
3124: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE));
3125: PetscCall(KSPSetOperators(ksp, mat, mat));
3126: /* Drop reference so the KSPSetOperators below will blow it away. */
3127: PetscCall(MatDestroy(&mat));
3128: }
3129: PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
3130: if (!ksp->setfromoptionscalled) PetscCall(KSPSetFromOptions(ksp));
3131: /* Disgusting trick to reuse work vectors */
3132: PetscCall(KSPGetOperators(ksp, &op, NULL));
3133: PetscCall(MatGetLocalSize(op, &m, &n));
3134: x->map->n = m;
3135: y->map->n = n;
3136: x->map->N = m;
3137: y->map->N = n;
3138: x->map->setupcalled = PETSC_FALSE;
3139: y->map->setupcalled = PETSC_FALSE;
3140: PetscCall(KSPSolve(ksp, x, y));
3141: PetscCall(KSPCheckSolve(ksp, pc, y));
3142: PetscCall(PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0));
3143: if (!patch->save_operators) {
3144: PC pc;
3145: PetscCall(KSPSetOperators(ksp, NULL, NULL));
3146: PetscCall(KSPGetPC(ksp, &pc));
3147: /* Destroy PC context too, otherwise the factored matrix hangs around. */
3148: PetscCall(PCReset(pc));
3149: }
3150: PetscFunctionReturn(PETSC_SUCCESS);
3151: }
3153: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
3154: {
3155: PC_PATCH *patch = (PC_PATCH *)pc->data;
3156: Mat multMat;
3157: PetscInt n, m;
3159: PetscFunctionBegin;
3160: if (patch->save_operators) {
3161: multMat = patch->matWithArtificial[i];
3162: } else {
3163: /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
3164: Mat matSquare;
3165: PetscInt dof;
3166: IS rowis;
3167: PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
3168: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));
3169: PetscCall(MatGetSize(matSquare, &dof, NULL));
3170: PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
3171: PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat));
3172: PetscCall(MatDestroy(&matSquare));
3173: PetscCall(ISDestroy(&rowis));
3174: }
3175: /* Disgusting trick to reuse work vectors */
3176: PetscCall(MatGetLocalSize(multMat, &m, &n));
3177: patch->patchUpdate->map->n = n;
3178: patch->patchRHSWithArtificial->map->n = m;
3179: patch->patchUpdate->map->N = n;
3180: patch->patchRHSWithArtificial->map->N = m;
3181: patch->patchUpdate->map->setupcalled = PETSC_FALSE;
3182: patch->patchRHSWithArtificial->map->setupcalled = PETSC_FALSE;
3183: PetscCall(MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial));
3184: PetscCall(VecScale(patch->patchRHSWithArtificial, -1.0));
3185: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL));
3186: if (!patch->save_operators) PetscCall(MatDestroy(&multMat));
3187: PetscFunctionReturn(PETSC_SUCCESS);
3188: }
3190: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
3191: {
3192: PC_PATCH *patch = (PC_PATCH *)pc->data;
3193: const PetscScalar *globalRHS = NULL;
3194: PetscScalar *localRHS = NULL;
3195: PetscScalar *globalUpdate = NULL;
3196: const PetscInt *bcNodes = NULL;
3197: PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1;
3198: PetscInt start[2] = {0, 0};
3199: PetscInt end[2] = {-1, -1};
3200: const PetscInt inc[2] = {1, -1};
3201: const PetscScalar *localUpdate;
3202: const PetscInt *iterationSet;
3203: PetscInt pStart, numBcs, n, sweep, bc, j;
3205: PetscFunctionBegin;
3206: PetscCall(PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0));
3207: PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
3208: /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
3209: end[0] = patch->npatch;
3210: start[1] = patch->npatch - 1;
3211: if (patch->user_patches) {
3212: PetscCall(ISGetLocalSize(patch->iterationSet, &end[0]));
3213: start[1] = end[0] - 1;
3214: PetscCall(ISGetIndices(patch->iterationSet, &iterationSet));
3215: }
3216: /* Scatter from global space into overlapped local spaces */
3217: PetscCall(VecGetArrayRead(x, &globalRHS));
3218: PetscCall(VecGetArray(patch->localRHS, &localRHS));
3219: PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
3220: PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
3221: PetscCall(VecRestoreArrayRead(x, &globalRHS));
3222: PetscCall(VecRestoreArray(patch->localRHS, &localRHS));
3224: PetscCall(VecSet(patch->localUpdate, 0.0));
3225: PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
3226: PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
3227: for (sweep = 0; sweep < nsweep; sweep++) {
3228: for (j = start[sweep]; j * inc[sweep] < end[sweep] * inc[sweep]; j += inc[sweep]) {
3229: PetscInt i = patch->user_patches ? iterationSet[j] : j;
3230: PetscInt start, len;
3232: PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &len));
3233: PetscCall(PetscSectionGetOffset(patch->gtolCounts, i + pStart, &start));
3234: /* TODO: Squash out these guys in the setup as well. */
3235: if (len <= 0) continue;
3236: /* TODO: Do we need different scatters for X and Y? */
3237: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->localRHS, patch->patchRHS, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR));
3238: PetscCall((*patch->applysolver)(pc, i, patch->patchRHS, patch->patchUpdate));
3239: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchUpdate, patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
3240: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall((*patch->updatemultiplicative)(pc, i, pStart));
3241: }
3242: }
3243: PetscCall(PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0));
3244: if (patch->user_patches) PetscCall(ISRestoreIndices(patch->iterationSet, &iterationSet));
3245: /* XXX: should we do this on the global vector? */
3246: if (patch->partition_of_unity) PetscCall(VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights));
3247: /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
3248: PetscCall(VecSet(y, 0.0));
3249: PetscCall(VecGetArray(y, &globalUpdate));
3250: PetscCall(VecGetArrayRead(patch->localUpdate, &localUpdate));
3251: PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM));
3252: PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM));
3253: PetscCall(VecRestoreArrayRead(patch->localUpdate, &localUpdate));
3255: /* Now we need to send the global BC values through */
3256: PetscCall(VecGetArrayRead(x, &globalRHS));
3257: PetscCall(ISGetSize(patch->globalBcNodes, &numBcs));
3258: PetscCall(ISGetIndices(patch->globalBcNodes, &bcNodes));
3259: PetscCall(VecGetLocalSize(x, &n));
3260: for (bc = 0; bc < numBcs; ++bc) {
3261: const PetscInt idx = bcNodes[bc];
3262: if (idx < n) globalUpdate[idx] = globalRHS[idx];
3263: }
3265: PetscCall(ISRestoreIndices(patch->globalBcNodes, &bcNodes));
3266: PetscCall(VecRestoreArrayRead(x, &globalRHS));
3267: PetscCall(VecRestoreArray(y, &globalUpdate));
3269: PetscCall(PetscOptionsPopCreateViewerOff());
3270: PetscCall(PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0));
3271: PetscFunctionReturn(PETSC_SUCCESS);
3272: }
3274: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
3275: {
3276: PC_PATCH *patch = (PC_PATCH *)pc->data;
3277: PetscInt i;
3279: PetscFunctionBegin;
3280: if (patch->solver) {
3281: for (i = 0; i < patch->npatch; ++i) PetscCall(KSPReset((KSP)patch->solver[i]));
3282: }
3283: PetscFunctionReturn(PETSC_SUCCESS);
3284: }
3286: static PetscErrorCode PCReset_PATCH(PC pc)
3287: {
3288: PC_PATCH *patch = (PC_PATCH *)pc->data;
3289: PetscInt i;
3291: PetscFunctionBegin;
3292: PetscCall(PetscSFDestroy(&patch->sectionSF));
3293: PetscCall(PetscSectionDestroy(&patch->cellCounts));
3294: PetscCall(PetscSectionDestroy(&patch->pointCounts));
3295: PetscCall(PetscSectionDestroy(&patch->cellNumbering));
3296: PetscCall(PetscSectionDestroy(&patch->gtolCounts));
3297: PetscCall(ISDestroy(&patch->gtol));
3298: PetscCall(ISDestroy(&patch->cells));
3299: PetscCall(ISDestroy(&patch->points));
3300: PetscCall(ISDestroy(&patch->dofs));
3301: PetscCall(ISDestroy(&patch->offs));
3302: PetscCall(PetscSectionDestroy(&patch->patchSection));
3303: PetscCall(ISDestroy(&patch->ghostBcNodes));
3304: PetscCall(ISDestroy(&patch->globalBcNodes));
3305: PetscCall(PetscSectionDestroy(&patch->gtolCountsWithArtificial));
3306: PetscCall(ISDestroy(&patch->gtolWithArtificial));
3307: PetscCall(ISDestroy(&patch->dofsWithArtificial));
3308: PetscCall(ISDestroy(&patch->offsWithArtificial));
3309: PetscCall(PetscSectionDestroy(&patch->gtolCountsWithAll));
3310: PetscCall(ISDestroy(&patch->gtolWithAll));
3311: PetscCall(ISDestroy(&patch->dofsWithAll));
3312: PetscCall(ISDestroy(&patch->offsWithAll));
3313: PetscCall(VecDestroy(&patch->cellMats));
3314: PetscCall(VecDestroy(&patch->intFacetMats));
3315: PetscCall(ISDestroy(&patch->allCells));
3316: PetscCall(ISDestroy(&patch->intFacets));
3317: PetscCall(ISDestroy(&patch->extFacets));
3318: PetscCall(ISDestroy(&patch->intFacetsToPatchCell));
3319: PetscCall(ISDestroy(&patch->extFacetsToPatchCell));
3320: PetscCall(PetscSectionDestroy(&patch->intFacetCounts));
3321: PetscCall(PetscSectionDestroy(&patch->extFacetCounts));
3323: if (patch->dofSection)
3324: for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscSectionDestroy(&patch->dofSection[i]));
3325: PetscCall(PetscFree(patch->dofSection));
3326: PetscCall(PetscFree(patch->bs));
3327: PetscCall(PetscFree(patch->nodesPerCell));
3328: if (patch->cellNodeMap)
3329: for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscFree(patch->cellNodeMap[i]));
3330: PetscCall(PetscFree(patch->cellNodeMap));
3331: PetscCall(PetscFree(patch->subspaceOffsets));
3333: PetscCall((*patch->resetsolver)(pc));
3335: PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));
3337: PetscCall(VecDestroy(&patch->localRHS));
3338: PetscCall(VecDestroy(&patch->localUpdate));
3339: PetscCall(VecDestroy(&patch->patchRHS));
3340: PetscCall(VecDestroy(&patch->patchUpdate));
3341: PetscCall(VecDestroy(&patch->dof_weights));
3342: if (patch->patch_dof_weights) {
3343: for (i = 0; i < patch->npatch; ++i) PetscCall(VecDestroy(&patch->patch_dof_weights[i]));
3344: PetscCall(PetscFree(patch->patch_dof_weights));
3345: }
3346: if (patch->mat) {
3347: for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->mat[i]));
3348: PetscCall(PetscFree(patch->mat));
3349: }
3350: if (patch->matWithArtificial && !patch->isNonlinear) {
3351: for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->matWithArtificial[i]));
3352: PetscCall(PetscFree(patch->matWithArtificial));
3353: }
3354: PetscCall(VecDestroy(&patch->patchRHSWithArtificial));
3355: if (patch->dofMappingWithoutToWithArtificial) {
3356: for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]));
3357: PetscCall(PetscFree(patch->dofMappingWithoutToWithArtificial));
3358: }
3359: if (patch->dofMappingWithoutToWithAll) {
3360: for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithAll[i]));
3361: PetscCall(PetscFree(patch->dofMappingWithoutToWithAll));
3362: }
3363: PetscCall(PetscFree(patch->sub_mat_type));
3364: if (patch->userIS) {
3365: for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->userIS[i]));
3366: PetscCall(PetscFree(patch->userIS));
3367: }
3368: PetscCall(PetscFree(patch->precomputedTensorLocations));
3369: PetscCall(PetscFree(patch->precomputedIntFacetTensorLocations));
3371: patch->bs = NULL;
3372: patch->cellNodeMap = NULL;
3373: patch->nsubspaces = 0;
3374: PetscCall(ISDestroy(&patch->iterationSet));
3376: PetscCall(PetscViewerDestroy(&patch->viewerCells));
3377: PetscCall(PetscViewerDestroy(&patch->viewerIntFacets));
3378: PetscCall(PetscViewerDestroy(&patch->viewerPoints));
3379: PetscCall(PetscViewerDestroy(&patch->viewerSection));
3380: PetscCall(PetscViewerDestroy(&patch->viewerMatrix));
3381: PetscFunctionReturn(PETSC_SUCCESS);
3382: }
3384: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
3385: {
3386: PC_PATCH *patch = (PC_PATCH *)pc->data;
3387: PetscInt i;
3389: PetscFunctionBegin;
3390: if (patch->solver) {
3391: for (i = 0; i < patch->npatch; ++i) PetscCall(KSPDestroy((KSP *)&patch->solver[i]));
3392: PetscCall(PetscFree(patch->solver));
3393: }
3394: PetscFunctionReturn(PETSC_SUCCESS);
3395: }
3397: static PetscErrorCode PCDestroy_PATCH(PC pc)
3398: {
3399: PC_PATCH *patch = (PC_PATCH *)pc->data;
3401: PetscFunctionBegin;
3402: PetscCall(PCReset_PATCH(pc));
3403: PetscCall((*patch->destroysolver)(pc));
3404: PetscCall(PetscFree(pc->data));
3405: PetscFunctionReturn(PETSC_SUCCESS);
3406: }
3408: static PetscErrorCode PCSetFromOptions_PATCH(PC pc, PetscOptionItems PetscOptionsObject)
3409: {
3410: PC_PATCH *patch = (PC_PATCH *)pc->data;
3411: PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
3412: char sub_mat_type[PETSC_MAX_PATH_LEN];
3413: char option[PETSC_MAX_PATH_LEN];
3414: const char *prefix;
3415: PetscBool flg, dimflg, codimflg;
3416: MPI_Comm comm;
3417: PetscInt *ifields, nfields, k;
3418: PCCompositeType loctype = PC_COMPOSITE_ADDITIVE;
3420: PetscFunctionBegin;
3421: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
3422: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
3423: PetscOptionsHeadBegin(PetscOptionsObject, "Patch solver options");
3425: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname));
3426: PetscCall(PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg));
3428: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_precompute_element_tensors", patch->classname));
3429: PetscCall(PetscOptionsBool(option, "Compute each element tensor only once?", "PCPatchSetPrecomputeElementTensors", patch->precomputeElementTensors, &patch->precomputeElementTensors, &flg));
3430: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname));
3431: PetscCall(PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg));
3433: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname));
3434: PetscCall(PetscOptionsEnum(option, "Type of local solver composition (additive or multiplicative)", "PCPatchSetLocalComposition", PCCompositeTypes, (PetscEnum)loctype, (PetscEnum *)&loctype, &flg));
3435: if (flg) PetscCall(PCPatchSetLocalComposition(pc, loctype));
3436: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_dense_inverse", patch->classname));
3437: PetscCall(PetscOptionsBool(option, "Compute inverses of patch matrices and apply directly? Ignores KSP/PC settings on patch.", "PCPatchSetDenseInverse", patch->denseinverse, &patch->denseinverse, &flg));
3438: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname));
3439: PetscCall(PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg));
3440: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname));
3441: PetscCall(PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg));
3442: PetscCheck(!dimflg || !codimflg, comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");
3444: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname));
3445: PetscCall(PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum)patchConstructionType, (PetscEnum *)&patchConstructionType, &flg));
3446: if (flg) PetscCall(PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL));
3448: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname));
3449: PetscCall(PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg));
3451: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname));
3452: PetscCall(PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg));
3454: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_pardecomp_overlap", patch->classname));
3455: PetscCall(PetscOptionsInt(option, "What overlap should we use in construct type pardecomp?", "PCPATCH", patch->pardecomp_overlap, &patch->pardecomp_overlap, &flg));
3457: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname));
3458: PetscCall(PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg));
3459: if (flg) PetscCall(PCPatchSetSubMatType(pc, sub_mat_type));
3461: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname));
3462: PetscCall(PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg));
3464: /* If the user has set the number of subspaces, use that for the buffer size,
3465: otherwise use a large number */
3466: if (patch->nsubspaces <= 0) {
3467: nfields = 128;
3468: } else {
3469: nfields = patch->nsubspaces;
3470: }
3471: PetscCall(PetscMalloc1(nfields, &ifields));
3472: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname));
3473: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, option, ifields, &nfields, &flg));
3474: PetscCheck(!flg || !(patchConstructionType == PC_PATCH_USER), comm, PETSC_ERR_ARG_INCOMP, "We cannot support excluding a subspace with user patches because we do not index patches with a mesh point");
3475: if (flg) {
3476: PetscCall(PetscHSetIClear(patch->subspaces_to_exclude));
3477: for (k = 0; k < nfields; k++) PetscCall(PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]));
3478: }
3479: PetscCall(PetscFree(ifields));
3481: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname));
3482: PetscCall(PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg));
3483: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname));
3484: PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells));
3485: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname));
3486: PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets));
3487: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname));
3488: PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets));
3489: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname));
3490: PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints));
3491: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname));
3492: PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection));
3493: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname));
3494: PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix));
3495: PetscOptionsHeadEnd();
3496: patch->optionsSet = PETSC_TRUE;
3497: PetscFunctionReturn(PETSC_SUCCESS);
3498: }
3500: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3501: {
3502: PC_PATCH *patch = (PC_PATCH *)pc->data;
3503: KSPConvergedReason reason;
3504: PetscInt i;
3506: PetscFunctionBegin;
3507: if (!patch->save_operators) {
3508: /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3509: PetscFunctionReturn(PETSC_SUCCESS);
3510: }
3511: if (patch->denseinverse) {
3512: /* No solvers */
3513: PetscFunctionReturn(PETSC_SUCCESS);
3514: }
3515: for (i = 0; i < patch->npatch; ++i) {
3516: if (!((KSP)patch->solver[i])->setfromoptionscalled) PetscCall(KSPSetFromOptions((KSP)patch->solver[i]));
3517: PetscCall(KSPSetUp((KSP)patch->solver[i]));
3518: PetscCall(KSPGetConvergedReason((KSP)patch->solver[i], &reason));
3519: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3520: }
3521: PetscFunctionReturn(PETSC_SUCCESS);
3522: }
3524: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3525: {
3526: PC_PATCH *patch = (PC_PATCH *)pc->data;
3527: PetscViewer sviewer;
3528: PetscBool isascii;
3529: PetscMPIInt rank;
3531: PetscFunctionBegin;
3532: /* TODO Redo tabbing with set tbas in new style */
3533: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
3534: if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);
3535: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
3536: PetscCall(PetscViewerASCIIPushTab(viewer));
3537: PetscCall(PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %" PetscInt_FMT " patches\n", patch->npatch));
3538: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3539: PetscCall(PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n"));
3540: } else {
3541: PetscCall(PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n"));
3542: }
3543: if (patch->partition_of_unity) PetscCall(PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n"));
3544: else PetscCall(PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n"));
3545: if (patch->symmetrise_sweep) PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n"));
3546: else PetscCall(PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n"));
3547: if (!patch->precomputeElementTensors) PetscCall(PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n"));
3548: else PetscCall(PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n"));
3549: if (!patch->save_operators) PetscCall(PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n"));
3550: else PetscCall(PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n"));
3551: if (patch->patchconstructop == PCPatchConstruct_Star) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n"));
3552: else if (patch->patchconstructop == PCPatchConstruct_Vanka) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n"));
3553: else if (patch->patchconstructop == PCPatchConstruct_User) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n"));
3554: else PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n"));
3556: if (patch->denseinverse) {
3557: PetscCall(PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n"));
3558: } else {
3559: if (patch->isNonlinear) {
3560: PetscCall(PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n"));
3561: } else {
3562: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n"));
3563: }
3564: if (patch->solver) {
3565: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
3566: if (rank == 0) {
3567: PetscCall(PetscViewerASCIIPushTab(sviewer));
3568: PetscCall(PetscObjectView(patch->solver[0], sviewer));
3569: PetscCall(PetscViewerASCIIPopTab(sviewer));
3570: }
3571: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
3572: } else {
3573: PetscCall(PetscViewerASCIIPushTab(viewer));
3574: PetscCall(PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n"));
3575: PetscCall(PetscViewerASCIIPopTab(viewer));
3576: }
3577: }
3578: PetscCall(PetscViewerASCIIPopTab(viewer));
3579: PetscFunctionReturn(PETSC_SUCCESS);
3580: }
3582: /*MC
3583: PCPATCH - A `PC` object that encapsulates flexible definition of blocks for overlapping and non-overlapping
3584: small block additive preconditioners. Block definition is based on topology from
3585: a `DM` and equation numbering from a `PetscSection`.
3587: Options Database Keys:
3588: + -pc_patch_cells_view - Views the process local cell numbers for each patch
3589: . -pc_patch_points_view - Views the process local mesh point numbers for each patch
3590: . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch
3591: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3592: - -pc_patch_sub_mat_view - Views the matrix associated with each patch
3594: Level: intermediate
3596: .seealso: [](ch_ksp), `PCType`, `PCCreate()`, `PCSetType()`, `PCASM`, `PCJACOBI`, `PCPBJACOBI`, `PCVPBJACOBI`, `SNESPATCH`
3597: M*/
3598: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3599: {
3600: PC_PATCH *patch;
3602: PetscFunctionBegin;
3603: PetscCall(PetscCitationsRegister(PCPatchCitation, &PCPatchcite));
3604: PetscCall(PetscNew(&patch));
3606: PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));
3607: PetscCall(PetscHSetICreate(&patch->subspaces_to_exclude));
3609: patch->classname = "pc";
3610: patch->isNonlinear = PETSC_FALSE;
3612: /* Set some defaults */
3613: patch->combined = PETSC_FALSE;
3614: patch->save_operators = PETSC_TRUE;
3615: patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3616: patch->precomputeElementTensors = PETSC_FALSE;
3617: patch->partition_of_unity = PETSC_FALSE;
3618: patch->codim = -1;
3619: patch->dim = -1;
3620: patch->vankadim = -1;
3621: patch->ignoredim = -1;
3622: patch->pardecomp_overlap = 0;
3623: patch->patchconstructop = PCPatchConstruct_Star;
3624: patch->symmetrise_sweep = PETSC_FALSE;
3625: patch->npatch = 0;
3626: patch->userIS = NULL;
3627: patch->optionsSet = PETSC_FALSE;
3628: patch->iterationSet = NULL;
3629: patch->user_patches = PETSC_FALSE;
3630: PetscCall(PetscStrallocpy(MATDENSE, (char **)&patch->sub_mat_type));
3631: patch->viewPatches = PETSC_FALSE;
3632: patch->viewCells = PETSC_FALSE;
3633: patch->viewPoints = PETSC_FALSE;
3634: patch->viewSection = PETSC_FALSE;
3635: patch->viewMatrix = PETSC_FALSE;
3636: patch->densesolve = NULL;
3637: patch->setupsolver = PCSetUp_PATCH_Linear;
3638: patch->applysolver = PCApply_PATCH_Linear;
3639: patch->resetsolver = PCReset_PATCH_Linear;
3640: patch->destroysolver = PCDestroy_PATCH_Linear;
3641: patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3642: patch->dofMappingWithoutToWithArtificial = NULL;
3643: patch->dofMappingWithoutToWithAll = NULL;
3645: pc->data = (void *)patch;
3646: pc->ops->apply = PCApply_PATCH;
3647: pc->ops->applytranspose = NULL; /* PCApplyTranspose_PATCH; */
3648: pc->ops->setup = PCSetUp_PATCH;
3649: pc->ops->reset = PCReset_PATCH;
3650: pc->ops->destroy = PCDestroy_PATCH;
3651: pc->ops->setfromoptions = PCSetFromOptions_PATCH;
3652: pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH;
3653: pc->ops->view = PCView_PATCH;
3654: pc->ops->applyrichardson = NULL;
3655: PetscFunctionReturn(PETSC_SUCCESS);
3656: }