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