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: }
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /* The user's already set the patches in patch->userIS. Build the hash tables */
166: static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
167: {
168: PC_PATCH *patch = (PC_PATCH *)vpatch;
169: IS patchis = patch->userIS[point];
170: PetscInt n;
171: const PetscInt *patchdata;
172: PetscInt pStart, pEnd, i;
174: PetscFunctionBegin;
175: PetscCall(PetscHSetIClear(ht));
176: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
177: PetscCall(ISGetLocalSize(patchis, &n));
178: PetscCall(ISGetIndices(patchis, &patchdata));
179: for (i = 0; i < n; ++i) {
180: const PetscInt ownedpoint = patchdata[i];
182: 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);
183: PetscCall(PetscHSetIAdd(ht, ownedpoint));
184: }
185: PetscCall(ISRestoreIndices(patchis, &patchdata));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs)
190: {
191: PC_PATCH *patch = (PC_PATCH *)pc->data;
192: PetscInt i;
194: PetscFunctionBegin;
195: if (n == 1 && bs[0] == 1) {
196: patch->sectionSF = sf[0];
197: PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
198: } else {
199: PetscInt allRoots = 0, allLeaves = 0;
200: PetscInt leafOffset = 0;
201: PetscInt *ilocal = NULL;
202: PetscSFNode *iremote = NULL;
203: PetscInt *remoteOffsets = NULL;
204: PetscInt index = 0;
205: PetscHMapI rankToIndex;
206: PetscInt numRanks = 0;
207: PetscSFNode *remote = NULL;
208: PetscSF rankSF;
209: PetscInt *ranks = NULL;
210: PetscInt *offsets = NULL;
211: MPI_Datatype contig;
212: PetscHSetI ranksUniq;
214: /* First figure out how many dofs there are in the concatenated numbering.
215: allRoots: number of owned global dofs;
216: allLeaves: number of visible dofs (global + ghosted).
217: */
218: for (i = 0; i < n; ++i) {
219: PetscInt nroots, nleaves;
221: PetscCall(PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL));
222: allRoots += nroots * bs[i];
223: allLeaves += nleaves * bs[i];
224: }
225: PetscCall(PetscMalloc1(allLeaves, &ilocal));
226: PetscCall(PetscMalloc1(allLeaves, &iremote));
227: /* Now build an SF that just contains process connectivity. */
228: PetscCall(PetscHSetICreate(&ranksUniq));
229: for (i = 0; i < n; ++i) {
230: const PetscMPIInt *ranks = NULL;
231: PetscInt nranks, j;
233: PetscCall(PetscSFSetUp(sf[i]));
234: PetscCall(PetscSFGetRootRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL));
235: /* These are all the ranks who communicate with me. */
236: for (j = 0; j < nranks; ++j) PetscCall(PetscHSetIAdd(ranksUniq, (PetscInt)ranks[j]));
237: }
238: PetscCall(PetscHSetIGetSize(ranksUniq, &numRanks));
239: PetscCall(PetscMalloc1(numRanks, &remote));
240: PetscCall(PetscMalloc1(numRanks, &ranks));
241: PetscCall(PetscHSetIGetElems(ranksUniq, &index, ranks));
243: PetscCall(PetscHMapICreate(&rankToIndex));
244: for (i = 0; i < numRanks; ++i) {
245: remote[i].rank = ranks[i];
246: remote[i].index = 0;
247: PetscCall(PetscHMapISet(rankToIndex, ranks[i], i));
248: }
249: PetscCall(PetscFree(ranks));
250: PetscCall(PetscHSetIDestroy(&ranksUniq));
251: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &rankSF));
252: PetscCall(PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
253: PetscCall(PetscSFSetUp(rankSF));
254: /* OK, use it to communicate the root offset on the remote processes for each subspace. */
255: PetscCall(PetscMalloc1(n, &offsets));
256: PetscCall(PetscMalloc1(n * numRanks, &remoteOffsets));
258: offsets[0] = 0;
259: for (i = 1; i < n; ++i) {
260: PetscInt nroots;
262: PetscCall(PetscSFGetGraph(sf[i - 1], &nroots, NULL, NULL, NULL));
263: offsets[i] = offsets[i - 1] + nroots * bs[i - 1];
264: }
265: /* Offsets are the offsets on the current process of the global dof numbering for the subspaces. */
266: PetscCallMPI(MPI_Type_contiguous(n, 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 (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;
313: PetscFunctionBegin;
314: *dim = patch->ignoredim;
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /* TODO: Docs */
319: PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg)
320: {
321: PC_PATCH *patch = (PC_PATCH *)pc->data;
322: PetscFunctionBegin;
323: patch->save_operators = flg;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /* TODO: Docs */
328: PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg)
329: {
330: PC_PATCH *patch = (PC_PATCH *)pc->data;
331: PetscFunctionBegin;
332: *flg = patch->save_operators;
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /* TODO: Docs */
337: PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg)
338: {
339: PC_PATCH *patch = (PC_PATCH *)pc->data;
340: PetscFunctionBegin;
341: patch->precomputeElementTensors = flg;
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /* TODO: Docs */
346: PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg)
347: {
348: PC_PATCH *patch = (PC_PATCH *)pc->data;
349: PetscFunctionBegin;
350: *flg = patch->precomputeElementTensors;
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /* TODO: Docs */
355: PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg)
356: {
357: PC_PATCH *patch = (PC_PATCH *)pc->data;
358: PetscFunctionBegin;
359: patch->partition_of_unity = flg;
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: /* TODO: Docs */
364: PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg)
365: {
366: PC_PATCH *patch = (PC_PATCH *)pc->data;
367: PetscFunctionBegin;
368: *flg = patch->partition_of_unity;
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: /* TODO: Docs */
373: static PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
374: {
375: PC_PATCH *patch = (PC_PATCH *)pc->data;
376: PetscFunctionBegin;
377: PetscCheck(type == PC_COMPOSITE_ADDITIVE || type == PC_COMPOSITE_MULTIPLICATIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Only supports additive or multiplicative as the local type");
378: patch->local_composition_type = type;
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /* TODO: Docs */
383: PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type)
384: {
385: PC_PATCH *patch = (PC_PATCH *)pc->data;
387: PetscFunctionBegin;
388: if (patch->sub_mat_type) PetscCall(PetscFree(patch->sub_mat_type));
389: PetscCall(PetscStrallocpy(sub_mat_type, (char **)&patch->sub_mat_type));
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /* TODO: Docs */
394: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
395: {
396: PC_PATCH *patch = (PC_PATCH *)pc->data;
397: PetscFunctionBegin;
398: *sub_mat_type = patch->sub_mat_type;
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /* TODO: Docs */
403: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
404: {
405: PC_PATCH *patch = (PC_PATCH *)pc->data;
407: PetscFunctionBegin;
408: patch->cellNumbering = cellNumbering;
409: PetscCall(PetscObjectReference((PetscObject)cellNumbering));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /* TODO: Docs */
414: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
415: {
416: PC_PATCH *patch = (PC_PATCH *)pc->data;
417: PetscFunctionBegin;
418: *cellNumbering = patch->cellNumbering;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /* TODO: Docs */
423: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx)
424: {
425: PC_PATCH *patch = (PC_PATCH *)pc->data;
427: PetscFunctionBegin;
428: patch->ctype = ctype;
429: switch (ctype) {
430: case PC_PATCH_STAR:
431: patch->user_patches = PETSC_FALSE;
432: patch->patchconstructop = PCPatchConstruct_Star;
433: break;
434: case PC_PATCH_VANKA:
435: patch->user_patches = PETSC_FALSE;
436: patch->patchconstructop = PCPatchConstruct_Vanka;
437: break;
438: case PC_PATCH_PARDECOMP:
439: patch->user_patches = PETSC_FALSE;
440: patch->patchconstructop = PCPatchConstruct_Pardecomp;
441: break;
442: case PC_PATCH_USER:
443: case PC_PATCH_PYTHON:
444: patch->user_patches = PETSC_TRUE;
445: patch->patchconstructop = PCPatchConstruct_User;
446: if (func) {
447: patch->userpatchconstructionop = func;
448: patch->userpatchconstructctx = ctx;
449: }
450: break;
451: default:
452: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
453: }
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /* TODO: Docs */
458: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
459: {
460: PC_PATCH *patch = (PC_PATCH *)pc->data;
462: PetscFunctionBegin;
463: *ctype = patch->ctype;
464: switch (patch->ctype) {
465: case PC_PATCH_STAR:
466: case PC_PATCH_VANKA:
467: case PC_PATCH_PARDECOMP:
468: break;
469: case PC_PATCH_USER:
470: case PC_PATCH_PYTHON:
471: *func = patch->userpatchconstructionop;
472: *ctx = patch->userpatchconstructctx;
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 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)
482: {
483: PC_PATCH *patch = (PC_PATCH *)pc->data;
484: DM dm, plex;
485: PetscSF *sfs;
486: PetscInt cStart, cEnd, i, j;
488: PetscFunctionBegin;
489: PetscCall(PCGetDM(pc, &dm));
490: PetscCall(DMConvert(dm, DMPLEX, &plex));
491: dm = plex;
492: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
493: PetscCall(PetscMalloc1(nsubspaces, &sfs));
494: PetscCall(PetscMalloc1(nsubspaces, &patch->dofSection));
495: PetscCall(PetscMalloc1(nsubspaces, &patch->bs));
496: PetscCall(PetscMalloc1(nsubspaces, &patch->nodesPerCell));
497: PetscCall(PetscMalloc1(nsubspaces, &patch->cellNodeMap));
498: PetscCall(PetscMalloc1(nsubspaces + 1, &patch->subspaceOffsets));
500: patch->nsubspaces = nsubspaces;
501: patch->totalDofsPerCell = 0;
502: for (i = 0; i < nsubspaces; ++i) {
503: PetscCall(DMGetLocalSection(dms[i], &patch->dofSection[i]));
504: PetscCall(PetscObjectReference((PetscObject)patch->dofSection[i]));
505: PetscCall(DMGetSectionSF(dms[i], &sfs[i]));
506: patch->bs[i] = bs[i];
507: patch->nodesPerCell[i] = nodesPerCell[i];
508: patch->totalDofsPerCell += nodesPerCell[i] * bs[i];
509: PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
510: for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
511: patch->subspaceOffsets[i] = subspaceOffsets[i];
512: }
513: PetscCall(PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs));
514: PetscCall(PetscFree(sfs));
516: patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
517: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
518: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
519: PetscCall(DMDestroy(&dm));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: /* TODO: Docs */
524: static PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
525: {
526: PC_PATCH *patch = (PC_PATCH *)pc->data;
527: PetscInt cStart, cEnd, i, j;
529: PetscFunctionBegin;
530: patch->combined = PETSC_TRUE;
531: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
532: PetscCall(DMGetNumFields(dm, &patch->nsubspaces));
533: PetscCall(PetscCalloc1(patch->nsubspaces, &patch->dofSection));
534: PetscCall(PetscMalloc1(patch->nsubspaces, &patch->bs));
535: PetscCall(PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell));
536: PetscCall(PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap));
537: PetscCall(PetscCalloc1(patch->nsubspaces + 1, &patch->subspaceOffsets));
538: PetscCall(DMGetLocalSection(dm, &patch->dofSection[0]));
539: PetscCall(PetscObjectReference((PetscObject)patch->dofSection[0]));
540: PetscCall(PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]));
541: patch->totalDofsPerCell = 0;
542: for (i = 0; i < patch->nsubspaces; ++i) {
543: patch->bs[i] = 1;
544: patch->nodesPerCell[i] = nodesPerCell[i];
545: patch->totalDofsPerCell += nodesPerCell[i];
546: PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
547: for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
548: }
549: PetscCall(DMGetSectionSF(dm, &patch->sectionSF));
550: PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
551: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
552: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: /*@C
557: PCPatchSetComputeFunction - Set the callback function used to compute patch residuals
559: Logically Collective
561: Input Parameters:
562: + pc - The `PC`
563: . func - The callback function
564: - ctx - The user context
566: Calling sequence of `func`:
567: + pc - The `PC`
568: . point - The point
569: . x - The input solution (not used in linear problems)
570: . f - The patch residual vector
571: . cellIS - An array of the cell numbers
572: . n - The size of `dofsArray`
573: . dofsArray - The dofmap for the dofs to be solved for
574: . dofsArrayWithAll - The dofmap for all dofs on the patch
575: - ctx - The user context
577: Level: advanced
579: Note:
580: The entries of `f` (the output residual vector) have been set to zero before the call.
582: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunctionInteriorFacets()`
583: @*/
584: PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Vec f, IS cellIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
585: {
586: PC_PATCH *patch = (PC_PATCH *)pc->data;
588: PetscFunctionBegin;
589: patch->usercomputef = func;
590: patch->usercomputefctx = ctx;
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: /*@C
595: PCPatchSetComputeFunctionInteriorFacets - Set the callback function used to compute facet integrals for patch residuals
597: Logically Collective
599: Input Parameters:
600: + pc - The `PC`
601: . func - The callback function
602: - ctx - The user context
604: Calling sequence of `func`:
605: + pc - The `PC`
606: . point - The point
607: . x - The input solution (not used in linear problems)
608: . f - The patch residual vector
609: . facetIS - An array of the facet numbers
610: . n - The size of `dofsArray`
611: . dofsArray - The dofmap for the dofs to be solved for
612: . dofsArrayWithAll - The dofmap for all dofs on the patch
613: - ctx - The user context
615: Level: advanced
617: Note:
618: The entries of `f` (the output residual vector) have been set to zero before the call.
620: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunction()`
621: @*/
622: PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Vec f, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
623: {
624: PC_PATCH *patch = (PC_PATCH *)pc->data;
626: PetscFunctionBegin;
627: patch->usercomputefintfacet = func;
628: patch->usercomputefintfacetctx = ctx;
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: /*@C
633: PCPatchSetComputeOperator - Set the callback function used to compute patch matrices
635: Logically Collective
637: Input Parameters:
638: + pc - The `PC`
639: . func - The callback function
640: - ctx - The user context
642: Calling sequence of `func`:
643: + pc - The `PC`
644: . point - The point
645: . x - The input solution (not used in linear problems)
646: . mat - The patch matrix
647: . facetIS - An array of the cell numbers
648: . n - The size of `dofsArray`
649: . dofsArray - The dofmap for the dofs to be solved for
650: . dofsArrayWithAll - The dofmap for all dofs on the patch
651: - ctx - The user context
653: Level: advanced
655: Note:
656: The matrix entries have been set to zero before the call.
658: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
659: @*/
660: PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
661: {
662: PC_PATCH *patch = (PC_PATCH *)pc->data;
664: PetscFunctionBegin;
665: patch->usercomputeop = func;
666: patch->usercomputeopctx = ctx;
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: /*@C
672: PCPatchSetComputeOperatorInteriorFacets - Set the callback function used to compute facet integrals for patch matrices
674: Logically Collective
676: Input Parameters:
677: + pc - The `PC`
678: . func - The callback function
679: - ctx - The user context
681: Calling sequence of `func`:
682: + pc - The `PC`
683: . point - The point
684: . x - The input solution (not used in linear problems)
685: . mat - The patch matrix
686: . facetIS - An array of the facet numbers
687: . n - The size of `dofsArray`
688: . dofsArray - The dofmap for the dofs to be solved for
689: . dofsArrayWithAll - The dofmap for all dofs on the patch
690: - ctx - The user context
692: Level: advanced
694: Note:
695: The matrix entries have been set to zero before the call.
697: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
698: @*/
699: PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
700: {
701: PC_PATCH *patch = (PC_PATCH *)pc->data;
703: PetscFunctionBegin;
704: patch->usercomputeopintfacet = func;
705: patch->usercomputeopintfacetctx = ctx;
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: /* On entry, ht contains the topological entities whose dofs we are responsible for solving for;
710: on exit, cht contains all the topological entities we need to compute their residuals.
711: In full generality this should incorporate knowledge of the sparsity pattern of the matrix;
712: here we assume a standard FE sparsity pattern.*/
713: /* TODO: Use DMPlexGetAdjacency() */
714: static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht)
715: {
716: DM dm, plex;
717: PC_PATCH *patch = (PC_PATCH *)pc->data;
718: PetscHashIter hi;
719: PetscInt point;
720: PetscInt *star = NULL, *closure = NULL;
721: PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci;
722: PetscInt *fStar = NULL, *fClosure = NULL;
723: PetscInt fBegin, fEnd, fsi, fci, fStarSize, fClosureSize;
725: PetscFunctionBegin;
726: PetscCall(PCGetDM(pc, &dm));
727: PetscCall(DMConvert(dm, DMPLEX, &plex));
728: dm = plex;
729: PetscCall(DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd));
730: PetscCall(PCPatchGetIgnoreDim(pc, &ignoredim));
731: if (ignoredim >= 0) PetscCall(DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd));
732: PetscCall(PetscHSetIClear(cht));
733: PetscHashIterBegin(ht, hi);
734: while (!PetscHashIterAtEnd(ht, hi)) {
735: PetscHashIterGetKey(ht, hi, point);
736: PetscHashIterNext(ht, hi);
738: /* Loop over all the cells that this point connects to */
739: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
740: for (si = 0; si < starSize * 2; si += 2) {
741: const PetscInt ownedpoint = star[si];
742: /* TODO Check for point in cht before running through closure again */
743: /* now loop over all entities in the closure of that cell */
744: PetscCall(DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure));
745: for (ci = 0; ci < closureSize * 2; ci += 2) {
746: const PetscInt seenpoint = closure[ci];
747: if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue;
748: PetscCall(PetscHSetIAdd(cht, seenpoint));
749: /* Facet integrals couple dofs across facets, so in that case for each of
750: the facets we need to add all dofs on the other side of the facet to
751: the seen dofs. */
752: if (patch->usercomputeopintfacet) {
753: if (fBegin <= seenpoint && seenpoint < fEnd) {
754: PetscCall(DMPlexGetTransitiveClosure(dm, seenpoint, PETSC_FALSE, &fStarSize, &fStar));
755: for (fsi = 0; fsi < fStarSize * 2; fsi += 2) {
756: PetscCall(DMPlexGetTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, &fClosureSize, &fClosure));
757: for (fci = 0; fci < fClosureSize * 2; fci += 2) PetscCall(PetscHSetIAdd(cht, fClosure[fci]));
758: PetscCall(DMPlexRestoreTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, NULL, &fClosure));
759: }
760: PetscCall(DMPlexRestoreTransitiveClosure(dm, seenpoint, PETSC_FALSE, NULL, &fStar));
761: }
762: }
763: }
764: PetscCall(DMPlexRestoreTransitiveClosure(dm, ownedpoint, PETSC_TRUE, NULL, &closure));
765: }
766: PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, NULL, &star));
767: }
768: PetscCall(DMDestroy(&dm));
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
773: {
774: PetscFunctionBegin;
775: if (combined) {
776: if (f < 0) {
777: if (dof) PetscCall(PetscSectionGetDof(dofSection[0], p, dof));
778: if (off) PetscCall(PetscSectionGetOffset(dofSection[0], p, off));
779: } else {
780: if (dof) PetscCall(PetscSectionGetFieldDof(dofSection[0], p, f, dof));
781: if (off) PetscCall(PetscSectionGetFieldOffset(dofSection[0], p, f, off));
782: }
783: } else {
784: if (f < 0) {
785: PC_PATCH *patch = (PC_PATCH *)pc->data;
786: PetscInt fdof, g;
788: if (dof) {
789: *dof = 0;
790: for (g = 0; g < patch->nsubspaces; ++g) {
791: PetscCall(PetscSectionGetDof(dofSection[g], p, &fdof));
792: *dof += fdof;
793: }
794: }
795: if (off) {
796: *off = 0;
797: for (g = 0; g < patch->nsubspaces; ++g) {
798: PetscCall(PetscSectionGetOffset(dofSection[g], p, &fdof));
799: *off += fdof;
800: }
801: }
802: } else {
803: if (dof) PetscCall(PetscSectionGetDof(dofSection[f], p, dof));
804: if (off) PetscCall(PetscSectionGetOffset(dofSection[f], p, off));
805: }
806: }
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: /* Given a hash table with a set of topological entities (pts), compute the degrees of
811: freedom in global concatenated numbering on those entities.
812: For Vanka smoothing, this needs to do something special: ignore dofs of the
813: constraint subspace on entities that aren't the base entity we're building the patch
814: around. */
815: static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI *subspaces_to_exclude)
816: {
817: PC_PATCH *patch = (PC_PATCH *)pc->data;
818: PetscHashIter hi;
819: PetscInt ldof, loff;
820: PetscInt k, p;
822: PetscFunctionBegin;
823: PetscCall(PetscHSetIClear(dofs));
824: for (k = 0; k < patch->nsubspaces; ++k) {
825: PetscInt subspaceOffset = patch->subspaceOffsets[k];
826: PetscInt bs = patch->bs[k];
827: PetscInt j, l;
829: if (subspaces_to_exclude != NULL) {
830: PetscBool should_exclude_k = PETSC_FALSE;
831: PetscCall(PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k));
832: if (should_exclude_k) {
833: /* only get this subspace dofs at the base entity, not any others */
834: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff));
835: if (0 == ldof) continue;
836: for (j = loff; j < ldof + loff; ++j) {
837: for (l = 0; l < bs; ++l) {
838: PetscInt dof = bs * j + l + subspaceOffset;
839: PetscCall(PetscHSetIAdd(dofs, dof));
840: }
841: }
842: continue; /* skip the other dofs of this subspace */
843: }
844: }
846: PetscHashIterBegin(pts, hi);
847: while (!PetscHashIterAtEnd(pts, hi)) {
848: PetscHashIterGetKey(pts, hi, p);
849: PetscHashIterNext(pts, hi);
850: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff));
851: if (0 == ldof) continue;
852: for (j = loff; j < ldof + loff; ++j) {
853: for (l = 0; l < bs; ++l) {
854: PetscInt dof = bs * j + l + subspaceOffset;
855: PetscCall(PetscHSetIAdd(dofs, dof));
856: }
857: }
858: }
859: }
860: PetscFunctionReturn(PETSC_SUCCESS);
861: }
863: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
864: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
865: {
866: PetscHashIter hi;
867: PetscInt key;
868: PetscBool flg;
870: PetscFunctionBegin;
871: PetscCall(PetscHSetIClear(C));
872: PetscHashIterBegin(B, hi);
873: while (!PetscHashIterAtEnd(B, hi)) {
874: PetscHashIterGetKey(B, hi, key);
875: PetscHashIterNext(B, hi);
876: PetscCall(PetscHSetIHas(A, key, &flg));
877: if (!flg) PetscCall(PetscHSetIAdd(C, key));
878: }
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }
882: // PetscClangLinter pragma disable: -fdoc-sowing-chars
883: /*
884: PCPatchCreateCellPatches - create patches.
886: Input Parameter:
887: . dm - The DMPlex object defining the mesh
889: Output Parameters:
890: + cellCounts - Section with counts of cells around each vertex
891: . cells - IS of the cell point indices of cells in each patch
892: . pointCounts - Section with counts of cells around each vertex
893: - point - IS of the cell point indices of cells in each patch
894: */
895: static PetscErrorCode PCPatchCreateCellPatches(PC pc)
896: {
897: PC_PATCH *patch = (PC_PATCH *)pc->data;
898: DMLabel ghost = NULL;
899: DM dm, plex;
900: PetscHSetI ht = NULL, cht = NULL;
901: PetscSection cellCounts, pointCounts, intFacetCounts, extFacetCounts;
902: PetscInt *cellsArray, *pointsArray, *intFacetsArray, *extFacetsArray, *intFacetsToPatchCell;
903: PetscInt numCells, numPoints, numIntFacets, numExtFacets;
904: const PetscInt *leaves;
905: PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, v;
906: PetscBool isFiredrake;
908: PetscFunctionBegin;
909: /* Used to keep track of the cells in the patch. */
910: PetscCall(PetscHSetICreate(&ht));
911: PetscCall(PetscHSetICreate(&cht));
913: PetscCall(PCGetDM(pc, &dm));
914: PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC");
915: PetscCall(DMConvert(dm, DMPLEX, &plex));
916: dm = plex;
917: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
918: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
920: if (patch->user_patches) {
921: PetscCall(patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx));
922: vStart = 0;
923: vEnd = patch->npatch;
924: } else if (patch->ctype == PC_PATCH_PARDECOMP) {
925: vStart = 0;
926: vEnd = 1;
927: } else if (patch->codim < 0) {
928: if (patch->dim < 0) PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
929: else PetscCall(DMPlexGetDepthStratum(dm, patch->dim, &vStart, &vEnd));
930: } else PetscCall(DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd));
931: patch->npatch = vEnd - vStart;
933: /* These labels mark the owned points. We only create patches around points that this process owns. */
934: PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
935: if (isFiredrake) {
936: PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
937: PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
938: } else {
939: PetscSF sf;
941: PetscCall(DMGetPointSF(dm, &sf));
942: PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
943: nleaves = PetscMax(nleaves, 0);
944: }
946: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts));
947: PetscCall(PetscObjectSetName((PetscObject)patch->cellCounts, "Patch Cell Layout"));
948: cellCounts = patch->cellCounts;
949: PetscCall(PetscSectionSetChart(cellCounts, vStart, vEnd));
950: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts));
951: PetscCall(PetscObjectSetName((PetscObject)patch->pointCounts, "Patch Point Layout"));
952: pointCounts = patch->pointCounts;
953: PetscCall(PetscSectionSetChart(pointCounts, vStart, vEnd));
954: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts));
955: PetscCall(PetscObjectSetName((PetscObject)patch->extFacetCounts, "Patch Exterior Facet Layout"));
956: extFacetCounts = patch->extFacetCounts;
957: PetscCall(PetscSectionSetChart(extFacetCounts, vStart, vEnd));
958: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts));
959: PetscCall(PetscObjectSetName((PetscObject)patch->intFacetCounts, "Patch Interior Facet Layout"));
960: intFacetCounts = patch->intFacetCounts;
961: PetscCall(PetscSectionSetChart(intFacetCounts, vStart, vEnd));
962: /* Count cells and points in the patch surrounding each entity */
963: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
964: for (v = vStart; v < vEnd; ++v) {
965: PetscHashIter hi;
966: PetscInt chtSize, loc = -1;
967: PetscBool flg;
969: if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
970: if (ghost) PetscCall(DMLabelHasPoint(ghost, v, &flg));
971: else {
972: PetscCall(PetscFindInt(v, nleaves, leaves, &loc));
973: flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
974: }
975: /* Not an owned entity, don't make a cell patch. */
976: if (flg) continue;
977: }
979: PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
980: PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
981: PetscCall(PetscHSetIGetSize(cht, &chtSize));
982: /* empty patch, continue */
983: if (chtSize == 0) continue;
985: /* safe because size(cht) > 0 from above */
986: PetscHashIterBegin(cht, hi);
987: while (!PetscHashIterAtEnd(cht, hi)) {
988: PetscInt point, pdof;
990: PetscHashIterGetKey(cht, hi, point);
991: if (fStart <= point && point < fEnd) {
992: const PetscInt *support;
993: PetscInt supportSize, p;
994: PetscBool interior = PETSC_TRUE;
995: PetscCall(DMPlexGetSupport(dm, point, &support));
996: PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
997: if (supportSize == 1) {
998: interior = PETSC_FALSE;
999: } else {
1000: for (p = 0; p < supportSize; p++) {
1001: PetscBool found;
1002: /* FIXME: can I do this while iterating over cht? */
1003: PetscCall(PetscHSetIHas(cht, support[p], &found));
1004: if (!found) {
1005: interior = PETSC_FALSE;
1006: break;
1007: }
1008: }
1009: }
1010: if (interior) {
1011: PetscCall(PetscSectionAddDof(intFacetCounts, v, 1));
1012: } else {
1013: PetscCall(PetscSectionAddDof(extFacetCounts, v, 1));
1014: }
1015: }
1016: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
1017: if (pdof) PetscCall(PetscSectionAddDof(pointCounts, v, 1));
1018: if (point >= cStart && point < cEnd) PetscCall(PetscSectionAddDof(cellCounts, v, 1));
1019: PetscHashIterNext(cht, hi);
1020: }
1021: }
1022: if (isFiredrake) PetscCall(DMLabelDestroyIndex(ghost));
1024: PetscCall(PetscSectionSetUp(cellCounts));
1025: PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1026: PetscCall(PetscMalloc1(numCells, &cellsArray));
1027: PetscCall(PetscSectionSetUp(pointCounts));
1028: PetscCall(PetscSectionGetStorageSize(pointCounts, &numPoints));
1029: PetscCall(PetscMalloc1(numPoints, &pointsArray));
1031: PetscCall(PetscSectionSetUp(intFacetCounts));
1032: PetscCall(PetscSectionSetUp(extFacetCounts));
1033: PetscCall(PetscSectionGetStorageSize(intFacetCounts, &numIntFacets));
1034: PetscCall(PetscSectionGetStorageSize(extFacetCounts, &numExtFacets));
1035: PetscCall(PetscMalloc1(numIntFacets, &intFacetsArray));
1036: PetscCall(PetscMalloc1(numIntFacets * 2, &intFacetsToPatchCell));
1037: PetscCall(PetscMalloc1(numExtFacets, &extFacetsArray));
1039: /* Now that we know how much space we need, run through again and actually remember the cells. */
1040: for (v = vStart; v < vEnd; v++) {
1041: PetscHashIter hi;
1042: PetscInt dof, off, cdof, coff, efdof, efoff, ifdof, ifoff, pdof, n = 0, cn = 0, ifn = 0, efn = 0;
1044: PetscCall(PetscSectionGetDof(pointCounts, v, &dof));
1045: PetscCall(PetscSectionGetOffset(pointCounts, v, &off));
1046: PetscCall(PetscSectionGetDof(cellCounts, v, &cdof));
1047: PetscCall(PetscSectionGetOffset(cellCounts, v, &coff));
1048: PetscCall(PetscSectionGetDof(intFacetCounts, v, &ifdof));
1049: PetscCall(PetscSectionGetOffset(intFacetCounts, v, &ifoff));
1050: PetscCall(PetscSectionGetDof(extFacetCounts, v, &efdof));
1051: PetscCall(PetscSectionGetOffset(extFacetCounts, v, &efoff));
1052: if (dof <= 0) continue;
1053: PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
1054: PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
1055: PetscHashIterBegin(cht, hi);
1056: while (!PetscHashIterAtEnd(cht, hi)) {
1057: PetscInt point;
1059: PetscHashIterGetKey(cht, hi, point);
1060: if (fStart <= point && point < fEnd) {
1061: const PetscInt *support;
1062: PetscInt supportSize, p;
1063: PetscBool interior = PETSC_TRUE;
1064: PetscCall(DMPlexGetSupport(dm, point, &support));
1065: PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
1066: if (supportSize == 1) {
1067: interior = PETSC_FALSE;
1068: } else {
1069: for (p = 0; p < supportSize; p++) {
1070: PetscBool found;
1071: /* FIXME: can I do this while iterating over cht? */
1072: PetscCall(PetscHSetIHas(cht, support[p], &found));
1073: if (!found) {
1074: interior = PETSC_FALSE;
1075: break;
1076: }
1077: }
1078: }
1079: if (interior) {
1080: intFacetsToPatchCell[2 * (ifoff + ifn)] = support[0];
1081: intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = support[1];
1082: intFacetsArray[ifoff + ifn++] = point;
1083: } else {
1084: extFacetsArray[efoff + efn++] = point;
1085: }
1086: }
1087: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
1088: if (pdof) pointsArray[off + n++] = point;
1089: if (point >= cStart && point < cEnd) cellsArray[coff + cn++] = point;
1090: PetscHashIterNext(cht, hi);
1091: }
1092: 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);
1093: 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);
1094: 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);
1095: 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);
1097: for (ifn = 0; ifn < ifdof; ifn++) {
1098: PetscInt cell0 = intFacetsToPatchCell[2 * (ifoff + ifn)];
1099: PetscInt cell1 = intFacetsToPatchCell[2 * (ifoff + ifn) + 1];
1100: PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
1101: for (n = 0; n < cdof; n++) {
1102: if (!found0 && cell0 == cellsArray[coff + n]) {
1103: intFacetsToPatchCell[2 * (ifoff + ifn)] = n;
1104: found0 = PETSC_TRUE;
1105: }
1106: if (!found1 && cell1 == cellsArray[coff + n]) {
1107: intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = n;
1108: found1 = PETSC_TRUE;
1109: }
1110: if (found0 && found1) break;
1111: }
1112: PetscCheck(found0 && found1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support");
1113: }
1114: }
1115: PetscCall(PetscHSetIDestroy(&ht));
1116: PetscCall(PetscHSetIDestroy(&cht));
1118: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells));
1119: PetscCall(PetscObjectSetName((PetscObject)patch->cells, "Patch Cells"));
1120: if (patch->viewCells) {
1121: PetscCall(ObjectView((PetscObject)patch->cellCounts, patch->viewerCells, patch->formatCells));
1122: PetscCall(ObjectView((PetscObject)patch->cells, patch->viewerCells, patch->formatCells));
1123: }
1124: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets));
1125: PetscCall(PetscObjectSetName((PetscObject)patch->intFacets, "Patch Interior Facets"));
1126: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2 * numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell));
1127: PetscCall(PetscObjectSetName((PetscObject)patch->intFacetsToPatchCell, "Patch Interior Facets local support"));
1128: if (patch->viewIntFacets) {
1129: PetscCall(ObjectView((PetscObject)patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets));
1130: PetscCall(ObjectView((PetscObject)patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets));
1131: PetscCall(ObjectView((PetscObject)patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets));
1132: }
1133: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets));
1134: PetscCall(PetscObjectSetName((PetscObject)patch->extFacets, "Patch Exterior Facets"));
1135: if (patch->viewExtFacets) {
1136: PetscCall(ObjectView((PetscObject)patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets));
1137: PetscCall(ObjectView((PetscObject)patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets));
1138: }
1139: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points));
1140: PetscCall(PetscObjectSetName((PetscObject)patch->points, "Patch Points"));
1141: if (patch->viewPoints) {
1142: PetscCall(ObjectView((PetscObject)patch->pointCounts, patch->viewerPoints, patch->formatPoints));
1143: PetscCall(ObjectView((PetscObject)patch->points, patch->viewerPoints, patch->formatPoints));
1144: }
1145: PetscCall(DMDestroy(&dm));
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: /*
1150: PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches
1152: Input Parameters:
1153: + dm - The DMPlex object defining the mesh
1154: . cellCounts - Section with counts of cells around each vertex
1155: . cells - IS of the cell point indices of cells in each patch
1156: . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1157: . nodesPerCell - number of nodes per cell.
1158: - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)
1160: Output Parameters:
1161: + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering
1162: . gtolCounts - Section with counts of dofs per cell patch
1163: - gtol - IS mapping from global dofs to local dofs for each patch.
1164: */
1165: static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc)
1166: {
1167: PC_PATCH *patch = (PC_PATCH *)pc->data;
1168: PetscSection cellCounts = patch->cellCounts;
1169: PetscSection pointCounts = patch->pointCounts;
1170: PetscSection gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL;
1171: IS cells = patch->cells;
1172: IS points = patch->points;
1173: PetscSection cellNumbering = patch->cellNumbering;
1174: PetscInt Nf = patch->nsubspaces;
1175: PetscInt numCells, numPoints;
1176: PetscInt numDofs;
1177: PetscInt numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll;
1178: PetscInt totalDofsPerCell = patch->totalDofsPerCell;
1179: PetscInt vStart, vEnd, v;
1180: const PetscInt *cellsArray, *pointsArray;
1181: PetscInt *newCellsArray = NULL;
1182: PetscInt *dofsArray = NULL;
1183: PetscInt *dofsArrayWithArtificial = NULL;
1184: PetscInt *dofsArrayWithAll = NULL;
1185: PetscInt *offsArray = NULL;
1186: PetscInt *offsArrayWithArtificial = NULL;
1187: PetscInt *offsArrayWithAll = NULL;
1188: PetscInt *asmArray = NULL;
1189: PetscInt *asmArrayWithArtificial = NULL;
1190: PetscInt *asmArrayWithAll = NULL;
1191: PetscInt *globalDofsArray = NULL;
1192: PetscInt *globalDofsArrayWithArtificial = NULL;
1193: PetscInt *globalDofsArrayWithAll = NULL;
1194: PetscInt globalIndex = 0;
1195: PetscInt key = 0;
1196: PetscInt asmKey = 0;
1197: DM dm = NULL, plex;
1198: const PetscInt *bcNodes = NULL;
1199: PetscHMapI ht;
1200: PetscHMapI htWithArtificial;
1201: PetscHMapI htWithAll;
1202: PetscHSetI globalBcs;
1203: PetscInt numBcs;
1204: PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs;
1205: PetscInt pStart, pEnd, p, i;
1206: char option[PETSC_MAX_PATH_LEN];
1207: PetscBool isNonlinear;
1209: PetscFunctionBegin;
1211: PetscCall(PCGetDM(pc, &dm));
1212: PetscCall(DMConvert(dm, DMPLEX, &plex));
1213: dm = plex;
1214: /* dofcounts section is cellcounts section * dofPerCell */
1215: PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1216: PetscCall(PetscSectionGetStorageSize(patch->pointCounts, &numPoints));
1217: numDofs = numCells * totalDofsPerCell;
1218: PetscCall(PetscMalloc1(numDofs, &dofsArray));
1219: PetscCall(PetscMalloc1(numPoints * Nf, &offsArray));
1220: PetscCall(PetscMalloc1(numDofs, &asmArray));
1221: PetscCall(PetscMalloc1(numCells, &newCellsArray));
1222: PetscCall(PetscSectionGetChart(cellCounts, &vStart, &vEnd));
1223: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts));
1224: gtolCounts = patch->gtolCounts;
1225: PetscCall(PetscSectionSetChart(gtolCounts, vStart, vEnd));
1226: PetscCall(PetscObjectSetName((PetscObject)patch->gtolCounts, "Patch Global Index Section"));
1228: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1229: PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithArtificial));
1230: PetscCall(PetscMalloc1(numDofs, &asmArrayWithArtificial));
1231: PetscCall(PetscMalloc1(numDofs, &dofsArrayWithArtificial));
1232: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial));
1233: gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1234: PetscCall(PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd));
1235: PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs"));
1236: }
1238: isNonlinear = patch->isNonlinear;
1239: if (isNonlinear) {
1240: PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithAll));
1241: PetscCall(PetscMalloc1(numDofs, &asmArrayWithAll));
1242: PetscCall(PetscMalloc1(numDofs, &dofsArrayWithAll));
1243: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll));
1244: gtolCountsWithAll = patch->gtolCountsWithAll;
1245: PetscCall(PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd));
1246: PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs"));
1247: }
1249: /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1250: conditions */
1251: PetscCall(PetscHSetICreate(&globalBcs));
1252: PetscCall(ISGetIndices(patch->ghostBcNodes, &bcNodes));
1253: PetscCall(ISGetSize(patch->ghostBcNodes, &numBcs));
1254: for (i = 0; i < numBcs; ++i) { PetscCall(PetscHSetIAdd(globalBcs, bcNodes[i])); /* these are already in concatenated numbering */ }
1255: PetscCall(ISRestoreIndices(patch->ghostBcNodes, &bcNodes));
1256: PetscCall(ISDestroy(&patch->ghostBcNodes)); /* memory optimisation */
1258: /* Hash tables for artificial BC construction */
1259: PetscCall(PetscHSetICreate(&ownedpts));
1260: PetscCall(PetscHSetICreate(&seenpts));
1261: PetscCall(PetscHSetICreate(&owneddofs));
1262: PetscCall(PetscHSetICreate(&seendofs));
1263: PetscCall(PetscHSetICreate(&artificialbcs));
1265: PetscCall(ISGetIndices(cells, &cellsArray));
1266: PetscCall(ISGetIndices(points, &pointsArray));
1267: PetscCall(PetscHMapICreate(&ht));
1268: PetscCall(PetscHMapICreate(&htWithArtificial));
1269: PetscCall(PetscHMapICreate(&htWithAll));
1270: for (v = vStart; v < vEnd; ++v) {
1271: PetscInt localIndex = 0;
1272: PetscInt localIndexWithArtificial = 0;
1273: PetscInt localIndexWithAll = 0;
1274: PetscInt dof, off, i, j, k, l;
1276: PetscCall(PetscHMapIClear(ht));
1277: PetscCall(PetscHMapIClear(htWithArtificial));
1278: PetscCall(PetscHMapIClear(htWithAll));
1279: PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1280: PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1281: if (dof <= 0) continue;
1283: /* Calculate the global numbers of the artificial BC dofs here first */
1284: PetscCall(patch->patchconstructop((void *)patch, dm, v, ownedpts));
1285: PetscCall(PCPatchCompleteCellPatch(pc, ownedpts, seenpts));
1286: PetscCall(PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude));
1287: PetscCall(PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL));
1288: PetscCall(PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs));
1289: if (patch->viewPatches) {
1290: PetscHSetI globalbcdofs;
1291: PetscHashIter hi;
1292: MPI_Comm comm = PetscObjectComm((PetscObject)pc);
1294: PetscCall(PetscHSetICreate(&globalbcdofs));
1295: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": owned dofs:\n", v));
1296: PetscHashIterBegin(owneddofs, hi);
1297: while (!PetscHashIterAtEnd(owneddofs, hi)) {
1298: PetscInt globalDof;
1300: PetscHashIterGetKey(owneddofs, hi, globalDof);
1301: PetscHashIterNext(owneddofs, hi);
1302: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1303: }
1304: PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1305: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": seen dofs:\n", v));
1306: PetscHashIterBegin(seendofs, hi);
1307: while (!PetscHashIterAtEnd(seendofs, hi)) {
1308: PetscInt globalDof;
1309: PetscBool flg;
1311: PetscHashIterGetKey(seendofs, hi, globalDof);
1312: PetscHashIterNext(seendofs, hi);
1313: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1315: PetscCall(PetscHSetIHas(globalBcs, globalDof, &flg));
1316: if (flg) PetscCall(PetscHSetIAdd(globalbcdofs, globalDof));
1317: }
1318: PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1319: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": global BCs:\n", v));
1320: PetscCall(PetscHSetIGetSize(globalbcdofs, &numBcs));
1321: if (numBcs > 0) {
1322: PetscHashIterBegin(globalbcdofs, hi);
1323: while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1324: PetscInt globalDof;
1325: PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1326: PetscHashIterNext(globalbcdofs, hi);
1327: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1328: }
1329: }
1330: PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1331: PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": artificial BCs:\n", v));
1332: PetscCall(PetscHSetIGetSize(artificialbcs, &numBcs));
1333: if (numBcs > 0) {
1334: PetscHashIterBegin(artificialbcs, hi);
1335: while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1336: PetscInt globalDof;
1337: PetscHashIterGetKey(artificialbcs, hi, globalDof);
1338: PetscHashIterNext(artificialbcs, hi);
1339: PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1340: }
1341: }
1342: PetscCall(PetscSynchronizedPrintf(comm, "\n\n"));
1343: PetscCall(PetscHSetIDestroy(&globalbcdofs));
1344: }
1345: for (k = 0; k < patch->nsubspaces; ++k) {
1346: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1347: PetscInt nodesPerCell = patch->nodesPerCell[k];
1348: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1349: PetscInt bs = patch->bs[k];
1351: for (i = off; i < off + dof; ++i) {
1352: /* Walk over the cells in this patch. */
1353: const PetscInt c = cellsArray[i];
1354: PetscInt cell = c;
1356: /* TODO Change this to an IS */
1357: if (cellNumbering) {
1358: PetscCall(PetscSectionGetDof(cellNumbering, c, &cell));
1359: PetscCheck(cell > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " doesn't appear in cell numbering map", c);
1360: PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1361: }
1362: newCellsArray[i] = cell;
1363: for (j = 0; j < nodesPerCell; ++j) {
1364: /* For each global dof, map it into contiguous local storage. */
1365: const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + subspaceOffset;
1366: /* finally, loop over block size */
1367: for (l = 0; l < bs; ++l) {
1368: PetscInt localDof;
1369: PetscBool isGlobalBcDof, isArtificialBcDof;
1371: /* first, check if this is either a globally enforced or locally enforced BC dof */
1372: PetscCall(PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof));
1373: PetscCall(PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof));
1375: /* if it's either, don't ever give it a local dof number */
1376: if (isGlobalBcDof || isArtificialBcDof) {
1377: dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1378: } else {
1379: PetscCall(PetscHMapIGet(ht, globalDof + l, &localDof));
1380: if (localDof == -1) {
1381: localDof = localIndex++;
1382: PetscCall(PetscHMapISet(ht, globalDof + l, localDof));
1383: }
1384: PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1385: /* And store. */
1386: dofsArray[globalIndex] = localDof;
1387: }
1389: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1390: if (isGlobalBcDof) {
1391: dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1392: } else {
1393: PetscCall(PetscHMapIGet(htWithArtificial, globalDof + l, &localDof));
1394: if (localDof == -1) {
1395: localDof = localIndexWithArtificial++;
1396: PetscCall(PetscHMapISet(htWithArtificial, globalDof + l, localDof));
1397: }
1398: PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1399: /* And store.*/
1400: dofsArrayWithArtificial[globalIndex] = localDof;
1401: }
1402: }
1404: if (isNonlinear) {
1405: /* Build the dofmap for the function space with _all_ dofs,
1406: including those in any kind of boundary condition */
1407: PetscCall(PetscHMapIGet(htWithAll, globalDof + l, &localDof));
1408: if (localDof == -1) {
1409: localDof = localIndexWithAll++;
1410: PetscCall(PetscHMapISet(htWithAll, globalDof + l, localDof));
1411: }
1412: PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1413: /* And store.*/
1414: dofsArrayWithAll[globalIndex] = localDof;
1415: }
1416: globalIndex++;
1417: }
1418: }
1419: }
1420: }
1421: /*How many local dofs in this patch? */
1422: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1423: PetscCall(PetscHMapIGetSize(htWithArtificial, &dof));
1424: PetscCall(PetscSectionSetDof(gtolCountsWithArtificial, v, dof));
1425: }
1426: if (isNonlinear) {
1427: PetscCall(PetscHMapIGetSize(htWithAll, &dof));
1428: PetscCall(PetscSectionSetDof(gtolCountsWithAll, v, dof));
1429: }
1430: PetscCall(PetscHMapIGetSize(ht, &dof));
1431: PetscCall(PetscSectionSetDof(gtolCounts, v, dof));
1432: }
1434: PetscCall(DMDestroy(&dm));
1435: 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);
1436: PetscCall(PetscSectionSetUp(gtolCounts));
1437: PetscCall(PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs));
1438: PetscCall(PetscMalloc1(numGlobalDofs, &globalDofsArray));
1440: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1441: PetscCall(PetscSectionSetUp(gtolCountsWithArtificial));
1442: PetscCall(PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial));
1443: PetscCall(PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial));
1444: }
1445: if (isNonlinear) {
1446: PetscCall(PetscSectionSetUp(gtolCountsWithAll));
1447: PetscCall(PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll));
1448: PetscCall(PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll));
1449: }
1450: /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */
1451: for (v = vStart; v < vEnd; ++v) {
1452: PetscHashIter hi;
1453: PetscInt dof, off, Np, ooff, i, j, k, l;
1455: PetscCall(PetscHMapIClear(ht));
1456: PetscCall(PetscHMapIClear(htWithArtificial));
1457: PetscCall(PetscHMapIClear(htWithAll));
1458: PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1459: PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1460: PetscCall(PetscSectionGetDof(pointCounts, v, &Np));
1461: PetscCall(PetscSectionGetOffset(pointCounts, v, &ooff));
1462: if (dof <= 0) continue;
1464: for (k = 0; k < patch->nsubspaces; ++k) {
1465: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1466: PetscInt nodesPerCell = patch->nodesPerCell[k];
1467: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1468: PetscInt bs = patch->bs[k];
1469: PetscInt goff;
1471: for (i = off; i < off + dof; ++i) {
1472: /* Reconstruct mapping of global-to-local on this patch. */
1473: const PetscInt c = cellsArray[i];
1474: PetscInt cell = c;
1476: if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1477: for (j = 0; j < nodesPerCell; ++j) {
1478: for (l = 0; l < bs; ++l) {
1479: const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1480: const PetscInt localDof = dofsArray[key];
1481: if (localDof >= 0) PetscCall(PetscHMapISet(ht, globalDof, localDof));
1482: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1483: const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1484: if (localDofWithArtificial >= 0) PetscCall(PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial));
1485: }
1486: if (isNonlinear) {
1487: const PetscInt localDofWithAll = dofsArrayWithAll[key];
1488: if (localDofWithAll >= 0) PetscCall(PetscHMapISet(htWithAll, globalDof, localDofWithAll));
1489: }
1490: key++;
1491: }
1492: }
1493: }
1495: /* Shove it in the output data structure. */
1496: PetscCall(PetscSectionGetOffset(gtolCounts, v, &goff));
1497: PetscHashIterBegin(ht, hi);
1498: while (!PetscHashIterAtEnd(ht, hi)) {
1499: PetscInt globalDof, localDof;
1501: PetscHashIterGetKey(ht, hi, globalDof);
1502: PetscHashIterGetVal(ht, hi, localDof);
1503: if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1504: PetscHashIterNext(ht, hi);
1505: }
1507: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1508: PetscCall(PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff));
1509: PetscHashIterBegin(htWithArtificial, hi);
1510: while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1511: PetscInt globalDof, localDof;
1512: PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1513: PetscHashIterGetVal(htWithArtificial, hi, localDof);
1514: if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1515: PetscHashIterNext(htWithArtificial, hi);
1516: }
1517: }
1518: if (isNonlinear) {
1519: PetscCall(PetscSectionGetOffset(gtolCountsWithAll, v, &goff));
1520: PetscHashIterBegin(htWithAll, hi);
1521: while (!PetscHashIterAtEnd(htWithAll, hi)) {
1522: PetscInt globalDof, localDof;
1523: PetscHashIterGetKey(htWithAll, hi, globalDof);
1524: PetscHashIterGetVal(htWithAll, hi, localDof);
1525: if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1526: PetscHashIterNext(htWithAll, hi);
1527: }
1528: }
1530: for (p = 0; p < Np; ++p) {
1531: const PetscInt point = pointsArray[ooff + p];
1532: PetscInt globalDof, localDof;
1534: PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof));
1535: PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
1536: offsArray[(ooff + p) * Nf + k] = localDof;
1537: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1538: PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
1539: offsArrayWithArtificial[(ooff + p) * Nf + k] = localDof;
1540: }
1541: if (isNonlinear) {
1542: PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
1543: offsArrayWithAll[(ooff + p) * Nf + k] = localDof;
1544: }
1545: }
1546: }
1548: PetscCall(PetscHSetIDestroy(&globalBcs));
1549: PetscCall(PetscHSetIDestroy(&ownedpts));
1550: PetscCall(PetscHSetIDestroy(&seenpts));
1551: PetscCall(PetscHSetIDestroy(&owneddofs));
1552: PetscCall(PetscHSetIDestroy(&seendofs));
1553: PetscCall(PetscHSetIDestroy(&artificialbcs));
1555: /* At this point, we have a hash table ht built that maps globalDof -> localDof.
1556: We need to create the dof table laid out cellwise first, then by subspace,
1557: as the assembler assembles cell-wise and we need to stuff the different
1558: contributions of the different function spaces to the right places. So we loop
1559: over cells, then over subspaces. */
1560: if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */
1561: for (i = off; i < off + dof; ++i) {
1562: const PetscInt c = cellsArray[i];
1563: PetscInt cell = c;
1565: if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1566: for (k = 0; k < patch->nsubspaces; ++k) {
1567: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1568: PetscInt nodesPerCell = patch->nodesPerCell[k];
1569: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1570: PetscInt bs = patch->bs[k];
1572: for (j = 0; j < nodesPerCell; ++j) {
1573: for (l = 0; l < bs; ++l) {
1574: const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1575: PetscInt localDof;
1577: PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
1578: /* If it's not in the hash table, i.e. is a BC dof,
1579: then the PetscHSetIMap above gives -1, which matches
1580: exactly the convention for PETSc's matrix assembly to
1581: ignore the dof. So we don't need to do anything here */
1582: asmArray[asmKey] = localDof;
1583: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1584: PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
1585: asmArrayWithArtificial[asmKey] = localDof;
1586: }
1587: if (isNonlinear) {
1588: PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
1589: asmArrayWithAll[asmKey] = localDof;
1590: }
1591: asmKey++;
1592: }
1593: }
1594: }
1595: }
1596: }
1597: }
1598: if (1 == patch->nsubspaces) {
1599: PetscCall(PetscArraycpy(asmArray, dofsArray, numDofs));
1600: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs));
1601: if (isNonlinear) PetscCall(PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs));
1602: }
1604: PetscCall(PetscHMapIDestroy(&ht));
1605: PetscCall(PetscHMapIDestroy(&htWithArtificial));
1606: PetscCall(PetscHMapIDestroy(&htWithAll));
1607: PetscCall(ISRestoreIndices(cells, &cellsArray));
1608: PetscCall(ISRestoreIndices(points, &pointsArray));
1609: PetscCall(PetscFree(dofsArray));
1610: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscFree(dofsArrayWithArtificial));
1611: if (isNonlinear) PetscCall(PetscFree(dofsArrayWithAll));
1612: /* Create placeholder section for map from points to patch dofs */
1613: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection));
1614: PetscCall(PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces));
1615: if (patch->combined) {
1616: PetscInt numFields;
1617: PetscCall(PetscSectionGetNumFields(patch->dofSection[0], &numFields));
1618: 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);
1619: PetscCall(PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd));
1620: PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
1621: for (p = pStart; p < pEnd; ++p) {
1622: PetscInt dof, fdof, f;
1624: PetscCall(PetscSectionGetDof(patch->dofSection[0], p, &dof));
1625: PetscCall(PetscSectionSetDof(patch->patchSection, p, dof));
1626: for (f = 0; f < patch->nsubspaces; ++f) {
1627: PetscCall(PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof));
1628: PetscCall(PetscSectionSetFieldDof(patch->patchSection, p, f, fdof));
1629: }
1630: }
1631: } else {
1632: PetscInt pStartf, pEndf, f;
1633: pStart = PETSC_MAX_INT;
1634: pEnd = PETSC_MIN_INT;
1635: for (f = 0; f < patch->nsubspaces; ++f) {
1636: PetscCall(PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf));
1637: pStart = PetscMin(pStart, pStartf);
1638: pEnd = PetscMax(pEnd, pEndf);
1639: }
1640: PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
1641: for (f = 0; f < patch->nsubspaces; ++f) {
1642: PetscCall(PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf));
1643: for (p = pStartf; p < pEndf; ++p) {
1644: PetscInt fdof;
1645: PetscCall(PetscSectionGetDof(patch->dofSection[f], p, &fdof));
1646: PetscCall(PetscSectionAddDof(patch->patchSection, p, fdof));
1647: PetscCall(PetscSectionSetFieldDof(patch->patchSection, p, f, fdof));
1648: }
1649: }
1650: }
1651: PetscCall(PetscSectionSetUp(patch->patchSection));
1652: PetscCall(PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE));
1653: /* Replace cell indices with firedrake-numbered ones. */
1654: PetscCall(ISGeneralSetIndices(cells, numCells, (const PetscInt *)newCellsArray, PETSC_OWN_POINTER));
1655: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol));
1656: PetscCall(PetscObjectSetName((PetscObject)patch->gtol, "Global Indices"));
1657: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname));
1658: PetscCall(PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject)pc, option));
1659: PetscCall(ISViewFromOptions(patch->gtol, (PetscObject)pc, option));
1660: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs));
1661: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArray, PETSC_OWN_POINTER, &patch->offs));
1662: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1663: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial));
1664: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial));
1665: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial));
1666: }
1667: if (isNonlinear) {
1668: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll));
1669: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll));
1670: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll));
1671: }
1672: PetscFunctionReturn(PETSC_SUCCESS);
1673: }
1675: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1676: {
1677: PC_PATCH *patch = (PC_PATCH *)pc->data;
1678: PetscBool flg;
1679: PetscInt csize, rsize;
1680: const char *prefix = NULL;
1682: PetscFunctionBegin;
1683: if (withArtificial) {
1684: /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1685: PetscInt pStart;
1686: PetscCall(PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL));
1687: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize));
1688: csize = rsize;
1689: } else {
1690: PetscInt pStart;
1691: PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
1692: PetscCall(PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize));
1693: csize = rsize;
1694: }
1696: PetscCall(MatCreate(PETSC_COMM_SELF, mat));
1697: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1698: PetscCall(MatSetOptionsPrefix(*mat, prefix));
1699: PetscCall(MatAppendOptionsPrefix(*mat, "pc_patch_sub_"));
1700: if (patch->sub_mat_type) PetscCall(MatSetType(*mat, patch->sub_mat_type));
1701: else if (!patch->sub_mat_type) PetscCall(MatSetType(*mat, MATDENSE));
1702: PetscCall(MatSetSizes(*mat, rsize, csize, rsize, csize));
1703: PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATDENSE, &flg));
1704: if (!flg) PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg));
1705: /* Sparse patch matrices */
1706: if (!flg) {
1707: PetscBT bt;
1708: PetscInt *dnnz = NULL;
1709: const PetscInt *dofsArray = NULL;
1710: PetscInt pStart, pEnd, ncell, offset, c, i, j;
1712: if (withArtificial) {
1713: PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
1714: } else {
1715: PetscCall(ISGetIndices(patch->dofs, &dofsArray));
1716: }
1717: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
1718: point += pStart;
1719: PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
1720: PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
1721: PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
1722: PetscCall(PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0));
1723: /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1724: * patch. This is probably OK if the patches are not too big,
1725: * but uses too much memory. We therefore switch based on rsize. */
1726: if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1727: PetscScalar *zeroes;
1728: PetscInt rows;
1730: PetscCall(PetscCalloc1(rsize, &dnnz));
1731: PetscCall(PetscBTCreate(rsize * rsize, &bt));
1732: for (c = 0; c < ncell; ++c) {
1733: const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1734: for (i = 0; i < patch->totalDofsPerCell; ++i) {
1735: const PetscInt row = idx[i];
1736: if (row < 0) continue;
1737: for (j = 0; j < patch->totalDofsPerCell; ++j) {
1738: const PetscInt col = idx[j];
1739: const PetscInt key = row * rsize + col;
1740: if (col < 0) continue;
1741: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1742: }
1743: }
1744: }
1746: if (patch->usercomputeopintfacet) {
1747: const PetscInt *intFacetsArray = NULL;
1748: PetscInt i, numIntFacets, intFacetOffset;
1749: const PetscInt *facetCells = NULL;
1751: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1752: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1753: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1754: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1755: for (i = 0; i < numIntFacets; i++) {
1756: const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
1757: const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
1758: PetscInt celli, cellj;
1760: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1761: const PetscInt row = dofsArray[(offset + cell0) * patch->totalDofsPerCell + celli];
1762: if (row < 0) continue;
1763: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1764: const PetscInt col = dofsArray[(offset + cell1) * patch->totalDofsPerCell + cellj];
1765: const PetscInt key = row * rsize + col;
1766: if (col < 0) continue;
1767: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1768: }
1769: }
1771: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1772: const PetscInt row = dofsArray[(offset + cell1) * patch->totalDofsPerCell + celli];
1773: if (row < 0) continue;
1774: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1775: const PetscInt col = dofsArray[(offset + cell0) * patch->totalDofsPerCell + cellj];
1776: const PetscInt key = row * rsize + col;
1777: if (col < 0) continue;
1778: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1779: }
1780: }
1781: }
1782: }
1783: PetscCall(PetscBTDestroy(&bt));
1784: PetscCall(MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL));
1785: PetscCall(PetscFree(dnnz));
1787: PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &zeroes));
1788: for (c = 0; c < ncell; ++c) {
1789: const PetscInt *idx = &dofsArray[(offset + c) * patch->totalDofsPerCell];
1790: PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES));
1791: }
1792: PetscCall(MatGetLocalSize(*mat, &rows, NULL));
1793: for (i = 0; i < rows; ++i) PetscCall(MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES));
1795: if (patch->usercomputeopintfacet) {
1796: const PetscInt *intFacetsArray = NULL;
1797: PetscInt i, numIntFacets, intFacetOffset;
1798: const PetscInt *facetCells = NULL;
1800: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1801: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1802: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1803: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1804: for (i = 0; i < numIntFacets; i++) {
1805: const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
1806: const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
1807: const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1808: const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1809: PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES));
1810: PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES));
1811: }
1812: }
1814: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
1815: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
1817: PetscCall(PetscFree(zeroes));
1819: } else { /* rsize too big, use MATPREALLOCATOR */
1820: Mat preallocator;
1821: PetscScalar *vals;
1823: PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &vals));
1824: PetscCall(MatCreate(PETSC_COMM_SELF, &preallocator));
1825: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
1826: PetscCall(MatSetSizes(preallocator, rsize, rsize, rsize, rsize));
1827: PetscCall(MatSetUp(preallocator));
1829: for (c = 0; c < ncell; ++c) {
1830: const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1831: PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES));
1832: }
1834: if (patch->usercomputeopintfacet) {
1835: const PetscInt *intFacetsArray = NULL;
1836: PetscInt i, numIntFacets, intFacetOffset;
1837: const PetscInt *facetCells = NULL;
1839: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1840: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1841: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1842: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1843: for (i = 0; i < numIntFacets; i++) {
1844: const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
1845: const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
1846: const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1847: const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1848: PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES));
1849: PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES));
1850: }
1851: }
1853: PetscCall(PetscFree(vals));
1854: PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
1855: PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
1856: PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat));
1857: PetscCall(MatDestroy(&preallocator));
1858: }
1859: PetscCall(PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0));
1860: if (withArtificial) {
1861: PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
1862: } else {
1863: PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
1864: }
1865: }
1866: PetscCall(MatSetUp(*mat));
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: static PetscErrorCode PCPatchComputeFunction_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Vec F, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx)
1871: {
1872: PC_PATCH *patch = (PC_PATCH *)pc->data;
1873: DM dm, plex;
1874: PetscSection s;
1875: const PetscInt *parray, *oarray;
1876: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
1878: PetscFunctionBegin;
1879: PetscCheck(!patch->precomputeElementTensors, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
1880: PetscCall(PCGetDM(pc, &dm));
1881: PetscCall(DMConvert(dm, DMPLEX, &plex));
1882: dm = plex;
1883: PetscCall(DMGetLocalSection(dm, &s));
1884: /* Set offset into patch */
1885: PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
1886: PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
1887: PetscCall(ISGetIndices(patch->points, &parray));
1888: PetscCall(ISGetIndices(patch->offs, &oarray));
1889: for (f = 0; f < Nf; ++f) {
1890: for (p = 0; p < Np; ++p) {
1891: const PetscInt point = parray[poff + p];
1892: PetscInt dof;
1894: PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
1895: PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
1896: if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
1897: else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
1898: }
1899: }
1900: PetscCall(ISRestoreIndices(patch->points, &parray));
1901: PetscCall(ISRestoreIndices(patch->offs, &oarray));
1902: if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
1903: PetscCall(DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx));
1904: PetscCall(DMDestroy(&dm));
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1908: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1909: {
1910: PC_PATCH *patch = (PC_PATCH *)pc->data;
1911: const PetscInt *dofsArray;
1912: const PetscInt *dofsArrayWithAll;
1913: const PetscInt *cellsArray;
1914: PetscInt ncell, offset, pStart, pEnd;
1916: PetscFunctionBegin;
1917: PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
1918: PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
1919: PetscCall(ISGetIndices(patch->dofs, &dofsArray));
1920: PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
1921: PetscCall(ISGetIndices(patch->cells, &cellsArray));
1922: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
1924: point += pStart;
1925: PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
1927: PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
1928: PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
1929: if (ncell <= 0) {
1930: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
1931: PetscFunctionReturn(PETSC_SUCCESS);
1932: }
1933: PetscCall(VecSet(F, 0.0));
1934: /* Cannot reuse the same IS because the geometry info is being cached in it */
1935: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS));
1936: PetscCallBack("PCPatch callback", patch->usercomputef(pc, point, x, F, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, dofsArrayWithAll + offset * patch->totalDofsPerCell, patch->usercomputefctx));
1937: PetscCall(ISDestroy(&patch->cellIS));
1938: PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
1939: PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
1940: PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
1941: if (patch->viewMatrix) {
1942: char name[PETSC_MAX_PATH_LEN];
1944: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch vector for Point %" PetscInt_FMT, point));
1945: PetscCall(PetscObjectSetName((PetscObject)F, name));
1946: PetscCall(ObjectView((PetscObject)F, patch->viewerMatrix, patch->formatMatrix));
1947: }
1948: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
1949: PetscFunctionReturn(PETSC_SUCCESS);
1950: }
1952: static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx)
1953: {
1954: PC_PATCH *patch = (PC_PATCH *)pc->data;
1955: DM dm, plex;
1956: PetscSection s;
1957: const PetscInt *parray, *oarray;
1958: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
1960: PetscFunctionBegin;
1961: PetscCall(PCGetDM(pc, &dm));
1962: PetscCall(DMConvert(dm, DMPLEX, &plex));
1963: dm = plex;
1964: PetscCall(DMGetLocalSection(dm, &s));
1965: /* Set offset into patch */
1966: PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
1967: PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
1968: PetscCall(ISGetIndices(patch->points, &parray));
1969: PetscCall(ISGetIndices(patch->offs, &oarray));
1970: for (f = 0; f < Nf; ++f) {
1971: for (p = 0; p < Np; ++p) {
1972: const PetscInt point = parray[poff + p];
1973: PetscInt dof;
1975: PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
1976: PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
1977: if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
1978: else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
1979: }
1980: }
1981: PetscCall(ISRestoreIndices(patch->points, &parray));
1982: PetscCall(ISRestoreIndices(patch->offs, &oarray));
1983: if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
1984: /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
1985: PetscCall(DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx));
1986: PetscCall(DMDestroy(&dm));
1987: PetscFunctionReturn(PETSC_SUCCESS);
1988: }
1990: /* This function zeros mat on entry */
1991: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
1992: {
1993: PC_PATCH *patch = (PC_PATCH *)pc->data;
1994: const PetscInt *dofsArray;
1995: const PetscInt *dofsArrayWithAll = NULL;
1996: const PetscInt *cellsArray;
1997: PetscInt ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
1998: PetscBool isNonlinear;
2000: PetscFunctionBegin;
2001: PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
2002: isNonlinear = patch->isNonlinear;
2003: PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
2004: if (withArtificial) {
2005: PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
2006: } else {
2007: PetscCall(ISGetIndices(patch->dofs, &dofsArray));
2008: }
2009: if (isNonlinear) PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
2010: PetscCall(ISGetIndices(patch->cells, &cellsArray));
2011: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
2013: point += pStart;
2014: PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
2016: PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
2017: PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
2018: if (ncell <= 0) {
2019: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2020: PetscFunctionReturn(PETSC_SUCCESS);
2021: }
2022: PetscCall(MatZeroEntries(mat));
2023: if (patch->precomputeElementTensors) {
2024: PetscInt i;
2025: PetscInt ndof = patch->totalDofsPerCell;
2026: const PetscScalar *elementTensors;
2028: PetscCall(VecGetArrayRead(patch->cellMats, &elementTensors));
2029: for (i = 0; i < ncell; i++) {
2030: const PetscInt cell = cellsArray[i + offset];
2031: const PetscInt *idx = dofsArray + (offset + i) * ndof;
2032: const PetscScalar *v = elementTensors + patch->precomputedTensorLocations[cell] * ndof * ndof;
2033: PetscCall(MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES));
2034: }
2035: PetscCall(VecRestoreArrayRead(patch->cellMats, &elementTensors));
2036: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2037: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2038: } else {
2039: /* Cannot reuse the same IS because the geometry info is being cached in it */
2040: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS));
2041: PetscCallBack("PCPatch callback",
2042: patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset * patch->totalDofsPerCell : NULL, patch->usercomputeopctx));
2043: }
2044: if (patch->usercomputeopintfacet) {
2045: PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
2046: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
2047: if (numIntFacets > 0) {
2048: /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2049: PetscInt *facetDofs = NULL, *facetDofsWithAll = NULL;
2050: const PetscInt *intFacetsArray = NULL;
2051: PetscInt idx = 0;
2052: PetscInt i, c, d;
2053: PetscInt fStart;
2054: DM dm, plex;
2055: IS facetIS = NULL;
2056: const PetscInt *facetCells = NULL;
2058: PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
2059: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2060: PetscCall(PCGetDM(pc, &dm));
2061: PetscCall(DMConvert(dm, DMPLEX, &plex));
2062: dm = plex;
2063: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, NULL));
2064: /* FIXME: Pull this malloc out. */
2065: PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs));
2066: if (dofsArrayWithAll) PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll));
2067: if (patch->precomputeElementTensors) {
2068: PetscInt nFacetDof = 2 * patch->totalDofsPerCell;
2069: const PetscScalar *elementTensors;
2071: PetscCall(VecGetArrayRead(patch->intFacetMats, &elementTensors));
2073: for (i = 0; i < numIntFacets; i++) {
2074: const PetscInt facet = intFacetsArray[i + intFacetOffset];
2075: const PetscScalar *v = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart] * nFacetDof * nFacetDof;
2076: idx = 0;
2077: /*
2078: 0--1
2079: |\-|
2080: |+\|
2081: 2--3
2082: [0, 2, 3, 0, 1, 3]
2083: */
2084: for (c = 0; c < 2; c++) {
2085: const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
2086: for (d = 0; d < patch->totalDofsPerCell; d++) {
2087: facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2088: idx++;
2089: }
2090: }
2091: PetscCall(MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES));
2092: }
2093: PetscCall(VecRestoreArrayRead(patch->intFacetMats, &elementTensors));
2094: } else {
2095: /*
2096: 0--1
2097: |\-|
2098: |+\|
2099: 2--3
2100: [0, 2, 3, 0, 1, 3]
2101: */
2102: for (i = 0; i < numIntFacets; i++) {
2103: for (c = 0; c < 2; c++) {
2104: const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
2105: for (d = 0; d < patch->totalDofsPerCell; d++) {
2106: facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2107: if (dofsArrayWithAll) facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell) * patch->totalDofsPerCell + d];
2108: idx++;
2109: }
2110: }
2111: }
2112: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS));
2113: PetscCall(patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2 * numIntFacets * patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx));
2114: PetscCall(ISDestroy(&facetIS));
2115: }
2116: PetscCall(ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells));
2117: PetscCall(ISRestoreIndices(patch->intFacets, &intFacetsArray));
2118: PetscCall(PetscFree(facetDofs));
2119: PetscCall(PetscFree(facetDofsWithAll));
2120: PetscCall(DMDestroy(&dm));
2121: }
2122: }
2124: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2125: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2127: if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2128: MatFactorInfo info;
2129: PetscBool flg;
2130: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg));
2131: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Invalid Mat type for dense inverse");
2132: PetscCall(MatFactorInfoInitialize(&info));
2133: PetscCall(MatLUFactor(mat, NULL, NULL, &info));
2134: PetscCall(MatSeqDenseInvertFactors_Private(mat));
2135: }
2136: PetscCall(ISDestroy(&patch->cellIS));
2137: if (withArtificial) {
2138: PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
2139: } else {
2140: PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
2141: }
2142: if (isNonlinear) PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
2143: PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
2144: if (patch->viewMatrix) {
2145: char name[PETSC_MAX_PATH_LEN];
2147: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch matrix for Point %" PetscInt_FMT, point));
2148: PetscCall(PetscObjectSetName((PetscObject)mat, name));
2149: PetscCall(ObjectView((PetscObject)mat, patch->viewerMatrix, patch->formatMatrix));
2150: }
2151: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2152: PetscFunctionReturn(PETSC_SUCCESS);
2153: }
2155: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2156: {
2157: Vec data;
2158: PetscScalar *array;
2159: PetscInt bs, nz, i, j, cell;
2161: PetscCall(MatShellGetContext(mat, &data));
2162: PetscCall(VecGetBlockSize(data, &bs));
2163: PetscCall(VecGetSize(data, &nz));
2164: PetscCall(VecGetArray(data, &array));
2165: PetscCheck(m == n, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2166: cell = (PetscInt)(idxm[0] / bs); /* use the fact that this is called once per cell */
2167: for (i = 0; i < m; i++) {
2168: PetscCheck(idxm[i] == idxn[i], PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2169: for (j = 0; j < n; j++) {
2170: const PetscScalar v_ = v[i * bs + j];
2171: /* Indexing is special to the data structure we have! */
2172: if (addv == INSERT_VALUES) {
2173: array[cell * bs * bs + i * bs + j] = v_;
2174: } else {
2175: array[cell * bs * bs + i * bs + j] += v_;
2176: }
2177: }
2178: }
2179: PetscCall(VecRestoreArray(data, &array));
2180: PetscFunctionReturn(PETSC_SUCCESS);
2181: }
2183: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2184: {
2185: PC_PATCH *patch = (PC_PATCH *)pc->data;
2186: const PetscInt *cellsArray;
2187: PetscInt ncell, offset;
2188: const PetscInt *dofMapArray;
2189: PetscInt i, j;
2190: IS dofMap;
2191: IS cellIS;
2192: const PetscInt ndof = patch->totalDofsPerCell;
2193: Mat vecMat;
2194: PetscInt cStart, cEnd;
2195: DM dm, plex;
2197: PetscCall(ISGetSize(patch->cells, &ncell));
2198: if (!ncell) { /* No cells to assemble over -> skip */
2199: PetscFunctionReturn(PETSC_SUCCESS);
2200: }
2202: PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
2204: PetscCall(PCGetDM(pc, &dm));
2205: PetscCall(DMConvert(dm, DMPLEX, &plex));
2206: dm = plex;
2207: if (!patch->allCells) {
2208: PetscHSetI cells;
2209: PetscHashIter hi;
2210: PetscInt pStart, pEnd;
2211: PetscInt *allCells = NULL;
2212: PetscCall(PetscHSetICreate(&cells));
2213: PetscCall(ISGetIndices(patch->cells, &cellsArray));
2214: PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
2215: for (i = pStart; i < pEnd; i++) {
2216: PetscCall(PetscSectionGetDof(patch->cellCounts, i, &ncell));
2217: PetscCall(PetscSectionGetOffset(patch->cellCounts, i, &offset));
2218: if (ncell <= 0) continue;
2219: for (j = 0; j < ncell; j++) PetscCall(PetscHSetIAdd(cells, cellsArray[offset + j]));
2220: }
2221: PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
2222: PetscCall(PetscHSetIGetSize(cells, &ncell));
2223: PetscCall(PetscMalloc1(ncell, &allCells));
2224: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2225: PetscCall(PetscMalloc1(cEnd - cStart, &patch->precomputedTensorLocations));
2226: i = 0;
2227: PetscHashIterBegin(cells, hi);
2228: while (!PetscHashIterAtEnd(cells, hi)) {
2229: PetscHashIterGetKey(cells, hi, allCells[i]);
2230: patch->precomputedTensorLocations[allCells[i]] = i;
2231: PetscHashIterNext(cells, hi);
2232: i++;
2233: }
2234: PetscCall(PetscHSetIDestroy(&cells));
2235: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells));
2236: }
2237: PetscCall(ISGetSize(patch->allCells, &ncell));
2238: if (!patch->cellMats) {
2239: PetscCall(VecCreateSeq(PETSC_COMM_SELF, ncell * ndof * ndof, &patch->cellMats));
2240: PetscCall(VecSetBlockSize(patch->cellMats, ndof));
2241: }
2242: PetscCall(VecSet(patch->cellMats, 0));
2244: PetscCall(MatCreateShell(PETSC_COMM_SELF, ncell * ndof, ncell * ndof, ncell * ndof, ncell * ndof, (void *)patch->cellMats, &vecMat));
2245: PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void (*)(void)) & MatSetValues_PCPatch_Private));
2246: PetscCall(ISGetSize(patch->allCells, &ncell));
2247: PetscCall(ISCreateStride(PETSC_COMM_SELF, ndof * ncell, 0, 1, &dofMap));
2248: PetscCall(ISGetIndices(dofMap, &dofMapArray));
2249: PetscCall(ISGetIndices(patch->allCells, &cellsArray));
2250: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS));
2251: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2252: PetscCallBack("PCPatch callback", patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof * ncell, dofMapArray, NULL, patch->usercomputeopctx));
2253: PetscCall(ISDestroy(&cellIS));
2254: PetscCall(MatDestroy(&vecMat));
2255: PetscCall(ISRestoreIndices(patch->allCells, &cellsArray));
2256: PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
2257: PetscCall(ISDestroy(&dofMap));
2259: if (patch->usercomputeopintfacet) {
2260: PetscInt nIntFacets;
2261: IS intFacetsIS;
2262: const PetscInt *intFacetsArray = NULL;
2263: if (!patch->allIntFacets) {
2264: PetscHSetI facets;
2265: PetscHashIter hi;
2266: PetscInt pStart, pEnd, fStart, fEnd;
2267: PetscInt *allIntFacets = NULL;
2268: PetscCall(PetscHSetICreate(&facets));
2269: PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2270: PetscCall(PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd));
2271: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2272: for (i = pStart; i < pEnd; i++) {
2273: PetscCall(PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets));
2274: PetscCall(PetscSectionGetOffset(patch->intFacetCounts, i, &offset));
2275: if (nIntFacets <= 0) continue;
2276: for (j = 0; j < nIntFacets; j++) PetscCall(PetscHSetIAdd(facets, intFacetsArray[offset + j]));
2277: }
2278: PetscCall(ISRestoreIndices(patch->intFacets, &intFacetsArray));
2279: PetscCall(PetscHSetIGetSize(facets, &nIntFacets));
2280: PetscCall(PetscMalloc1(nIntFacets, &allIntFacets));
2281: PetscCall(PetscMalloc1(fEnd - fStart, &patch->precomputedIntFacetTensorLocations));
2282: i = 0;
2283: PetscHashIterBegin(facets, hi);
2284: while (!PetscHashIterAtEnd(facets, hi)) {
2285: PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2286: patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2287: PetscHashIterNext(facets, hi);
2288: i++;
2289: }
2290: PetscCall(PetscHSetIDestroy(&facets));
2291: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets));
2292: }
2293: PetscCall(ISGetSize(patch->allIntFacets, &nIntFacets));
2294: if (!patch->intFacetMats) {
2295: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nIntFacets * ndof * ndof * 4, &patch->intFacetMats));
2296: PetscCall(VecSetBlockSize(patch->intFacetMats, ndof * 2));
2297: }
2298: PetscCall(VecSet(patch->intFacetMats, 0));
2300: PetscCall(MatCreateShell(PETSC_COMM_SELF, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, (void *)patch->intFacetMats, &vecMat));
2301: PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void (*)(void)) & MatSetValues_PCPatch_Private));
2302: PetscCall(ISCreateStride(PETSC_COMM_SELF, 2 * ndof * nIntFacets, 0, 1, &dofMap));
2303: PetscCall(ISGetIndices(dofMap, &dofMapArray));
2304: PetscCall(ISGetIndices(patch->allIntFacets, &intFacetsArray));
2305: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS));
2306: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2307: PetscCallBack("PCPatch callback (interior facets)", patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2 * ndof * nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx));
2308: PetscCall(ISDestroy(&intFacetsIS));
2309: PetscCall(MatDestroy(&vecMat));
2310: PetscCall(ISRestoreIndices(patch->allIntFacets, &intFacetsArray));
2311: PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
2312: PetscCall(ISDestroy(&dofMap));
2313: }
2314: PetscCall(DMDestroy(&dm));
2315: PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2317: PetscFunctionReturn(PETSC_SUCCESS);
2318: }
2320: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2321: {
2322: PC_PATCH *patch = (PC_PATCH *)pc->data;
2323: const PetscScalar *xArray = NULL;
2324: PetscScalar *yArray = NULL;
2325: const PetscInt *gtolArray = NULL;
2326: PetscInt dof, offset, lidx;
2328: PetscFunctionBeginHot;
2329: PetscCall(VecGetArrayRead(x, &xArray));
2330: PetscCall(VecGetArray(y, &yArray));
2331: if (scattertype == SCATTER_WITHARTIFICIAL) {
2332: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof));
2333: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset));
2334: PetscCall(ISGetIndices(patch->gtolWithArtificial, >olArray));
2335: } else if (scattertype == SCATTER_WITHALL) {
2336: PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof));
2337: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset));
2338: PetscCall(ISGetIndices(patch->gtolWithAll, >olArray));
2339: } else {
2340: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
2341: PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2342: PetscCall(ISGetIndices(patch->gtol, >olArray));
2343: }
2344: PetscCheck(mode != INSERT_VALUES || scat == SCATTER_FORWARD, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward");
2345: PetscCheck(mode != ADD_VALUES || scat == SCATTER_REVERSE, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse");
2346: for (lidx = 0; lidx < dof; ++lidx) {
2347: const PetscInt gidx = gtolArray[offset + lidx];
2349: if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
2350: else yArray[gidx] += xArray[lidx]; /* Reverse */
2351: }
2352: if (scattertype == SCATTER_WITHARTIFICIAL) {
2353: PetscCall(ISRestoreIndices(patch->gtolWithArtificial, >olArray));
2354: } else if (scattertype == SCATTER_WITHALL) {
2355: PetscCall(ISRestoreIndices(patch->gtolWithAll, >olArray));
2356: } else {
2357: PetscCall(ISRestoreIndices(patch->gtol, >olArray));
2358: }
2359: PetscCall(VecRestoreArrayRead(x, &xArray));
2360: PetscCall(VecRestoreArray(y, &yArray));
2361: PetscFunctionReturn(PETSC_SUCCESS);
2362: }
2364: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2365: {
2366: PC_PATCH *patch = (PC_PATCH *)pc->data;
2367: const char *prefix;
2368: PetscInt i;
2370: PetscFunctionBegin;
2371: if (!pc->setupcalled) {
2372: PetscCheck(patch->save_operators || !patch->denseinverse, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
2373: if (!patch->denseinverse) {
2374: PetscCall(PetscMalloc1(patch->npatch, &patch->solver));
2375: PetscCall(PCGetOptionsPrefix(pc, &prefix));
2376: for (i = 0; i < patch->npatch; ++i) {
2377: KSP ksp;
2378: PC subpc;
2380: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
2381: PetscCall(KSPSetNestLevel(ksp, pc->kspnestlevel));
2382: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
2383: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
2384: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
2385: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
2386: PetscCall(KSPGetPC(ksp, &subpc));
2387: PetscCall(PetscObjectIncrementTabLevel((PetscObject)subpc, (PetscObject)pc, 1));
2388: patch->solver[i] = (PetscObject)ksp;
2389: }
2390: }
2391: }
2392: if (patch->save_operators) {
2393: if (patch->precomputeElementTensors) PetscCall(PCPatchPrecomputePatchTensors_Private(pc));
2394: for (i = 0; i < patch->npatch; ++i) {
2395: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE));
2396: if (!patch->denseinverse) {
2397: PetscCall(KSPSetOperators((KSP)patch->solver[i], patch->mat[i], patch->mat[i]));
2398: } else if (patch->mat[i] && !patch->densesolve) {
2399: /* Setup matmult callback */
2400: PetscCall(MatGetOperation(patch->mat[i], MATOP_MULT, (void (**)(void)) & patch->densesolve));
2401: }
2402: }
2403: }
2404: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2405: for (i = 0; i < patch->npatch; ++i) {
2406: /* Instead of padding patch->patchUpdate with zeros to get */
2407: /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2408: /* just get rid of the columns that correspond to the dofs with */
2409: /* artificial bcs. That's of course fairly inefficient, hopefully we */
2410: /* can just assemble the rectangular matrix in the first place. */
2411: Mat matSquare;
2412: IS rowis;
2413: PetscInt dof;
2415: PetscCall(MatGetSize(patch->mat[i], &dof, NULL));
2416: if (dof == 0) {
2417: patch->matWithArtificial[i] = NULL;
2418: continue;
2419: }
2421: PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
2422: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));
2424: PetscCall(MatGetSize(matSquare, &dof, NULL));
2425: PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
2426: if (pc->setupcalled) {
2427: PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]));
2428: } else {
2429: PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]));
2430: }
2431: PetscCall(ISDestroy(&rowis));
2432: PetscCall(MatDestroy(&matSquare));
2433: }
2434: }
2435: PetscFunctionReturn(PETSC_SUCCESS);
2436: }
2438: static PetscErrorCode PCSetUp_PATCH(PC pc)
2439: {
2440: PC_PATCH *patch = (PC_PATCH *)pc->data;
2441: PetscInt i;
2442: PetscBool isNonlinear;
2443: PetscInt maxDof = -1, maxDofWithArtificial = -1;
2445: PetscFunctionBegin;
2446: if (!pc->setupcalled) {
2447: PetscInt pStart, pEnd, p;
2448: PetscInt localSize;
2450: PetscCall(PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0));
2452: isNonlinear = patch->isNonlinear;
2453: if (!patch->nsubspaces) {
2454: DM dm, plex;
2455: PetscSection s;
2456: PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;
2458: PetscCall(PCGetDM(pc, &dm));
2459: PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2460: PetscCall(DMConvert(dm, DMPLEX, &plex));
2461: dm = plex;
2462: PetscCall(DMGetLocalSection(dm, &s));
2463: PetscCall(PetscSectionGetNumFields(s, &Nf));
2464: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2465: for (p = pStart; p < pEnd; ++p) {
2466: PetscInt cdof;
2467: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2468: numGlobalBcs += cdof;
2469: }
2470: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2471: PetscCall(PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs));
2472: for (f = 0; f < Nf; ++f) {
2473: PetscFE fe;
2474: PetscDualSpace sp;
2475: PetscInt cdoff = 0;
2477: PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
2478: /* PetscCall(PetscFEGetNumComponents(fe, &Nc[f])); */
2479: PetscCall(PetscFEGetDualSpace(fe, &sp));
2480: PetscCall(PetscDualSpaceGetDimension(sp, &Nb[f]));
2482: PetscCall(PetscMalloc1((cEnd - cStart) * Nb[f], &cellDofs[f]));
2483: for (c = cStart; c < cEnd; ++c) {
2484: PetscInt *closure = NULL;
2485: PetscInt clSize = 0, cl;
2487: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2488: for (cl = 0; cl < clSize * 2; cl += 2) {
2489: const PetscInt p = closure[cl];
2490: PetscInt fdof, d, foff;
2492: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2493: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2494: for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2495: }
2496: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2497: }
2498: 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]);
2499: }
2500: numGlobalBcs = 0;
2501: for (p = pStart; p < pEnd; ++p) {
2502: const PetscInt *ind;
2503: PetscInt off, cdof, d;
2505: PetscCall(PetscSectionGetOffset(s, p, &off));
2506: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2507: PetscCall(PetscSectionGetConstraintIndices(s, p, &ind));
2508: for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2509: }
2511: PetscCall(PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **)cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs));
2512: for (f = 0; f < Nf; ++f) PetscCall(PetscFree(cellDofs[f]));
2513: PetscCall(PetscFree3(Nb, cellDofs, globalBcs));
2514: PetscCall(PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL));
2515: PetscCall(PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL));
2516: PetscCall(DMDestroy(&dm));
2517: }
2519: localSize = patch->subspaceOffsets[patch->nsubspaces];
2520: PetscCall(VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS));
2521: PetscCall(VecSetUp(patch->localRHS));
2522: PetscCall(VecDuplicate(patch->localRHS, &patch->localUpdate));
2523: PetscCall(PCPatchCreateCellPatches(pc));
2524: PetscCall(PCPatchCreateCellPatchDiscretisationInfo(pc));
2526: /* OK, now build the work vectors */
2527: PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd));
2529: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial));
2530: if (isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll));
2531: for (p = pStart; p < pEnd; ++p) {
2532: PetscInt dof;
2534: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
2535: maxDof = PetscMax(maxDof, dof);
2536: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2537: const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
2538: PetscInt numPatchDofs, offset;
2539: PetscInt numPatchDofsWithArtificial, offsetWithArtificial;
2540: PetscInt dofWithoutArtificialCounter = 0;
2541: PetscInt *patchWithoutArtificialToWithArtificialArray;
2543: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof));
2544: maxDofWithArtificial = PetscMax(maxDofWithArtificial, dof);
2546: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
2547: /* the index in the patch with all dofs */
2548: PetscCall(ISGetIndices(patch->gtol, >olArray));
2550: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
2551: if (numPatchDofs == 0) {
2552: patch->dofMappingWithoutToWithArtificial[p - pStart] = NULL;
2553: continue;
2554: }
2556: PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2557: PetscCall(ISGetIndices(patch->gtolWithArtificial, >olArrayWithArtificial));
2558: PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial));
2559: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial));
2561: PetscCall(PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray));
2562: for (i = 0; i < numPatchDofsWithArtificial; i++) {
2563: if (gtolArrayWithArtificial[i + offsetWithArtificial] == gtolArray[offset + dofWithoutArtificialCounter]) {
2564: patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2565: dofWithoutArtificialCounter++;
2566: if (dofWithoutArtificialCounter == numPatchDofs) break;
2567: }
2568: }
2569: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p - pStart]));
2570: PetscCall(ISRestoreIndices(patch->gtol, >olArray));
2571: PetscCall(ISRestoreIndices(patch->gtolWithArtificial, >olArrayWithArtificial));
2572: }
2573: }
2574: for (p = pStart; p < pEnd; ++p) {
2575: if (isNonlinear) {
2576: const PetscInt *gtolArray, *gtolArrayWithAll = NULL;
2577: PetscInt numPatchDofs, offset;
2578: PetscInt numPatchDofsWithAll, offsetWithAll;
2579: PetscInt dofWithoutAllCounter = 0;
2580: PetscInt *patchWithoutAllToWithAllArray;
2582: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
2583: /* the index in the patch with all dofs */
2584: PetscCall(ISGetIndices(patch->gtol, >olArray));
2586: PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
2587: if (numPatchDofs == 0) {
2588: patch->dofMappingWithoutToWithAll[p - pStart] = NULL;
2589: continue;
2590: }
2592: PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2593: PetscCall(ISGetIndices(patch->gtolWithAll, >olArrayWithAll));
2594: PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll));
2595: PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll));
2597: PetscCall(PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray));
2599: for (i = 0; i < numPatchDofsWithAll; i++) {
2600: if (gtolArrayWithAll[i + offsetWithAll] == gtolArray[offset + dofWithoutAllCounter]) {
2601: patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2602: dofWithoutAllCounter++;
2603: if (dofWithoutAllCounter == numPatchDofs) break;
2604: }
2605: }
2606: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p - pStart]));
2607: PetscCall(ISRestoreIndices(patch->gtol, >olArray));
2608: PetscCall(ISRestoreIndices(patch->gtolWithAll, >olArrayWithAll));
2609: }
2610: }
2611: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2612: PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial));
2613: PetscCall(VecSetUp(patch->patchRHSWithArtificial));
2614: }
2615: PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS));
2616: PetscCall(VecSetUp(patch->patchRHS));
2617: PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate));
2618: PetscCall(VecSetUp(patch->patchUpdate));
2619: if (patch->save_operators) {
2620: PetscCall(PetscMalloc1(patch->npatch, &patch->mat));
2621: for (i = 0; i < patch->npatch; ++i) PetscCall(PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE));
2622: }
2623: PetscCall(PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0));
2625: /* If desired, calculate weights for dof multiplicity */
2626: if (patch->partition_of_unity) {
2627: PetscScalar *input = NULL;
2628: PetscScalar *output = NULL;
2629: Vec global;
2631: PetscCall(VecDuplicate(patch->localRHS, &patch->dof_weights));
2632: if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2633: for (i = 0; i < patch->npatch; ++i) {
2634: PetscInt dof;
2636: PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &dof));
2637: if (dof <= 0) continue;
2638: PetscCall(VecSet(patch->patchRHS, 1.0));
2639: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
2640: }
2641: } else {
2642: /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2643: PetscCall(VecSet(patch->dof_weights, 1.0));
2644: }
2646: PetscCall(VecDuplicate(patch->dof_weights, &global));
2647: PetscCall(VecSet(global, 0.));
2649: PetscCall(VecGetArray(patch->dof_weights, &input));
2650: PetscCall(VecGetArray(global, &output));
2651: PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
2652: PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
2653: PetscCall(VecRestoreArray(patch->dof_weights, &input));
2654: PetscCall(VecRestoreArray(global, &output));
2656: PetscCall(VecReciprocal(global));
2658: PetscCall(VecGetArray(patch->dof_weights, &output));
2659: PetscCall(VecGetArray(global, &input));
2660: PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
2661: PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
2662: PetscCall(VecRestoreArray(patch->dof_weights, &output));
2663: PetscCall(VecRestoreArray(global, &input));
2664: PetscCall(VecDestroy(&global));
2665: }
2666: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators && !patch->isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->matWithArtificial));
2667: }
2668: PetscCall((*patch->setupsolver)(pc));
2669: PetscFunctionReturn(PETSC_SUCCESS);
2670: }
2672: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2673: {
2674: PC_PATCH *patch = (PC_PATCH *)pc->data;
2675: KSP ksp;
2676: Mat op;
2677: PetscInt m, n;
2679: PetscFunctionBegin;
2680: if (patch->denseinverse) {
2681: PetscCall((*patch->densesolve)(patch->mat[i], x, y));
2682: PetscFunctionReturn(PETSC_SUCCESS);
2683: }
2684: ksp = (KSP)patch->solver[i];
2685: if (!patch->save_operators) {
2686: Mat mat;
2688: PetscCall(PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE));
2689: /* Populate operator here. */
2690: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE));
2691: PetscCall(KSPSetOperators(ksp, mat, mat));
2692: /* Drop reference so the KSPSetOperators below will blow it away. */
2693: PetscCall(MatDestroy(&mat));
2694: }
2695: PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
2696: if (!ksp->setfromoptionscalled) PetscCall(KSPSetFromOptions(ksp));
2697: /* Disgusting trick to reuse work vectors */
2698: PetscCall(KSPGetOperators(ksp, &op, NULL));
2699: PetscCall(MatGetLocalSize(op, &m, &n));
2700: x->map->n = m;
2701: y->map->n = n;
2702: x->map->N = m;
2703: y->map->N = n;
2704: PetscCall(KSPSolve(ksp, x, y));
2705: PetscCall(KSPCheckSolve(ksp, pc, y));
2706: PetscCall(PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0));
2707: if (!patch->save_operators) {
2708: PC pc;
2709: PetscCall(KSPSetOperators(ksp, NULL, NULL));
2710: PetscCall(KSPGetPC(ksp, &pc));
2711: /* Destroy PC context too, otherwise the factored matrix hangs around. */
2712: PetscCall(PCReset(pc));
2713: }
2714: PetscFunctionReturn(PETSC_SUCCESS);
2715: }
2717: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2718: {
2719: PC_PATCH *patch = (PC_PATCH *)pc->data;
2720: Mat multMat;
2721: PetscInt n, m;
2723: PetscFunctionBegin;
2725: if (patch->save_operators) {
2726: multMat = patch->matWithArtificial[i];
2727: } else {
2728: /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2729: Mat matSquare;
2730: PetscInt dof;
2731: IS rowis;
2732: PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
2733: PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));
2734: PetscCall(MatGetSize(matSquare, &dof, NULL));
2735: PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
2736: PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat));
2737: PetscCall(MatDestroy(&matSquare));
2738: PetscCall(ISDestroy(&rowis));
2739: }
2740: /* Disgusting trick to reuse work vectors */
2741: PetscCall(MatGetLocalSize(multMat, &m, &n));
2742: patch->patchUpdate->map->n = n;
2743: patch->patchRHSWithArtificial->map->n = m;
2744: patch->patchUpdate->map->N = n;
2745: patch->patchRHSWithArtificial->map->N = m;
2746: PetscCall(MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial));
2747: PetscCall(VecScale(patch->patchRHSWithArtificial, -1.0));
2748: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL));
2749: if (!patch->save_operators) PetscCall(MatDestroy(&multMat));
2750: PetscFunctionReturn(PETSC_SUCCESS);
2751: }
2753: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2754: {
2755: PC_PATCH *patch = (PC_PATCH *)pc->data;
2756: const PetscScalar *globalRHS = NULL;
2757: PetscScalar *localRHS = NULL;
2758: PetscScalar *globalUpdate = NULL;
2759: const PetscInt *bcNodes = NULL;
2760: PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1;
2761: PetscInt start[2] = {0, 0};
2762: PetscInt end[2] = {-1, -1};
2763: const PetscInt inc[2] = {1, -1};
2764: const PetscScalar *localUpdate;
2765: const PetscInt *iterationSet;
2766: PetscInt pStart, numBcs, n, sweep, bc, j;
2768: PetscFunctionBegin;
2769: PetscCall(PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0));
2770: PetscCall(PetscOptionsPushGetViewerOff(PETSC_TRUE));
2771: /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2772: end[0] = patch->npatch;
2773: start[1] = patch->npatch - 1;
2774: if (patch->user_patches) {
2775: PetscCall(ISGetLocalSize(patch->iterationSet, &end[0]));
2776: start[1] = end[0] - 1;
2777: PetscCall(ISGetIndices(patch->iterationSet, &iterationSet));
2778: }
2779: /* Scatter from global space into overlapped local spaces */
2780: PetscCall(VecGetArrayRead(x, &globalRHS));
2781: PetscCall(VecGetArray(patch->localRHS, &localRHS));
2782: PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
2783: PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
2784: PetscCall(VecRestoreArrayRead(x, &globalRHS));
2785: PetscCall(VecRestoreArray(patch->localRHS, &localRHS));
2787: PetscCall(VecSet(patch->localUpdate, 0.0));
2788: PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
2789: PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
2790: for (sweep = 0; sweep < nsweep; sweep++) {
2791: for (j = start[sweep]; j * inc[sweep] < end[sweep] * inc[sweep]; j += inc[sweep]) {
2792: PetscInt i = patch->user_patches ? iterationSet[j] : j;
2793: PetscInt start, len;
2795: PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &len));
2796: PetscCall(PetscSectionGetOffset(patch->gtolCounts, i + pStart, &start));
2797: /* TODO: Squash out these guys in the setup as well. */
2798: if (len <= 0) continue;
2799: /* TODO: Do we need different scatters for X and Y? */
2800: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->localRHS, patch->patchRHS, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR));
2801: PetscCall((*patch->applysolver)(pc, i, patch->patchRHS, patch->patchUpdate));
2802: PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchUpdate, patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
2803: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall((*patch->updatemultiplicative)(pc, i, pStart));
2804: }
2805: }
2806: PetscCall(PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0));
2807: if (patch->user_patches) PetscCall(ISRestoreIndices(patch->iterationSet, &iterationSet));
2808: /* XXX: should we do this on the global vector? */
2809: if (patch->partition_of_unity) PetscCall(VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights));
2810: /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
2811: PetscCall(VecSet(y, 0.0));
2812: PetscCall(VecGetArray(y, &globalUpdate));
2813: PetscCall(VecGetArrayRead(patch->localUpdate, &localUpdate));
2814: PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM));
2815: PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM));
2816: PetscCall(VecRestoreArrayRead(patch->localUpdate, &localUpdate));
2818: /* Now we need to send the global BC values through */
2819: PetscCall(VecGetArrayRead(x, &globalRHS));
2820: PetscCall(ISGetSize(patch->globalBcNodes, &numBcs));
2821: PetscCall(ISGetIndices(patch->globalBcNodes, &bcNodes));
2822: PetscCall(VecGetLocalSize(x, &n));
2823: for (bc = 0; bc < numBcs; ++bc) {
2824: const PetscInt idx = bcNodes[bc];
2825: if (idx < n) globalUpdate[idx] = globalRHS[idx];
2826: }
2828: PetscCall(ISRestoreIndices(patch->globalBcNodes, &bcNodes));
2829: PetscCall(VecRestoreArrayRead(x, &globalRHS));
2830: PetscCall(VecRestoreArray(y, &globalUpdate));
2832: PetscCall(PetscOptionsPopGetViewerOff());
2833: PetscCall(PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0));
2834: PetscFunctionReturn(PETSC_SUCCESS);
2835: }
2837: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2838: {
2839: PC_PATCH *patch = (PC_PATCH *)pc->data;
2840: PetscInt i;
2842: PetscFunctionBegin;
2843: if (patch->solver) {
2844: for (i = 0; i < patch->npatch; ++i) PetscCall(KSPReset((KSP)patch->solver[i]));
2845: }
2846: PetscFunctionReturn(PETSC_SUCCESS);
2847: }
2849: static PetscErrorCode PCReset_PATCH(PC pc)
2850: {
2851: PC_PATCH *patch = (PC_PATCH *)pc->data;
2852: PetscInt i;
2854: PetscFunctionBegin;
2856: PetscCall(PetscSFDestroy(&patch->sectionSF));
2857: PetscCall(PetscSectionDestroy(&patch->cellCounts));
2858: PetscCall(PetscSectionDestroy(&patch->pointCounts));
2859: PetscCall(PetscSectionDestroy(&patch->cellNumbering));
2860: PetscCall(PetscSectionDestroy(&patch->gtolCounts));
2861: PetscCall(ISDestroy(&patch->gtol));
2862: PetscCall(ISDestroy(&patch->cells));
2863: PetscCall(ISDestroy(&patch->points));
2864: PetscCall(ISDestroy(&patch->dofs));
2865: PetscCall(ISDestroy(&patch->offs));
2866: PetscCall(PetscSectionDestroy(&patch->patchSection));
2867: PetscCall(ISDestroy(&patch->ghostBcNodes));
2868: PetscCall(ISDestroy(&patch->globalBcNodes));
2869: PetscCall(PetscSectionDestroy(&patch->gtolCountsWithArtificial));
2870: PetscCall(ISDestroy(&patch->gtolWithArtificial));
2871: PetscCall(ISDestroy(&patch->dofsWithArtificial));
2872: PetscCall(ISDestroy(&patch->offsWithArtificial));
2873: PetscCall(PetscSectionDestroy(&patch->gtolCountsWithAll));
2874: PetscCall(ISDestroy(&patch->gtolWithAll));
2875: PetscCall(ISDestroy(&patch->dofsWithAll));
2876: PetscCall(ISDestroy(&patch->offsWithAll));
2877: PetscCall(VecDestroy(&patch->cellMats));
2878: PetscCall(VecDestroy(&patch->intFacetMats));
2879: PetscCall(ISDestroy(&patch->allCells));
2880: PetscCall(ISDestroy(&patch->intFacets));
2881: PetscCall(ISDestroy(&patch->extFacets));
2882: PetscCall(ISDestroy(&patch->intFacetsToPatchCell));
2883: PetscCall(PetscSectionDestroy(&patch->intFacetCounts));
2884: PetscCall(PetscSectionDestroy(&patch->extFacetCounts));
2886: if (patch->dofSection)
2887: for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscSectionDestroy(&patch->dofSection[i]));
2888: PetscCall(PetscFree(patch->dofSection));
2889: PetscCall(PetscFree(patch->bs));
2890: PetscCall(PetscFree(patch->nodesPerCell));
2891: if (patch->cellNodeMap)
2892: for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscFree(patch->cellNodeMap[i]));
2893: PetscCall(PetscFree(patch->cellNodeMap));
2894: PetscCall(PetscFree(patch->subspaceOffsets));
2896: PetscCall((*patch->resetsolver)(pc));
2898: if (patch->subspaces_to_exclude) PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));
2900: PetscCall(VecDestroy(&patch->localRHS));
2901: PetscCall(VecDestroy(&patch->localUpdate));
2902: PetscCall(VecDestroy(&patch->patchRHS));
2903: PetscCall(VecDestroy(&patch->patchUpdate));
2904: PetscCall(VecDestroy(&patch->dof_weights));
2905: if (patch->patch_dof_weights) {
2906: for (i = 0; i < patch->npatch; ++i) PetscCall(VecDestroy(&patch->patch_dof_weights[i]));
2907: PetscCall(PetscFree(patch->patch_dof_weights));
2908: }
2909: if (patch->mat) {
2910: for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->mat[i]));
2911: PetscCall(PetscFree(patch->mat));
2912: }
2913: if (patch->matWithArtificial && !patch->isNonlinear) {
2914: for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->matWithArtificial[i]));
2915: PetscCall(PetscFree(patch->matWithArtificial));
2916: }
2917: PetscCall(VecDestroy(&patch->patchRHSWithArtificial));
2918: if (patch->dofMappingWithoutToWithArtificial) {
2919: for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]));
2920: PetscCall(PetscFree(patch->dofMappingWithoutToWithArtificial));
2921: }
2922: if (patch->dofMappingWithoutToWithAll) {
2923: for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithAll[i]));
2924: PetscCall(PetscFree(patch->dofMappingWithoutToWithAll));
2925: }
2926: PetscCall(PetscFree(patch->sub_mat_type));
2927: if (patch->userIS) {
2928: for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->userIS[i]));
2929: PetscCall(PetscFree(patch->userIS));
2930: }
2931: PetscCall(PetscFree(patch->precomputedTensorLocations));
2932: PetscCall(PetscFree(patch->precomputedIntFacetTensorLocations));
2934: patch->bs = NULL;
2935: patch->cellNodeMap = NULL;
2936: patch->nsubspaces = 0;
2937: PetscCall(ISDestroy(&patch->iterationSet));
2939: PetscCall(PetscViewerDestroy(&patch->viewerSection));
2940: PetscFunctionReturn(PETSC_SUCCESS);
2941: }
2943: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
2944: {
2945: PC_PATCH *patch = (PC_PATCH *)pc->data;
2946: PetscInt i;
2948: PetscFunctionBegin;
2949: if (patch->solver) {
2950: for (i = 0; i < patch->npatch; ++i) PetscCall(KSPDestroy((KSP *)&patch->solver[i]));
2951: PetscCall(PetscFree(patch->solver));
2952: }
2953: PetscFunctionReturn(PETSC_SUCCESS);
2954: }
2956: static PetscErrorCode PCDestroy_PATCH(PC pc)
2957: {
2958: PC_PATCH *patch = (PC_PATCH *)pc->data;
2960: PetscFunctionBegin;
2961: PetscCall(PCReset_PATCH(pc));
2962: PetscCall((*patch->destroysolver)(pc));
2963: PetscCall(PetscFree(pc->data));
2964: PetscFunctionReturn(PETSC_SUCCESS);
2965: }
2967: static PetscErrorCode PCSetFromOptions_PATCH(PC pc, PetscOptionItems *PetscOptionsObject)
2968: {
2969: PC_PATCH *patch = (PC_PATCH *)pc->data;
2970: PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
2971: char sub_mat_type[PETSC_MAX_PATH_LEN];
2972: char option[PETSC_MAX_PATH_LEN];
2973: const char *prefix;
2974: PetscBool flg, dimflg, codimflg;
2975: MPI_Comm comm;
2976: PetscInt *ifields, nfields, k;
2977: PCCompositeType loctype = PC_COMPOSITE_ADDITIVE;
2979: PetscFunctionBegin;
2980: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
2981: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
2982: PetscOptionsHeadBegin(PetscOptionsObject, "Patch solver options");
2984: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname));
2985: PetscCall(PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg));
2987: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_precompute_element_tensors", patch->classname));
2988: PetscCall(PetscOptionsBool(option, "Compute each element tensor only once?", "PCPatchSetPrecomputeElementTensors", patch->precomputeElementTensors, &patch->precomputeElementTensors, &flg));
2989: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname));
2990: PetscCall(PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg));
2992: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname));
2993: PetscCall(PetscOptionsEnum(option, "Type of local solver composition (additive or multiplicative)", "PCPatchSetLocalComposition", PCCompositeTypes, (PetscEnum)loctype, (PetscEnum *)&loctype, &flg));
2994: if (flg) PetscCall(PCPatchSetLocalComposition(pc, loctype));
2995: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_dense_inverse", patch->classname));
2996: PetscCall(PetscOptionsBool(option, "Compute inverses of patch matrices and apply directly? Ignores KSP/PC settings on patch.", "PCPatchSetDenseInverse", patch->denseinverse, &patch->denseinverse, &flg));
2997: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname));
2998: PetscCall(PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg));
2999: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname));
3000: PetscCall(PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg));
3001: PetscCheck(!dimflg || !codimflg, comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");
3003: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname));
3004: PetscCall(PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum)patchConstructionType, (PetscEnum *)&patchConstructionType, &flg));
3005: if (flg) PetscCall(PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL));
3007: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname));
3008: PetscCall(PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg));
3010: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname));
3011: PetscCall(PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg));
3013: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_pardecomp_overlap", patch->classname));
3014: PetscCall(PetscOptionsInt(option, "What overlap should we use in construct type pardecomp?", "PCPATCH", patch->pardecomp_overlap, &patch->pardecomp_overlap, &flg));
3016: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname));
3017: PetscCall(PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg));
3018: if (flg) PetscCall(PCPatchSetSubMatType(pc, sub_mat_type));
3020: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname));
3021: PetscCall(PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg));
3023: /* If the user has set the number of subspaces, use that for the buffer size,
3024: otherwise use a large number */
3025: if (patch->nsubspaces <= 0) {
3026: nfields = 128;
3027: } else {
3028: nfields = patch->nsubspaces;
3029: }
3030: PetscCall(PetscMalloc1(nfields, &ifields));
3031: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname));
3032: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, option, ifields, &nfields, &flg));
3033: 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");
3034: if (flg) {
3035: PetscCall(PetscHSetIClear(patch->subspaces_to_exclude));
3036: for (k = 0; k < nfields; k++) PetscCall(PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]));
3037: }
3038: PetscCall(PetscFree(ifields));
3040: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname));
3041: PetscCall(PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg));
3042: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname));
3043: PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells));
3044: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname));
3045: PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets));
3046: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname));
3047: PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets));
3048: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname));
3049: PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints));
3050: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname));
3051: PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection));
3052: PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname));
3053: PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix));
3054: PetscOptionsHeadEnd();
3055: patch->optionsSet = PETSC_TRUE;
3056: PetscFunctionReturn(PETSC_SUCCESS);
3057: }
3059: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3060: {
3061: PC_PATCH *patch = (PC_PATCH *)pc->data;
3062: KSPConvergedReason reason;
3063: PetscInt i;
3065: PetscFunctionBegin;
3066: if (!patch->save_operators) {
3067: /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3068: PetscFunctionReturn(PETSC_SUCCESS);
3069: }
3070: if (patch->denseinverse) {
3071: /* No solvers */
3072: PetscFunctionReturn(PETSC_SUCCESS);
3073: }
3074: for (i = 0; i < patch->npatch; ++i) {
3075: if (!((KSP)patch->solver[i])->setfromoptionscalled) PetscCall(KSPSetFromOptions((KSP)patch->solver[i]));
3076: PetscCall(KSPSetUp((KSP)patch->solver[i]));
3077: PetscCall(KSPGetConvergedReason((KSP)patch->solver[i], &reason));
3078: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3079: }
3080: PetscFunctionReturn(PETSC_SUCCESS);
3081: }
3083: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3084: {
3085: PC_PATCH *patch = (PC_PATCH *)pc->data;
3086: PetscViewer sviewer;
3087: PetscBool isascii;
3088: PetscMPIInt rank;
3090: PetscFunctionBegin;
3091: /* TODO Redo tabbing with set tbas in new style */
3092: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
3093: if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);
3094: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
3095: PetscCall(PetscViewerASCIIPushTab(viewer));
3096: PetscCall(PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %" PetscInt_FMT " patches\n", patch->npatch));
3097: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3098: PetscCall(PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n"));
3099: } else {
3100: PetscCall(PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n"));
3101: }
3102: if (patch->partition_of_unity) PetscCall(PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n"));
3103: else PetscCall(PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n"));
3104: if (patch->symmetrise_sweep) PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n"));
3105: else PetscCall(PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n"));
3106: if (!patch->precomputeElementTensors) PetscCall(PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n"));
3107: else PetscCall(PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n"));
3108: if (!patch->save_operators) PetscCall(PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n"));
3109: else PetscCall(PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n"));
3110: if (patch->patchconstructop == PCPatchConstruct_Star) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n"));
3111: else if (patch->patchconstructop == PCPatchConstruct_Vanka) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n"));
3112: else if (patch->patchconstructop == PCPatchConstruct_User) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n"));
3113: else PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n"));
3115: if (patch->denseinverse) {
3116: PetscCall(PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n"));
3117: } else {
3118: if (patch->isNonlinear) {
3119: PetscCall(PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n"));
3120: } else {
3121: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n"));
3122: }
3123: if (patch->solver) {
3124: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
3125: if (rank == 0) {
3126: PetscCall(PetscViewerASCIIPushTab(sviewer));
3127: PetscCall(PetscObjectView(patch->solver[0], sviewer));
3128: PetscCall(PetscViewerASCIIPopTab(sviewer));
3129: }
3130: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
3131: } else {
3132: PetscCall(PetscViewerASCIIPushTab(viewer));
3133: PetscCall(PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n"));
3134: PetscCall(PetscViewerASCIIPopTab(viewer));
3135: }
3136: }
3137: PetscCall(PetscViewerASCIIPopTab(viewer));
3138: PetscFunctionReturn(PETSC_SUCCESS);
3139: }
3141: /*MC
3142: PCPATCH - A `PC` object that encapsulates flexible definition of blocks for overlapping and non-overlapping
3143: small block additive preconditioners. Block definition is based on topology from
3144: a `DM` and equation numbering from a `PetscSection`.
3146: Options Database Keys:
3147: + -pc_patch_cells_view - Views the process local cell numbers for each patch
3148: . -pc_patch_points_view - Views the process local mesh point numbers for each patch
3149: . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch
3150: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3151: - -pc_patch_sub_mat_view - Views the matrix associated with each patch
3153: Level: intermediate
3155: .seealso: [](ch_ksp), `PCType`, `PCCreate()`, `PCSetType()`, `PCASM`, `PCJACOBI`, `PCPBJACOBI`, `PCVPBJACOBI`, `SNESPATCH`
3156: M*/
3157: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3158: {
3159: PC_PATCH *patch;
3161: PetscFunctionBegin;
3162: PetscCall(PetscCitationsRegister(PCPatchCitation, &PCPatchcite));
3163: PetscCall(PetscNew(&patch));
3165: if (patch->subspaces_to_exclude) PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));
3166: PetscCall(PetscHSetICreate(&patch->subspaces_to_exclude));
3168: patch->classname = "pc";
3169: patch->isNonlinear = PETSC_FALSE;
3171: /* Set some defaults */
3172: patch->combined = PETSC_FALSE;
3173: patch->save_operators = PETSC_TRUE;
3174: patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3175: patch->precomputeElementTensors = PETSC_FALSE;
3176: patch->partition_of_unity = PETSC_FALSE;
3177: patch->codim = -1;
3178: patch->dim = -1;
3179: patch->vankadim = -1;
3180: patch->ignoredim = -1;
3181: patch->pardecomp_overlap = 0;
3182: patch->patchconstructop = PCPatchConstruct_Star;
3183: patch->symmetrise_sweep = PETSC_FALSE;
3184: patch->npatch = 0;
3185: patch->userIS = NULL;
3186: patch->optionsSet = PETSC_FALSE;
3187: patch->iterationSet = NULL;
3188: patch->user_patches = PETSC_FALSE;
3189: PetscCall(PetscStrallocpy(MATDENSE, (char **)&patch->sub_mat_type));
3190: patch->viewPatches = PETSC_FALSE;
3191: patch->viewCells = PETSC_FALSE;
3192: patch->viewPoints = PETSC_FALSE;
3193: patch->viewSection = PETSC_FALSE;
3194: patch->viewMatrix = PETSC_FALSE;
3195: patch->densesolve = NULL;
3196: patch->setupsolver = PCSetUp_PATCH_Linear;
3197: patch->applysolver = PCApply_PATCH_Linear;
3198: patch->resetsolver = PCReset_PATCH_Linear;
3199: patch->destroysolver = PCDestroy_PATCH_Linear;
3200: patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3201: patch->dofMappingWithoutToWithArtificial = NULL;
3202: patch->dofMappingWithoutToWithAll = NULL;
3204: pc->data = (void *)patch;
3205: pc->ops->apply = PCApply_PATCH;
3206: pc->ops->applytranspose = NULL; /* PCApplyTranspose_PATCH; */
3207: pc->ops->setup = PCSetUp_PATCH;
3208: pc->ops->reset = PCReset_PATCH;
3209: pc->ops->destroy = PCDestroy_PATCH;
3210: pc->ops->setfromoptions = PCSetFromOptions_PATCH;
3211: pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH;
3212: pc->ops->view = PCView_PATCH;
3213: pc->ops->applyrichardson = NULL;
3215: PetscFunctionReturn(PETSC_SUCCESS);
3216: }