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