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: }