Actual source code: patch.c

petsc-master 2019-09-17
Report Typos and Errors
  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, 0,0,0, &dof, 0,0,0,0, &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, 0, &M, &N, &P, &l, &m, &n, 0,0,0,0,0,0);
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: }