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