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