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