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