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