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,upper - the upper right corner and the lower left corner of the requested patch
29: - commz - the new communicator for the patch, MPI_COMM_NULL indicates that the given rank will not own a patch
31: Output Parameters:
32: + dmz - the patch DM
33: . sfz - the PetscSF mapping the patch+halo to the zoomed version (optional)
34: - sfzr - the PetscSF mapping the patch to the restricted zoomed version
36: Level: intermediate
38: .seealso: DMPatchSolve(), DMDACreatePatchIS()
39: @*/
40: PetscErrorCode DMPatchZoom(DM dm, MatStencil lower, MatStencil upper, MPI_Comm commz, DM *dmz, PetscSF *sfz, PetscSF *sfzr)
41: {
42: DMDAStencilType st;
43: MatStencil blower, bupper, loclower, locupper;
44: IS is;
45: const PetscInt *ranges, *indices;
46: PetscInt *localPoints = NULL;
47: PetscSFNode *remotePoints = NULL;
48: PetscInt dim, dof;
49: 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;
50: PetscMPIInt size;
51: PetscBool patchis_offproc = PETSC_TRUE;
52: PetscErrorCode ierr;
53: Vec X;
56: if (!sfz) halo = 0;
57: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
58: /* Create patch DM */
59: DMDAGetInfo(dm, &dim, &M, &N, &P, NULL,NULL,NULL, &dof, NULL,NULL,NULL,NULL, &st);
61: /* Get piece for rank r, expanded by halo */
62: bupper.i = PetscMin(M, upper.i + halo); blower.i = PetscMax(lower.i - halo, 0);
63: bupper.j = PetscMin(N, upper.j + halo); blower.j = PetscMax(lower.j - halo, 0);
64: bupper.k = PetscMin(P, upper.k + halo); blower.k = PetscMax(lower.k - halo, 0);
65: rM = bupper.i - blower.i;
66: rN = bupper.j - blower.j;
67: rP = bupper.k - blower.k;
69: if (commz != MPI_COMM_NULL) {
70: DMDACreate(commz, dmz);
71: DMSetDimension(*dmz, dim);
72: DMDASetSizes(*dmz, rM, rN, rP);
73: DMDASetNumProcs(*dmz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
74: DMDASetBoundaryType(*dmz, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE);
75: DMDASetDof(*dmz, dof);
76: DMDASetStencilType(*dmz, st);
77: DMDASetStencilWidth(*dmz, 0);
78: DMDASetOwnershipRanges(*dmz, NULL, NULL, NULL);
79: DMSetFromOptions(*dmz);
80: DMSetUp(*dmz);
81: DMDAGetCorners(*dmz, &sxb, &syb, &szb, &mxb, &myb, &mzb);
82: sxr = PetscMax(sxb, lower.i - blower.i);
83: syr = PetscMax(syb, lower.j - blower.j);
84: szr = PetscMax(szb, lower.k - blower.k);
85: exr = PetscMin(sxb+mxb, upper.i - blower.i);
86: eyr = PetscMin(syb+myb, upper.j - blower.j);
87: ezr = PetscMin(szb+mzb, upper.k - blower.k);
88: PetscMalloc2(dof*rM*rN*PetscMax(rP,1),&localPoints,dof*rM*rN*PetscMax(rP,1),&remotePoints);
89: } else {
90: sxr = syr = szr = exr = eyr = ezr = sxb = syb = szb = mxb = myb = mzb = 0;
91: }
93: /* Create SF for restricted map */
94: DMCreateGlobalVector(dm,&X);
95: VecGetOwnershipRanges(X,&ranges);
97: loclower.i = blower.i + sxr; locupper.i = blower.i + exr;
98: loclower.j = blower.j + syr; locupper.j = blower.j + eyr;
99: loclower.k = blower.k + szr; locupper.k = blower.k + ezr;
101: DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc);
102: ISGetIndices(is, &indices);
104: if (dim < 3) {mzb = 1; ezr = 1;}
105: q = 0;
106: for (k = szb; k < szb+mzb; ++k) {
107: if ((k < szr) || (k >= ezr)) continue;
108: for (j = syb; j < syb+myb; ++j) {
109: if ((j < syr) || (j >= eyr)) continue;
110: for (i = sxb; i < sxb+mxb; ++i) {
111: for (l=0; l<dof; l++) {
112: const PetscInt lp = l + dof*(((k-szb)*rN + (j-syb))*rM + i-sxb);
113: PetscInt r;
115: if ((i < sxr) || (i >= exr)) continue;
116: localPoints[q] = lp;
117: PetscFindInt(indices[q], size+1, ranges, &r);
119: remotePoints[q].rank = r < 0 ? -(r+1) - 1 : r;
120: remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank];
121: ++q;
122: }
123: }
124: }
125: }
126: ISRestoreIndices(is, &indices);
127: ISDestroy(&is);
128: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfzr);
129: PetscObjectSetName((PetscObject) *sfzr, "Restricted Map");
130: PetscSFSetGraph(*sfzr, dof*M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);
132: if (sfz) {
133: /* Create SF for buffered map */
134: loclower.i = blower.i + sxb; locupper.i = blower.i + sxb+mxb;
135: loclower.j = blower.j + syb; locupper.j = blower.j + syb+myb;
136: loclower.k = blower.k + szb; locupper.k = blower.k + szb+mzb;
138: DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc);
139: ISGetIndices(is, &indices);
141: q = 0;
142: for (k = szb; k < szb+mzb; ++k) {
143: for (j = syb; j < syb+myb; ++j) {
144: for (i = sxb; i < sxb+mxb; ++i, ++q) {
145: PetscInt r;
147: localPoints[q] = q;
148: PetscFindInt(indices[q], size+1, ranges, &r);
149: remotePoints[q].rank = r < 0 ? -(r+1) - 1 : r;
150: remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank];
151: }
152: }
153: }
154: ISRestoreIndices(is, &indices);
155: ISDestroy(&is);
156: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfz);
157: PetscObjectSetName((PetscObject) *sfz, "Buffered Map");
158: PetscSFSetGraph(*sfz, M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);
159: }
161: VecDestroy(&X);
162: PetscFree2(localPoints, remotePoints);
163: return(0);
164: }
166: typedef enum {PATCH_COMM_TYPE_WORLD = 0, PATCH_COMM_TYPE_SELF = 1} PatchCommType;
168: PetscErrorCode DMPatchSolve(DM dm)
169: {
170: MPI_Comm comm;
171: MPI_Comm commz;
172: DM dmc;
173: PetscSF sfz, sfzr;
174: Vec XC;
175: MatStencil patchSize, commSize, gridRank, lower, upper;
176: PetscInt M, N, P, i, j, k, l, m, n, p = 0;
177: PetscMPIInt rank, size;
178: PetscInt debug = 0;
182: PetscObjectGetComm((PetscObject)dm,&comm);
183: MPI_Comm_rank(comm, &rank);
184: MPI_Comm_size(comm, &size);
185: DMPatchGetCoarse(dm, &dmc);
186: DMPatchGetPatchSize(dm, &patchSize);
187: DMPatchGetCommSize(dm, &commSize);
188: DMPatchGetCommSize(dm, &commSize);
189: DMGetGlobalVector(dmc, &XC);
190: DMDAGetInfo(dmc, NULL, &M, &N, &P, &l, &m, &n, NULL,NULL,NULL,NULL,NULL,NULL);
191: M = PetscMax(M, 1); l = PetscMax(l, 1);
192: N = PetscMax(N, 1); m = PetscMax(m, 1);
193: P = PetscMax(P, 1); n = PetscMax(n, 1);
195: gridRank.i = rank % l;
196: gridRank.j = rank/l % m;
197: gridRank.k = rank/(l*m) % n;
199: if (commSize.i*commSize.j*commSize.k == size || commSize.i*commSize.j*commSize.k == 0) {
200: commSize.i = l; commSize.j = m; commSize.k = n;
201: commz = comm;
202: } else if (commSize.i*commSize.j*commSize.k == 1) {
203: commz = PETSC_COMM_SELF;
204: } else {
205: const PetscMPIInt newComm = ((gridRank.k/commSize.k)*(m/commSize.j) + gridRank.j/commSize.j)*(l/commSize.i) + (gridRank.i/commSize.i);
206: const PetscMPIInt newRank = ((gridRank.k%commSize.k)*commSize.j + gridRank.j%commSize.j)*commSize.i + (gridRank.i%commSize.i);
208: MPI_Comm_split(comm, newComm, newRank, &commz);
209: if (debug) {PetscPrintf(PETSC_COMM_SELF, "Rank %d color %d key %d commz %d\n", rank, newComm, newRank, *((PetscMPIInt*) &commz));}
210: }
211: /*
212: Assumptions:
213: - patchSize divides gridSize
214: - commSize divides gridSize
215: - commSize divides l,m,n
216: Ignore multiple patches per rank for now
218: Multiple ranks per patch:
219: - l,m,n divides patchSize
220: - commSize divides patchSize
221: */
222: for (k = 0; k < P; k += PetscMax(patchSize.k, 1)) {
223: for (j = 0; j < N; j += PetscMax(patchSize.j, 1)) {
224: for (i = 0; i < M; i += PetscMax(patchSize.i, 1), ++p) {
225: MPI_Comm commp = MPI_COMM_NULL;
226: DM dmz = NULL;
227: #if 0
228: DM dmf = NULL;
229: Mat interpz = NULL;
230: #endif
231: Vec XZ = NULL;
232: PetscScalar *xcarray = NULL;
233: PetscScalar *xzarray = NULL;
235: if ((gridRank.k/commSize.k == p/(l/commSize.i * m/commSize.j) % n/commSize.k) &&
236: (gridRank.j/commSize.j == p/(l/commSize.i) % m/commSize.j) &&
237: (gridRank.i/commSize.i == p % l/commSize.i)) {
238: if (debug) {PetscPrintf(PETSC_COMM_SELF, "Rank %d is accepting Patch %d\n", rank, p);}
239: commp = commz;
240: }
241: /* Zoom to coarse patch */
242: lower.i = i; lower.j = j; lower.k = k;
243: upper.i = i + patchSize.i; upper.j = j + patchSize.j; upper.k = k + patchSize.k;
244: DMPatchZoom(dmc, lower, upper, commp, &dmz, &sfz, &sfzr);
245: lower.c = 0; /* initialize member, otherwise compiler issues warnings */
246: upper.c = 0; /* initialize member, otherwise compiler issues warnings */
247: if (debug) {PetscPrintf(comm, "Patch %d: (%d, %d, %d)--(%d, %d, %d)\n", p, lower.i, lower.j, lower.k, upper.i, upper.j, upper.k);}
248: if (dmz) {DMView(dmz, PETSC_VIEWER_STDOUT_(commz));}
249: PetscSFView(sfz, PETSC_VIEWER_STDOUT_(comm));
250: PetscSFView(sfzr, PETSC_VIEWER_STDOUT_(comm));
251: /* Scatter Xcoarse -> Xzoom */
252: if (dmz) {DMGetGlobalVector(dmz, &XZ);}
253: if (XZ) {VecGetArray(XZ, &xzarray);}
254: VecGetArray(XC, &xcarray);
255: PetscSFBcastBegin(sfz, MPIU_SCALAR, xcarray, xzarray,MPI_REPLACE);
256: PetscSFBcastEnd(sfz, MPIU_SCALAR, xcarray, xzarray,MPI_REPLACE);
257: VecRestoreArray(XC, &xcarray);
258: if (XZ) {VecRestoreArray(XZ, &xzarray);}
259: #if 0
260: /* Interpolate Xzoom -> Xfine, note that this may be on subcomms */
261: DMRefine(dmz, MPI_COMM_NULL, &dmf);
262: DMCreateInterpolation(dmz, dmf, &interpz, NULL);
263: DMInterpolate(dmz, interpz, dmf);
264: /* Smooth Xfine using two-step smoother, normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine */
265: /* Compute residual Rfine */
266: /* Restrict Rfine to Rzoom_restricted */
267: #endif
268: /* Scatter Rzoom_restricted -> Rcoarse_restricted */
269: if (XZ) {VecGetArray(XZ, &xzarray);}
270: VecGetArray(XC, &xcarray);
271: PetscSFReduceBegin(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);
272: PetscSFReduceEnd(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);
273: VecRestoreArray(XC, &xcarray);
274: if (XZ) {VecRestoreArray(XZ, &xzarray);}
275: if (dmz) {DMRestoreGlobalVector(dmz, &XZ);}
276: /* Compute global residual Rcoarse */
277: /* TauCoarse = Rcoarse - Rcoarse_restricted */
279: PetscSFDestroy(&sfz);
280: PetscSFDestroy(&sfzr);
281: DMDestroy(&dmz);
282: }
283: }
284: }
285: DMRestoreGlobalVector(dmc, &XC);
286: return(0);
287: }
289: PetscErrorCode DMPatchView_ASCII(DM dm, PetscViewer viewer)
290: {
291: DM_Patch *mesh = (DM_Patch*) dm->data;
292: PetscViewerFormat format;
293: const char *name;
294: PetscErrorCode ierr;
297: PetscViewerGetFormat(viewer, &format);
298: /* if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) */
299: PetscObjectGetName((PetscObject) dm, &name);
300: PetscViewerASCIIPrintf(viewer, "Patch DM %s\n", name);
301: PetscViewerASCIIPushTab(viewer);
302: PetscViewerASCIIPrintf(viewer, "Coarse DM\n");
303: DMView(mesh->dmCoarse, viewer);
304: PetscViewerASCIIPopTab(viewer);
305: return(0);
306: }
308: PetscErrorCode DMView_Patch(DM dm, PetscViewer viewer)
309: {
310: PetscBool iascii, isbinary;
316: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
317: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
318: if (iascii) {
319: DMPatchView_ASCII(dm, viewer);
320: }
321: return(0);
322: }
324: PetscErrorCode DMDestroy_Patch(DM dm)
325: {
326: DM_Patch *mesh = (DM_Patch*) dm->data;
330: if (--mesh->refct > 0) return(0);
331: DMDestroy(&mesh->dmCoarse);
332: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
333: PetscFree(mesh);
334: return(0);
335: }
337: PetscErrorCode DMSetUp_Patch(DM dm)
338: {
339: DM_Patch *mesh = (DM_Patch*) dm->data;
344: DMSetUp(mesh->dmCoarse);
345: return(0);
346: }
348: PetscErrorCode DMCreateGlobalVector_Patch(DM dm, Vec *g)
349: {
350: DM_Patch *mesh = (DM_Patch*) dm->data;
355: DMCreateGlobalVector(mesh->dmCoarse, g);
356: return(0);
357: }
359: PetscErrorCode DMCreateLocalVector_Patch(DM dm, Vec *l)
360: {
361: DM_Patch *mesh = (DM_Patch*) dm->data;
366: DMCreateLocalVector(mesh->dmCoarse, l);
367: return(0);
368: }
370: PetscErrorCode DMCreateSubDM_Patch(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
371: {
372: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tell me to code this");
373: }
375: PetscErrorCode DMPatchGetCoarse(DM dm, DM *dmCoarse)
376: {
377: DM_Patch *mesh = (DM_Patch*) dm->data;
381: *dmCoarse = mesh->dmCoarse;
382: return(0);
383: }
385: PetscErrorCode DMPatchGetPatchSize(DM dm, MatStencil *patchSize)
386: {
387: DM_Patch *mesh = (DM_Patch*) dm->data;
392: *patchSize = mesh->patchSize;
393: return(0);
394: }
396: PetscErrorCode DMPatchSetPatchSize(DM dm, MatStencil patchSize)
397: {
398: DM_Patch *mesh = (DM_Patch*) dm->data;
402: mesh->patchSize = patchSize;
403: return(0);
404: }
406: PetscErrorCode DMPatchGetCommSize(DM dm, MatStencil *commSize)
407: {
408: DM_Patch *mesh = (DM_Patch*) dm->data;
413: *commSize = mesh->commSize;
414: return(0);
415: }
417: PetscErrorCode DMPatchSetCommSize(DM dm, MatStencil commSize)
418: {
419: DM_Patch *mesh = (DM_Patch*) dm->data;
423: mesh->commSize = commSize;
424: return(0);
425: }