Actual source code: patch.c
1: #include <petsc/private/dmpatchimpl.h>
2: #include <petscdmda.h>
3: #include <petscsf.h>
5: /*
6: Solver loop to update \tau:
8: DMZoom(dmc, &dmz)
9: DMRefine(dmz, &dmf),
10: Scatter Xcoarse -> Xzoom,
11: Interpolate Xzoom -> Xfine (note that this may be on subcomms),
12: Smooth Xfine using two-step smoother
13: normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine
14: Compute residual Rfine
15: Restrict Rfine to Rzoom_restricted
16: Scatter Rzoom_restricted -> Rcoarse_restricted
17: Compute global residual Rcoarse
18: TauCoarse = Rcoarse - Rcoarse_restricted
19: */
21: /*@C
22: DMPatchZoom - Create patches of a DMDA on subsets of processes, indicated by commz
24: Collective on dm
26: Input Parameters:
27: + dm - the DM
28: . lower - the lower left corner of the requested patch
29: . upper - the upper right corner of the requested patch
30: - commz - the new communicator for the patch, MPI_COMM_NULL indicates that the given rank will not own a patch
32: Output Parameters:
33: + dmz - the patch DM
34: . sfz - the PetscSF mapping the patch+halo to the zoomed version (optional)
35: - sfzr - the PetscSF mapping the patch to the restricted zoomed version
37: Level: intermediate
39: .seealso: `DMPatchSolve()`, `DMDACreatePatchIS()`
40: @*/
41: PetscErrorCode DMPatchZoom(DM dm, MatStencil lower, MatStencil upper, MPI_Comm commz, DM *dmz, PetscSF *sfz, PetscSF *sfzr)
42: {
43: DMDAStencilType st;
44: MatStencil blower, bupper, loclower, locupper;
45: IS is;
46: const PetscInt *ranges, *indices;
47: PetscInt *localPoints = NULL;
48: PetscSFNode *remotePoints = NULL;
49: PetscInt dim, dof;
50: PetscInt M, N, P, rM, rN, rP, halo = 1, sxb, syb, szb, sxr, syr, szr, exr, eyr, ezr, mxb, myb, mzb, i, j, k, l, q;
51: PetscMPIInt size;
52: PetscBool patchis_offproc = PETSC_TRUE;
53: Vec X;
55: if (!sfz) halo = 0;
56: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
57: /* Create patch DM */
58: DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, &st);
60: /* Get piece for rank r, expanded by halo */
61: bupper.i = PetscMin(M, upper.i + halo);
62: blower.i = PetscMax(lower.i - halo, 0);
63: bupper.j = PetscMin(N, upper.j + halo);
64: blower.j = PetscMax(lower.j - halo, 0);
65: bupper.k = PetscMin(P, upper.k + halo);
66: blower.k = PetscMax(lower.k - halo, 0);
67: rM = bupper.i - blower.i;
68: rN = bupper.j - blower.j;
69: rP = bupper.k - blower.k;
71: if (commz != MPI_COMM_NULL) {
72: DMDACreate(commz, dmz);
73: DMSetDimension(*dmz, dim);
74: DMDASetSizes(*dmz, rM, rN, rP);
75: DMDASetNumProcs(*dmz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
76: DMDASetBoundaryType(*dmz, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE);
77: DMDASetDof(*dmz, dof);
78: DMDASetStencilType(*dmz, st);
79: DMDASetStencilWidth(*dmz, 0);
80: DMDASetOwnershipRanges(*dmz, NULL, NULL, NULL);
81: DMSetFromOptions(*dmz);
82: DMSetUp(*dmz);
83: DMDAGetCorners(*dmz, &sxb, &syb, &szb, &mxb, &myb, &mzb);
84: sxr = PetscMax(sxb, lower.i - blower.i);
85: syr = PetscMax(syb, lower.j - blower.j);
86: szr = PetscMax(szb, lower.k - blower.k);
87: exr = PetscMin(sxb + mxb, upper.i - blower.i);
88: eyr = PetscMin(syb + myb, upper.j - blower.j);
89: ezr = PetscMin(szb + mzb, upper.k - blower.k);
90: PetscMalloc2(dof * rM * rN * PetscMax(rP, 1), &localPoints, dof * rM * rN * PetscMax(rP, 1), &remotePoints);
91: } else {
92: sxr = syr = szr = exr = eyr = ezr = sxb = syb = szb = mxb = myb = mzb = 0;
93: }
95: /* Create SF for restricted map */
96: DMCreateGlobalVector(dm, &X);
97: VecGetOwnershipRanges(X, &ranges);
99: loclower.i = blower.i + sxr;
100: locupper.i = blower.i + exr;
101: loclower.j = blower.j + syr;
102: locupper.j = blower.j + eyr;
103: loclower.k = blower.k + szr;
104: locupper.k = blower.k + ezr;
106: DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc);
107: ISGetIndices(is, &indices);
109: if (dim < 3) {
110: mzb = 1;
111: ezr = 1;
112: }
113: q = 0;
114: for (k = szb; k < szb + mzb; ++k) {
115: if ((k < szr) || (k >= ezr)) continue;
116: for (j = syb; j < syb + myb; ++j) {
117: if ((j < syr) || (j >= eyr)) continue;
118: for (i = sxb; i < sxb + mxb; ++i) {
119: for (l = 0; l < dof; l++) {
120: const PetscInt lp = l + dof * (((k - szb) * rN + (j - syb)) * rM + i - sxb);
121: PetscInt r;
123: if ((i < sxr) || (i >= exr)) continue;
124: localPoints[q] = lp;
125: PetscFindInt(indices[q], size + 1, ranges, &r);
127: remotePoints[q].rank = r < 0 ? -(r + 1) - 1 : r;
128: remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank];
129: ++q;
130: }
131: }
132: }
133: }
134: ISRestoreIndices(is, &indices);
135: ISDestroy(&is);
136: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfzr);
137: PetscObjectSetName((PetscObject)*sfzr, "Restricted Map");
138: PetscSFSetGraph(*sfzr, dof * M * N * P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);
140: if (sfz) {
141: /* Create SF for buffered map */
142: loclower.i = blower.i + sxb;
143: locupper.i = blower.i + sxb + mxb;
144: loclower.j = blower.j + syb;
145: locupper.j = blower.j + syb + myb;
146: loclower.k = blower.k + szb;
147: locupper.k = blower.k + szb + mzb;
149: DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc);
150: ISGetIndices(is, &indices);
152: q = 0;
153: for (k = szb; k < szb + mzb; ++k) {
154: for (j = syb; j < syb + myb; ++j) {
155: for (i = sxb; i < sxb + mxb; ++i, ++q) {
156: PetscInt r;
158: localPoints[q] = q;
159: PetscFindInt(indices[q], size + 1, ranges, &r);
160: remotePoints[q].rank = r < 0 ? -(r + 1) - 1 : r;
161: remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank];
162: }
163: }
164: }
165: ISRestoreIndices(is, &indices);
166: ISDestroy(&is);
167: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfz);
168: PetscObjectSetName((PetscObject)*sfz, "Buffered Map");
169: PetscSFSetGraph(*sfz, M * N * P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);
170: }
172: VecDestroy(&X);
173: PetscFree2(localPoints, remotePoints);
174: return 0;
175: }
177: typedef enum {
178: PATCH_COMM_TYPE_WORLD = 0,
179: PATCH_COMM_TYPE_SELF = 1
180: } PatchCommType;
182: PetscErrorCode DMPatchSolve(DM dm)
183: {
184: MPI_Comm comm;
185: MPI_Comm commz;
186: DM dmc;
187: PetscSF sfz, sfzr;
188: Vec XC;
189: MatStencil patchSize, commSize, gridRank, lower, upper;
190: PetscInt M, N, P, i, j, k, l, m, n, p = 0;
191: PetscMPIInt rank, size;
192: PetscInt debug = 0;
194: PetscObjectGetComm((PetscObject)dm, &comm);
195: MPI_Comm_rank(comm, &rank);
196: MPI_Comm_size(comm, &size);
197: DMPatchGetCoarse(dm, &dmc);
198: DMPatchGetPatchSize(dm, &patchSize);
199: DMPatchGetCommSize(dm, &commSize);
200: DMPatchGetCommSize(dm, &commSize);
201: DMGetGlobalVector(dmc, &XC);
202: DMDAGetInfo(dmc, NULL, &M, &N, &P, &l, &m, &n, NULL, NULL, NULL, NULL, NULL, NULL);
203: M = PetscMax(M, 1);
204: l = PetscMax(l, 1);
205: N = PetscMax(N, 1);
206: m = PetscMax(m, 1);
207: P = PetscMax(P, 1);
208: n = PetscMax(n, 1);
210: gridRank.i = rank % l;
211: gridRank.j = rank / l % m;
212: gridRank.k = rank / (l * m) % n;
214: if (commSize.i * commSize.j * commSize.k == size || commSize.i * commSize.j * commSize.k == 0) {
215: commSize.i = l;
216: commSize.j = m;
217: commSize.k = n;
218: commz = comm;
219: } else if (commSize.i * commSize.j * commSize.k == 1) {
220: commz = PETSC_COMM_SELF;
221: } else {
222: const PetscMPIInt newComm = ((gridRank.k / commSize.k) * (m / commSize.j) + gridRank.j / commSize.j) * (l / commSize.i) + (gridRank.i / commSize.i);
223: const PetscMPIInt newRank = ((gridRank.k % commSize.k) * commSize.j + gridRank.j % commSize.j) * commSize.i + (gridRank.i % commSize.i);
225: MPI_Comm_split(comm, newComm, newRank, &commz);
226: if (debug) PetscPrintf(PETSC_COMM_SELF, "Rank %d color %d key %d commz %p\n", rank, newComm, newRank, (void *)(MPI_Aint)commz);
227: }
228: /*
229: Assumptions:
230: - patchSize divides gridSize
231: - commSize divides gridSize
232: - commSize divides l,m,n
233: Ignore multiple patches per rank for now
235: Multiple ranks per patch:
236: - l,m,n divides patchSize
237: - commSize divides patchSize
238: */
239: for (k = 0; k < P; k += PetscMax(patchSize.k, 1)) {
240: for (j = 0; j < N; j += PetscMax(patchSize.j, 1)) {
241: for (i = 0; i < M; i += PetscMax(patchSize.i, 1), ++p) {
242: MPI_Comm commp = MPI_COMM_NULL;
243: DM dmz = NULL;
244: #if 0
245: DM dmf = NULL;
246: Mat interpz = NULL;
247: #endif
248: Vec XZ = NULL;
249: PetscScalar *xcarray = NULL;
250: PetscScalar *xzarray = NULL;
252: if ((gridRank.k / commSize.k == p / (l / commSize.i * m / commSize.j) % n / commSize.k) && (gridRank.j / commSize.j == p / (l / commSize.i) % m / commSize.j) && (gridRank.i / commSize.i == p % l / commSize.i)) {
253: if (debug) PetscPrintf(PETSC_COMM_SELF, "Rank %d is accepting Patch %" PetscInt_FMT "\n", rank, p);
254: commp = commz;
255: }
256: /* Zoom to coarse patch */
257: lower.i = i;
258: lower.j = j;
259: lower.k = k;
260: upper.i = i + patchSize.i;
261: upper.j = j + patchSize.j;
262: upper.k = k + patchSize.k;
263: DMPatchZoom(dmc, lower, upper, commp, &dmz, &sfz, &sfzr);
264: lower.c = 0; /* initialize member, otherwise compiler issues warnings */
265: upper.c = 0; /* initialize member, otherwise compiler issues warnings */
266: if (debug)
267: PetscPrintf(comm, "Patch %" PetscInt_FMT ": (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")--(%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", p, lower.i, lower.j, lower.k, upper.i, upper.j, upper.k);
268: if (dmz) DMView(dmz, PETSC_VIEWER_STDOUT_(commz));
269: PetscSFView(sfz, PETSC_VIEWER_STDOUT_(comm));
270: PetscSFView(sfzr, PETSC_VIEWER_STDOUT_(comm));
271: /* Scatter Xcoarse -> Xzoom */
272: if (dmz) DMGetGlobalVector(dmz, &XZ);
273: if (XZ) VecGetArray(XZ, &xzarray);
274: VecGetArray(XC, &xcarray);
275: PetscSFBcastBegin(sfz, MPIU_SCALAR, xcarray, xzarray, MPI_REPLACE);
276: PetscSFBcastEnd(sfz, MPIU_SCALAR, xcarray, xzarray, MPI_REPLACE);
277: VecRestoreArray(XC, &xcarray);
278: if (XZ) VecRestoreArray(XZ, &xzarray);
279: #if 0
280: /* Interpolate Xzoom -> Xfine, note that this may be on subcomms */
281: DMRefine(dmz, MPI_COMM_NULL, &dmf);
282: DMCreateInterpolation(dmz, dmf, &interpz, NULL);
283: DMInterpolate(dmz, interpz, dmf);
284: /* Smooth Xfine using two-step smoother, normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine */
285: /* Compute residual Rfine */
286: /* Restrict Rfine to Rzoom_restricted */
287: #endif
288: /* Scatter Rzoom_restricted -> Rcoarse_restricted */
289: if (XZ) VecGetArray(XZ, &xzarray);
290: VecGetArray(XC, &xcarray);
291: PetscSFReduceBegin(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);
292: PetscSFReduceEnd(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);
293: VecRestoreArray(XC, &xcarray);
294: if (XZ) VecRestoreArray(XZ, &xzarray);
295: if (dmz) DMRestoreGlobalVector(dmz, &XZ);
296: /* Compute global residual Rcoarse */
297: /* TauCoarse = Rcoarse - Rcoarse_restricted */
299: PetscSFDestroy(&sfz);
300: PetscSFDestroy(&sfzr);
301: DMDestroy(&dmz);
302: }
303: }
304: }
305: DMRestoreGlobalVector(dmc, &XC);
306: return 0;
307: }
309: PetscErrorCode DMPatchView_ASCII(DM dm, PetscViewer viewer)
310: {
311: DM_Patch *mesh = (DM_Patch *)dm->data;
312: PetscViewerFormat format;
313: const char *name;
315: PetscViewerGetFormat(viewer, &format);
316: /* if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) */
317: PetscObjectGetName((PetscObject)dm, &name);
318: PetscViewerASCIIPrintf(viewer, "Patch DM %s\n", name);
319: PetscViewerASCIIPushTab(viewer);
320: PetscViewerASCIIPrintf(viewer, "Coarse DM\n");
321: DMView(mesh->dmCoarse, viewer);
322: PetscViewerASCIIPopTab(viewer);
323: return 0;
324: }
326: PetscErrorCode DMView_Patch(DM dm, PetscViewer viewer)
327: {
328: PetscBool iascii, isbinary;
332: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
333: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
334: if (iascii) DMPatchView_ASCII(dm, viewer);
335: return 0;
336: }
338: PetscErrorCode DMDestroy_Patch(DM dm)
339: {
340: DM_Patch *mesh = (DM_Patch *)dm->data;
342: if (--mesh->refct > 0) return 0;
343: DMDestroy(&mesh->dmCoarse);
344: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
345: PetscFree(mesh);
346: return 0;
347: }
349: PetscErrorCode DMSetUp_Patch(DM dm)
350: {
351: DM_Patch *mesh = (DM_Patch *)dm->data;
354: DMSetUp(mesh->dmCoarse);
355: return 0;
356: }
358: PetscErrorCode DMCreateGlobalVector_Patch(DM dm, Vec *g)
359: {
360: DM_Patch *mesh = (DM_Patch *)dm->data;
363: DMCreateGlobalVector(mesh->dmCoarse, g);
364: return 0;
365: }
367: PetscErrorCode DMCreateLocalVector_Patch(DM dm, Vec *l)
368: {
369: DM_Patch *mesh = (DM_Patch *)dm->data;
372: DMCreateLocalVector(mesh->dmCoarse, l);
373: return 0;
374: }
376: PetscErrorCode DMCreateSubDM_Patch(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
377: {
378: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tell me to code this");
379: }
381: PetscErrorCode DMPatchGetCoarse(DM dm, DM *dmCoarse)
382: {
383: DM_Patch *mesh = (DM_Patch *)dm->data;
386: *dmCoarse = mesh->dmCoarse;
387: return 0;
388: }
390: PetscErrorCode DMPatchGetPatchSize(DM dm, MatStencil *patchSize)
391: {
392: DM_Patch *mesh = (DM_Patch *)dm->data;
396: *patchSize = mesh->patchSize;
397: return 0;
398: }
400: PetscErrorCode DMPatchSetPatchSize(DM dm, MatStencil patchSize)
401: {
402: DM_Patch *mesh = (DM_Patch *)dm->data;
405: mesh->patchSize = patchSize;
406: return 0;
407: }
409: PetscErrorCode DMPatchGetCommSize(DM dm, MatStencil *commSize)
410: {
411: DM_Patch *mesh = (DM_Patch *)dm->data;
415: *commSize = mesh->commSize;
416: return 0;
417: }
419: PetscErrorCode DMPatchSetCommSize(DM dm, MatStencil commSize)
420: {
421: DM_Patch *mesh = (DM_Patch *)dm->data;
424: mesh->commSize = commSize;
425: return 0;
426: }