Actual source code: patch.c
petsc-3.14.6 2021-03-30
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: /*
22: DMPatchZoom - Create a version of the coarse patch (identified by rank) with halo on communicator commz
24: Collective on dm
26: Input Parameters:
27: + dm - the DM
28: . rank - the rank which holds the given patch
29: - commz - the new communicator for the patch
31: Output Parameters:
32: + dmz - the patch DM
33: . sfz - the PetscSF mapping the patch+halo to the zoomed version
34: . sfzr - the PetscSF mapping the patch to the restricted zoomed version
36: Level: intermediate
38: Note: All processes in commz should have the same rank (could autosplit comm)
40: .seealso: DMPatchSolve()
41: */
42: PetscErrorCode DMPatchZoom(DM dm, Vec X, MatStencil lower, MatStencil upper, MPI_Comm commz, DM *dmz, PetscSF *sfz, PetscSF *sfzr)
43: {
44: DMDAStencilType st;
45: MatStencil blower, bupper, loclower, locupper;
46: IS is;
47: const PetscInt *ranges, *indices;
48: PetscInt *localPoints = NULL;
49: PetscSFNode *remotePoints = NULL;
50: PetscInt dim, dof;
51: PetscInt M, N, P, rM, rN, rP, halo = 1, sxb, syb, szb, sxr, syr, szr, exr, eyr, ezr, mxb, myb, mzb, i, j, k, q;
52: PetscMPIInt size;
53: PetscErrorCode ierr;
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); blower.i = PetscMax(lower.i - halo, 0);
62: bupper.j = PetscMin(N, upper.j + halo); blower.j = PetscMax(lower.j - halo, 0);
63: bupper.k = PetscMin(P, upper.k + halo); blower.k = PetscMax(lower.k - halo, 0);
64: rM = bupper.i - blower.i;
65: rN = bupper.j - blower.j;
66: rP = bupper.k - blower.k;
68: if (commz != MPI_COMM_NULL) {
69: DMDACreate(commz, dmz);
70: DMSetDimension(*dmz, dim);
71: DMDASetSizes(*dmz, rM, rN, rP);
72: DMDASetNumProcs(*dmz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
73: DMDASetBoundaryType(*dmz, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE);
74: DMDASetDof(*dmz, dof);
75: DMDASetStencilType(*dmz, st);
76: DMDASetStencilWidth(*dmz, 0);
77: DMDASetOwnershipRanges(*dmz, NULL, NULL, NULL);
78: DMSetFromOptions(*dmz);
79: DMSetUp(*dmz);
80: DMDAGetCorners(*dmz, &sxb, &syb, &szb, &mxb, &myb, &mzb);
81: sxr = PetscMax(sxb, lower.i - blower.i);
82: syr = PetscMax(syb, lower.j - blower.j);
83: szr = PetscMax(szb, lower.k - blower.k);
84: exr = PetscMin(sxb+mxb, upper.i - blower.i);
85: eyr = PetscMin(syb+myb, upper.j - blower.j);
86: ezr = PetscMin(szb+mzb, upper.k - blower.k);
87: PetscMalloc2(rM*rN*rP,&localPoints,rM*rN*rP,&remotePoints);
88: } else {
89: sxr = syr = szr = exr = eyr = ezr = sxb = syb = szb = mxb = myb = mzb = 0;
90: }
92: /* Create SF for restricted map */
93: VecGetOwnershipRanges(X,&ranges);
95: loclower.i = blower.i + sxr; locupper.i = blower.i + exr;
96: loclower.j = blower.j + syr; locupper.j = blower.j + eyr;
97: loclower.k = blower.k + szr; locupper.k = blower.k + ezr;
99: DMDACreatePatchIS(dm, &loclower, &locupper, &is);
100: ISGetIndices(is, &indices);
102: q = 0;
103: for (k = szb; k < szb+mzb; ++k) {
104: if ((k < szr) || (k >= ezr)) continue;
105: for (j = syb; j < syb+myb; ++j) {
106: if ((j < syr) || (j >= eyr)) continue;
107: for (i = sxb; i < sxb+mxb; ++i) {
108: const PetscInt lp = ((k-szb)*rN + (j-syb))*rM + i-sxb;
109: PetscInt r;
111: if ((i < sxr) || (i >= exr)) continue;
112: localPoints[q] = lp;
113: PetscFindInt(indices[q], size+1, ranges, &r);
115: remotePoints[q].rank = r < 0 ? -(r+1) - 1 : r;
116: remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank];
117: ++q;
118: }
119: }
120: }
121: ISRestoreIndices(is, &indices);
122: ISDestroy(&is);
123: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfzr);
124: PetscObjectSetName((PetscObject) *sfzr, "Restricted Map");
125: PetscSFSetGraph(*sfzr, M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);
127: /* Create SF for buffered map */
128: loclower.i = blower.i + sxb; locupper.i = blower.i + sxb+mxb;
129: loclower.j = blower.j + syb; locupper.j = blower.j + syb+myb;
130: loclower.k = blower.k + szb; locupper.k = blower.k + szb+mzb;
132: DMDACreatePatchIS(dm, &loclower, &locupper, &is);
133: ISGetIndices(is, &indices);
135: q = 0;
136: for (k = szb; k < szb+mzb; ++k) {
137: for (j = syb; j < syb+myb; ++j) {
138: for (i = sxb; i < sxb+mxb; ++i, ++q) {
139: PetscInt r;
141: localPoints[q] = q;
142: PetscFindInt(indices[q], size+1, ranges, &r);
143: remotePoints[q].rank = r < 0 ? -(r+1) - 1 : r;
144: remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank];
145: }
146: }
147: }
148: ISRestoreIndices(is, &indices);
149: ISDestroy(&is);
150: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfz);
151: PetscObjectSetName((PetscObject) *sfz, "Buffered Map");
152: PetscSFSetGraph(*sfz, M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);
154: PetscFree2(localPoints, remotePoints);
155: return(0);
156: }
158: typedef enum {PATCH_COMM_TYPE_WORLD = 0, PATCH_COMM_TYPE_SELF = 1} PatchCommType;
160: PetscErrorCode DMPatchSolve(DM dm)
161: {
162: MPI_Comm comm;
163: MPI_Comm commz;
164: DM dmc;
165: PetscSF sfz, sfzr;
166: Vec XC;
167: MatStencil patchSize, commSize, gridRank, lower, upper;
168: PetscInt M, N, P, i, j, k, l, m, n, p = 0;
169: PetscMPIInt rank, size;
170: PetscInt debug = 0;
174: PetscObjectGetComm((PetscObject)dm,&comm);
175: MPI_Comm_rank(comm, &rank);
176: MPI_Comm_size(comm, &size);
177: DMPatchGetCoarse(dm, &dmc);
178: DMPatchGetPatchSize(dm, &patchSize);
179: DMPatchGetCommSize(dm, &commSize);
180: DMPatchGetCommSize(dm, &commSize);
181: DMGetGlobalVector(dmc, &XC);
182: DMDAGetInfo(dmc, NULL, &M, &N, &P, &l, &m, &n, NULL,NULL,NULL,NULL,NULL,NULL);
183: M = PetscMax(M, 1); l = PetscMax(l, 1);
184: N = PetscMax(N, 1); m = PetscMax(m, 1);
185: P = PetscMax(P, 1); n = PetscMax(n, 1);
187: gridRank.i = rank % l;
188: gridRank.j = rank/l % m;
189: gridRank.k = rank/(l*m) % n;
191: if (commSize.i*commSize.j*commSize.k == size || commSize.i*commSize.j*commSize.k == 0) {
192: commSize.i = l; commSize.j = m; commSize.k = n;
193: commz = comm;
194: } else if (commSize.i*commSize.j*commSize.k == 1) {
195: commz = PETSC_COMM_SELF;
196: } else {
197: const PetscMPIInt newComm = ((gridRank.k/commSize.k)*(m/commSize.j) + gridRank.j/commSize.j)*(l/commSize.i) + (gridRank.i/commSize.i);
198: const PetscMPIInt newRank = ((gridRank.k%commSize.k)*commSize.j + gridRank.j%commSize.j)*commSize.i + (gridRank.i%commSize.i);
200: MPI_Comm_split(comm, newComm, newRank, &commz);
201: if (debug) {PetscPrintf(PETSC_COMM_SELF, "Rank %d color %d key %d commz %d\n", rank, newComm, newRank, *((PetscMPIInt*) &commz));}
202: }
203: /*
204: Assumptions:
205: - patchSize divides gridSize
206: - commSize divides gridSize
207: - commSize divides l,m,n
208: Ignore multiple patches per rank for now
210: Multiple ranks per patch:
211: - l,m,n divides patchSize
212: - commSize divides patchSize
213: */
214: for (k = 0; k < P; k += PetscMax(patchSize.k, 1)) {
215: for (j = 0; j < N; j += PetscMax(patchSize.j, 1)) {
216: for (i = 0; i < M; i += PetscMax(patchSize.i, 1), ++p) {
217: MPI_Comm commp = MPI_COMM_NULL;
218: DM dmz = NULL;
219: #if 0
220: DM dmf = NULL;
221: Mat interpz = NULL;
222: #endif
223: Vec XZ = NULL;
224: PetscScalar *xcarray = NULL;
225: PetscScalar *xzarray = NULL;
227: if ((gridRank.k/commSize.k == p/(l/commSize.i * m/commSize.j) % n/commSize.k) &&
228: (gridRank.j/commSize.j == p/(l/commSize.i) % m/commSize.j) &&
229: (gridRank.i/commSize.i == p % l/commSize.i)) {
230: if (debug) {PetscPrintf(PETSC_COMM_SELF, "Rank %d is accepting Patch %d\n", rank, p);}
231: commp = commz;
232: }
233: /* Zoom to coarse patch */
234: lower.i = i; lower.j = j; lower.k = k;
235: upper.i = i + patchSize.i; upper.j = j + patchSize.j; upper.k = k + patchSize.k;
236: DMPatchZoom(dmc, XC, lower, upper, commp, &dmz, &sfz, &sfzr);
237: lower.c = 0; /* initialize member, otherwise compiler issues warnings */
238: upper.c = 0; /* initialize member, otherwise compiler issues warnings */
239: /* Debug */
240: PetscPrintf(comm, "Patch %d: (%d, %d, %d)--(%d, %d, %d)\n", p, lower.i, lower.j, lower.k, upper.i, upper.j, upper.k);
241: if (dmz) {DMView(dmz, PETSC_VIEWER_STDOUT_(commz));}
242: PetscSFView(sfz, PETSC_VIEWER_STDOUT_(comm));
243: PetscSFView(sfzr, PETSC_VIEWER_STDOUT_(comm));
244: /* Scatter Xcoarse -> Xzoom */
245: if (dmz) {DMGetGlobalVector(dmz, &XZ);}
246: if (XZ) {VecGetArray(XZ, &xzarray);}
247: VecGetArray(XC, &xcarray);
248: PetscSFBcastBegin(sfz, MPIU_SCALAR, xcarray, xzarray);
249: PetscSFBcastEnd(sfz, MPIU_SCALAR, xcarray, xzarray);
250: VecRestoreArray(XC, &xcarray);
251: if (XZ) {VecRestoreArray(XZ, &xzarray);}
252: #if 0
253: /* Interpolate Xzoom -> Xfine, note that this may be on subcomms */
254: DMRefine(dmz, MPI_COMM_NULL, &dmf);
255: DMCreateInterpolation(dmz, dmf, &interpz, NULL);
256: DMInterpolate(dmz, interpz, dmf);
257: /* Smooth Xfine using two-step smoother, normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine */
258: /* Compute residual Rfine */
259: /* Restrict Rfine to Rzoom_restricted */
260: #endif
261: /* Scatter Rzoom_restricted -> Rcoarse_restricted */
262: if (XZ) {VecGetArray(XZ, &xzarray);}
263: VecGetArray(XC, &xcarray);
264: PetscSFReduceBegin(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);
265: PetscSFReduceEnd(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);
266: VecRestoreArray(XC, &xcarray);
267: if (XZ) {VecRestoreArray(XZ, &xzarray);}
268: if (dmz) {DMRestoreGlobalVector(dmz, &XZ);}
269: /* Compute global residual Rcoarse */
270: /* TauCoarse = Rcoarse - Rcoarse_restricted */
272: PetscSFDestroy(&sfz);
273: PetscSFDestroy(&sfzr);
274: DMDestroy(&dmz);
275: }
276: }
277: }
278: DMRestoreGlobalVector(dmc, &XC);
279: return(0);
280: }
282: PetscErrorCode DMPatchView_Ascii(DM dm, PetscViewer viewer)
283: {
284: DM_Patch *mesh = (DM_Patch*) dm->data;
285: PetscViewerFormat format;
286: const char *name;
287: PetscErrorCode ierr;
290: PetscViewerGetFormat(viewer, &format);
291: /* if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) */
292: PetscObjectGetName((PetscObject) dm, &name);
293: PetscViewerASCIIPrintf(viewer, "Patch DM %s\n", name);
294: PetscViewerASCIIPushTab(viewer);
295: PetscViewerASCIIPrintf(viewer, "Coarse DM\n");
296: DMView(mesh->dmCoarse, viewer);
297: PetscViewerASCIIPopTab(viewer);
298: return(0);
299: }
301: PetscErrorCode DMView_Patch(DM dm, PetscViewer viewer)
302: {
303: PetscBool iascii, isbinary;
309: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
310: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
311: if (iascii) {
312: DMPatchView_Ascii(dm, viewer);
313: #if 0
314: } else if (isbinary) {
315: DMPatchView_Binary(dm, viewer);
316: #endif
317: }
318: return(0);
319: }
321: PetscErrorCode DMDestroy_Patch(DM dm)
322: {
323: DM_Patch *mesh = (DM_Patch*) dm->data;
327: if (--mesh->refct > 0) return(0);
328: DMDestroy(&mesh->dmCoarse);
329: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
330: PetscFree(mesh);
331: return(0);
332: }
334: PetscErrorCode DMSetUp_Patch(DM dm)
335: {
336: DM_Patch *mesh = (DM_Patch*) dm->data;
341: DMSetUp(mesh->dmCoarse);
342: return(0);
343: }
345: PetscErrorCode DMCreateGlobalVector_Patch(DM dm, Vec *g)
346: {
347: DM_Patch *mesh = (DM_Patch*) dm->data;
352: DMCreateGlobalVector(mesh->dmCoarse, g);
353: return(0);
354: }
356: PetscErrorCode DMCreateLocalVector_Patch(DM dm, Vec *l)
357: {
358: DM_Patch *mesh = (DM_Patch*) dm->data;
363: DMCreateLocalVector(mesh->dmCoarse, l);
364: return(0);
365: }
367: PetscErrorCode DMCreateSubDM_Patch(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
368: {
369: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tell me to code this");
370: }
372: PetscErrorCode DMPatchGetCoarse(DM dm, DM *dmCoarse)
373: {
374: DM_Patch *mesh = (DM_Patch*) dm->data;
378: *dmCoarse = mesh->dmCoarse;
379: return(0);
380: }
382: PetscErrorCode DMPatchGetPatchSize(DM dm, MatStencil *patchSize)
383: {
384: DM_Patch *mesh = (DM_Patch*) dm->data;
389: *patchSize = mesh->patchSize;
390: return(0);
391: }
393: PetscErrorCode DMPatchSetPatchSize(DM dm, MatStencil patchSize)
394: {
395: DM_Patch *mesh = (DM_Patch*) dm->data;
399: mesh->patchSize = patchSize;
400: return(0);
401: }
403: PetscErrorCode DMPatchGetCommSize(DM dm, MatStencil *commSize)
404: {
405: DM_Patch *mesh = (DM_Patch*) dm->data;
410: *commSize = mesh->commSize;
411: return(0);
412: }
414: PetscErrorCode DMPatchSetCommSize(DM dm, MatStencil commSize)
415: {
416: DM_Patch *mesh = (DM_Patch*) dm->data;
420: mesh->commSize = commSize;
421: return(0);
422: }