Actual source code: dadd.c

  1: #include <petsc/private/dmdaimpl.h>

  3: /*@
  4:   DMDACreatePatchIS - Creates an index set corresponding to a patch of the `DMDA`.

  6:   Collective

  8:   Input Parameters:
  9: +  da - the `DMDA`
 10: .  lower - a matstencil with i, j and k corresponding to the lower corner of the patch
 11: .  upper - a matstencil with i, j and k corresponding to the upper corner of the patch
 12: -  offproc - indicate whether the returned IS will contain off process indices

 14:   Output Parameter:
 15: .  is - the `IS` corresponding to the patch

 17:   Level: developer

 19:   Notes:
 20:   This routine always returns an `IS` on the `DMDA` comm, if offproc is set to `PETSC_TRUE`,
 21:   the routine returns an `IS` with all the indices requested regardless of whether these indices
 22:   are present on the requesting rank or not. Thus, it is upon the caller to ensure that
 23:   the indices returned in this mode are appropriate. If offproc is set to `PETSC_FALSE`,
 24:   the `IS` only returns the subset of indices that are present on the requesting rank and there
 25:   is no duplication of indices.

 27: .seealso: `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMCreateDomainDecompositionScatters()`
 28: @*/
 29: PetscErrorCode DMDACreatePatchIS(DM da, MatStencil *lower, MatStencil *upper, IS *is, PetscBool offproc)
 30: {
 31:   PetscInt        ms = 0, ns = 0, ps = 0;
 32:   PetscInt        mw = 0, nw = 0, pw = 0;
 33:   PetscInt        me = 1, ne = 1, pe = 1;
 34:   PetscInt        mr = 0, nr = 0, pr = 0;
 35:   PetscInt        ii, jj, kk;
 36:   PetscInt        si, sj, sk;
 37:   PetscInt        i, j, k, l, idx = 0;
 38:   PetscInt        base;
 39:   PetscInt        xm = 1, ym = 1, zm = 1;
 40:   PetscInt        ox, oy, oz;
 41:   PetscInt        m, n, p, M, N, P, dof;
 42:   const PetscInt *lx, *ly, *lz;
 43:   PetscInt        nindices;
 44:   PetscInt       *indices;
 45:   DM_DA          *dd     = (DM_DA *)da->data;
 46:   PetscBool       skip_i = PETSC_TRUE, skip_j = PETSC_TRUE, skip_k = PETSC_TRUE;
 47:   PetscBool       valid_j = PETSC_FALSE, valid_k = PETSC_FALSE; /* DMDA has at least 1 dimension */

 49:   PetscFunctionBegin;
 50:   M   = dd->M;
 51:   N   = dd->N;
 52:   P   = dd->P;
 53:   m   = dd->m;
 54:   n   = dd->n;
 55:   p   = dd->p;
 56:   dof = dd->w;

 58:   nindices = -1;
 59:   if (PetscLikely(upper->i - lower->i)) {
 60:     nindices = nindices * (upper->i - lower->i);
 61:     skip_i   = PETSC_FALSE;
 62:   }
 63:   if (N > 1) {
 64:     valid_j = PETSC_TRUE;
 65:     if (PetscLikely(upper->j - lower->j)) {
 66:       nindices = nindices * (upper->j - lower->j);
 67:       skip_j   = PETSC_FALSE;
 68:     }
 69:   }
 70:   if (P > 1) {
 71:     valid_k = PETSC_TRUE;
 72:     if (PetscLikely(upper->k - lower->k)) {
 73:       nindices = nindices * (upper->k - lower->k);
 74:       skip_k   = PETSC_FALSE;
 75:     }
 76:   }
 77:   if (PetscLikely(nindices < 0)) {
 78:     if (PetscUnlikely(skip_i && skip_j && skip_k)) {
 79:       nindices = 0;
 80:     } else nindices = nindices * (-1);
 81:   } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Lower and Upper stencils are identical! Please check inputs.");

 83:   PetscCall(PetscMalloc1(nindices * dof, &indices));
 84:   PetscCall(DMDAGetOffset(da, &ox, &oy, &oz, NULL, NULL, NULL));

 86:   if (!valid_k) {
 87:     k        = 0;
 88:     upper->k = 0;
 89:     lower->k = 0;
 90:   }
 91:   if (!valid_j) {
 92:     j        = 0;
 93:     upper->j = 0;
 94:     lower->j = 0;
 95:   }

 97:   if (offproc) {
 98:     PetscCall(DMDAGetOwnershipRanges(da, &lx, &ly, &lz));
 99:     /* start at index 0 on processor 0 */
100:     mr = 0;
101:     nr = 0;
102:     pr = 0;
103:     ms = 0;
104:     ns = 0;
105:     ps = 0;
106:     if (lx) me = lx[0];
107:     if (ly) ne = ly[0];
108:     if (lz) pe = lz[0];
109:     /*
110:        If no indices are to be returned, create an empty is,
111:        this prevents hanging in while loops
112:     */
113:     if (skip_i && skip_j && skip_k) goto createis;
114:     /*
115:        do..while loops to ensure the block gets entered once,
116:        regardless of control condition being met, necessary for
117:        cases when a subset of skip_i/j/k is true
118:     */
119:     if (skip_k) k = upper->k - oz;
120:     else k = lower->k - oz;
121:     do {
122:       if (skip_j) j = upper->j - oy;
123:       else j = lower->j - oy;
124:       do {
125:         if (skip_i) i = upper->i - ox;
126:         else i = lower->i - ox;
127:         do {
128:           /* "actual" indices rather than ones outside of the domain */
129:           ii = i;
130:           jj = j;
131:           kk = k;
132:           if (ii < 0) ii = M + ii;
133:           if (jj < 0) jj = N + jj;
134:           if (kk < 0) kk = P + kk;
135:           if (ii > M - 1) ii = ii - M;
136:           if (jj > N - 1) jj = jj - N;
137:           if (kk > P - 1) kk = kk - P;
138:           /* gone out of processor range on x axis */
139:           while (ii > me - 1 || ii < ms) {
140:             if (mr == m - 1) {
141:               ms = 0;
142:               me = lx[0];
143:               mr = 0;
144:             } else {
145:               mr++;
146:               ms = me;
147:               me += lx[mr];
148:             }
149:           }
150:           /* gone out of processor range on y axis */
151:           while (jj > ne - 1 || jj < ns) {
152:             if (nr == n - 1) {
153:               ns = 0;
154:               ne = ly[0];
155:               nr = 0;
156:             } else {
157:               nr++;
158:               ns = ne;
159:               ne += ly[nr];
160:             }
161:           }
162:           /* gone out of processor range on z axis */
163:           while (kk > pe - 1 || kk < ps) {
164:             if (pr == p - 1) {
165:               ps = 0;
166:               pe = lz[0];
167:               pr = 0;
168:             } else {
169:               pr++;
170:               ps = pe;
171:               pe += lz[pr];
172:             }
173:           }
174:           /* compute the vector base on owning processor */
175:           xm   = me - ms;
176:           ym   = ne - ns;
177:           zm   = pe - ps;
178:           base = ms * ym * zm + ns * M * zm + ps * M * N;
179:           /* compute the local coordinates on owning processor */
180:           si = ii - ms;
181:           sj = jj - ns;
182:           sk = kk - ps;
183:           for (l = 0; l < dof; l++) {
184:             indices[idx] = l + dof * (base + si + xm * sj + xm * ym * sk);
185:             idx++;
186:           }
187:           i++;
188:         } while (i < upper->i - ox);
189:         j++;
190:       } while (j < upper->j - oy);
191:       k++;
192:     } while (k < upper->k - oz);
193:   }

195:   if (!offproc) {
196:     PetscCall(DMDAGetCorners(da, &ms, &ns, &ps, &mw, &nw, &pw));
197:     me = ms + mw;
198:     if (N > 1) ne = ns + nw;
199:     if (P > 1) pe = ps + pw;
200:     /* Account for DM offsets */
201:     ms = ms - ox;
202:     me = me - ox;
203:     ns = ns - oy;
204:     ne = ne - oy;
205:     ps = ps - oz;
206:     pe = pe - oz;

208:     /* compute the vector base on owning processor */
209:     xm   = me - ms;
210:     ym   = ne - ns;
211:     zm   = pe - ps;
212:     base = ms * ym * zm + ns * M * zm + ps * M * N;
213:     /*
214:        if no indices are to be returned, create an empty is,
215:        this prevents hanging in while loops
216:     */
217:     if (skip_i && skip_j && skip_k) goto createis;
218:     /*
219:        do..while loops to ensure the block gets entered once,
220:        regardless of control condition being met, necessary for
221:        cases when a subset of skip_i/j/k is true
222:     */
223:     if (skip_k) k = upper->k - oz;
224:     else k = lower->k - oz;
225:     do {
226:       if (skip_j) j = upper->j - oy;
227:       else j = lower->j - oy;
228:       do {
229:         if (skip_i) i = upper->i - ox;
230:         else i = lower->i - ox;
231:         do {
232:           if (k >= ps && k <= pe - 1) {
233:             if (j >= ns && j <= ne - 1) {
234:               if (i >= ms && i <= me - 1) {
235:                 /* compute the local coordinates on owning processor */
236:                 si = i - ms;
237:                 sj = j - ns;
238:                 sk = k - ps;
239:                 for (l = 0; l < dof; l++) {
240:                   indices[idx] = l + dof * (base + si + xm * sj + xm * ym * sk);
241:                   idx++;
242:                 }
243:               }
244:             }
245:           }
246:           i++;
247:         } while (i < upper->i - ox);
248:         j++;
249:       } while (j < upper->j - oy);
250:       k++;
251:     } while (k < upper->k - oz);

253:     PetscCall(PetscRealloc((size_t)(idx * sizeof(PetscInt)), (void *)&indices));
254:   }

256: createis:
257:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)da), idx, indices, PETSC_OWN_POINTER, is));
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm)
262: {
263:   DM           *da;
264:   PetscInt      dim, size, i, j, k, idx;
265:   DMDALocalInfo info;
266:   PetscInt      xsize, ysize, zsize;
267:   PetscInt      xo, yo, zo;
268:   PetscInt      xs, ys, zs;
269:   PetscInt      xm = 1, ym = 1, zm = 1;
270:   PetscInt      xol, yol, zol;
271:   PetscInt      m = 1, n = 1, p = 1;
272:   PetscInt      M, N, P;
273:   PetscInt      pm, mtmp;

275:   PetscFunctionBegin;
276:   PetscCall(DMDAGetLocalInfo(dm, &info));
277:   PetscCall(DMDAGetOverlap(dm, &xol, &yol, &zol));
278:   PetscCall(DMDAGetNumLocalSubDomains(dm, &size));
279:   PetscCall(PetscMalloc1(size, &da));

281:   if (nlocal) *nlocal = size;

283:   dim = info.dim;

285:   M = info.xm;
286:   N = info.ym;
287:   P = info.zm;

289:   if (dim == 1) {
290:     m = size;
291:   } else if (dim == 2) {
292:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N)));
293:     while (m > 0) {
294:       n = size / m;
295:       if (m * n * p == size) break;
296:       m--;
297:     }
298:   } else if (dim == 3) {
299:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N * N) * ((PetscReal)size) / ((PetscReal)P * M), (PetscReal)(1. / 3.)));
300:     if (!n) n = 1;
301:     while (n > 0) {
302:       pm = size / n;
303:       if (n * pm == size) break;
304:       n--;
305:     }
306:     if (!n) n = 1;
307:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
308:     if (!m) m = 1;
309:     while (m > 0) {
310:       p = size / (m * n);
311:       if (m * n * p == size) break;
312:       m--;
313:     }
314:     if (M > P && m < p) {
315:       mtmp = m;
316:       m    = p;
317:       p    = mtmp;
318:     }
319:   }

321:   zs  = info.zs;
322:   idx = 0;
323:   for (k = 0; k < p; k++) {
324:     ys = info.ys;
325:     for (j = 0; j < n; j++) {
326:       xs = info.xs;
327:       for (i = 0; i < m; i++) {
328:         if (dim == 1) {
329:           xm = M / m + ((M % m) > i);
330:         } else if (dim == 2) {
331:           xm = M / m + ((M % m) > i);
332:           ym = N / n + ((N % n) > j);
333:         } else if (dim == 3) {
334:           xm = M / m + ((M % m) > i);
335:           ym = N / n + ((N % n) > j);
336:           zm = P / p + ((P % p) > k);
337:         }

339:         xsize = xm;
340:         ysize = ym;
341:         zsize = zm;
342:         xo    = xs;
343:         yo    = ys;
344:         zo    = zs;

346:         PetscCall(DMDACreate(PETSC_COMM_SELF, &(da[idx])));
347:         PetscCall(DMSetOptionsPrefix(da[idx], "sub_"));
348:         PetscCall(DMSetDimension(da[idx], info.dim));
349:         PetscCall(DMDASetDof(da[idx], info.dof));

351:         PetscCall(DMDASetStencilType(da[idx], info.st));
352:         PetscCall(DMDASetStencilWidth(da[idx], info.sw));

354:         if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) {
355:           xsize += xol;
356:           xo -= xol;
357:         }
358:         if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) {
359:           ysize += yol;
360:           yo -= yol;
361:         }
362:         if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) {
363:           zsize += zol;
364:           zo -= zol;
365:         }

367:         if (info.bx == DM_BOUNDARY_PERIODIC || (xs + xm != info.mx)) xsize += xol;
368:         if (info.by == DM_BOUNDARY_PERIODIC || (ys + ym != info.my)) ysize += yol;
369:         if (info.bz == DM_BOUNDARY_PERIODIC || (zs + zm != info.mz)) zsize += zol;

371:         if (info.bx != DM_BOUNDARY_PERIODIC) {
372:           if (xo < 0) {
373:             xsize += xo;
374:             xo = 0;
375:           }
376:           if (xo + xsize > info.mx - 1) xsize -= xo + xsize - info.mx;
377:         }
378:         if (info.by != DM_BOUNDARY_PERIODIC) {
379:           if (yo < 0) {
380:             ysize += yo;
381:             yo = 0;
382:           }
383:           if (yo + ysize > info.my - 1) ysize -= yo + ysize - info.my;
384:         }
385:         if (info.bz != DM_BOUNDARY_PERIODIC) {
386:           if (zo < 0) {
387:             zsize += zo;
388:             zo = 0;
389:           }
390:           if (zo + zsize > info.mz - 1) zsize -= zo + zsize - info.mz;
391:         }

393:         PetscCall(DMDASetSizes(da[idx], xsize, ysize, zsize));
394:         PetscCall(DMDASetNumProcs(da[idx], 1, 1, 1));
395:         PetscCall(DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED));

397:         /* set up as a block instead */
398:         PetscCall(DMSetUp(da[idx]));

400:         /* nonoverlapping region */
401:         PetscCall(DMDASetNonOverlappingRegion(da[idx], xs, ys, zs, xm, ym, zm));

403:         /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */
404:         PetscCall(DMDASetOffset(da[idx], xo, yo, zo, info.mx, info.my, info.mz));
405:         xs += xm;
406:         idx++;
407:       }
408:       ys += ym;
409:     }
410:     zs += zm;
411:   }
412:   *sdm = da;
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: /*
417:    Fills the local vector problem on the subdomain from the global problem.

419:    Right now this assumes one subdomain per processor.

421: */
422: PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm, PetscInt nsubdms, DM *subdms, VecScatter **iscat, VecScatter **oscat, VecScatter **lscat)
423: {
424:   DMDALocalInfo info, subinfo;
425:   DM            subdm;
426:   MatStencil    upper, lower;
427:   IS            idis, isis, odis, osis, gdis;
428:   Vec           svec, dvec, slvec;
429:   PetscInt      xm, ym, zm, xs, ys, zs;
430:   PetscInt      i;
431:   PetscBool     patchis_offproc = PETSC_TRUE;

433:   PetscFunctionBegin;
434:   /* allocate the arrays of scatters */
435:   if (iscat) PetscCall(PetscMalloc1(nsubdms, iscat));
436:   if (oscat) PetscCall(PetscMalloc1(nsubdms, oscat));
437:   if (lscat) PetscCall(PetscMalloc1(nsubdms, lscat));

439:   PetscCall(DMDAGetLocalInfo(dm, &info));
440:   for (i = 0; i < nsubdms; i++) {
441:     subdm = subdms[i];
442:     PetscCall(DMDAGetLocalInfo(subdm, &subinfo));
443:     PetscCall(DMDAGetNonOverlappingRegion(subdm, &xs, &ys, &zs, &xm, &ym, &zm));

445:     /* create the global and subdomain index sets for the inner domain */
446:     lower.i = xs;
447:     lower.j = ys;
448:     lower.k = zs;
449:     upper.i = xs + xm;
450:     upper.j = ys + ym;
451:     upper.k = zs + zm;
452:     PetscCall(DMDACreatePatchIS(dm, &lower, &upper, &idis, patchis_offproc));
453:     PetscCall(DMDACreatePatchIS(subdm, &lower, &upper, &isis, patchis_offproc));

455:     /* create the global and subdomain index sets for the outer subdomain */
456:     lower.i = subinfo.xs;
457:     lower.j = subinfo.ys;
458:     lower.k = subinfo.zs;
459:     upper.i = subinfo.xs + subinfo.xm;
460:     upper.j = subinfo.ys + subinfo.ym;
461:     upper.k = subinfo.zs + subinfo.zm;
462:     PetscCall(DMDACreatePatchIS(dm, &lower, &upper, &odis, patchis_offproc));
463:     PetscCall(DMDACreatePatchIS(subdm, &lower, &upper, &osis, patchis_offproc));

465:     /* global and subdomain ISes for the local indices of the subdomain */
466:     /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */
467:     lower.i = subinfo.gxs;
468:     lower.j = subinfo.gys;
469:     lower.k = subinfo.gzs;
470:     upper.i = subinfo.gxs + subinfo.gxm;
471:     upper.j = subinfo.gys + subinfo.gym;
472:     upper.k = subinfo.gzs + subinfo.gzm;
473:     PetscCall(DMDACreatePatchIS(dm, &lower, &upper, &gdis, patchis_offproc));

475:     /* form the scatter */
476:     PetscCall(DMGetGlobalVector(dm, &dvec));
477:     PetscCall(DMGetGlobalVector(subdm, &svec));
478:     PetscCall(DMGetLocalVector(subdm, &slvec));

480:     if (iscat) PetscCall(VecScatterCreate(dvec, idis, svec, isis, &(*iscat)[i]));
481:     if (oscat) PetscCall(VecScatterCreate(dvec, odis, svec, osis, &(*oscat)[i]));
482:     if (lscat) PetscCall(VecScatterCreate(dvec, gdis, slvec, NULL, &(*lscat)[i]));

484:     PetscCall(DMRestoreGlobalVector(dm, &dvec));
485:     PetscCall(DMRestoreGlobalVector(subdm, &svec));
486:     PetscCall(DMRestoreLocalVector(subdm, &slvec));

488:     PetscCall(ISDestroy(&idis));
489:     PetscCall(ISDestroy(&isis));

491:     PetscCall(ISDestroy(&odis));
492:     PetscCall(ISDestroy(&osis));

494:     PetscCall(ISDestroy(&gdis));
495:   }
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: }

499: PetscErrorCode DMDASubDomainIS_Private(DM dm, PetscInt n, DM *subdm, IS **iis, IS **ois)
500: {
501:   PetscInt      i;
502:   DMDALocalInfo info, subinfo;
503:   MatStencil    lower, upper;
504:   PetscBool     patchis_offproc = PETSC_TRUE;

506:   PetscFunctionBegin;
507:   PetscCall(DMDAGetLocalInfo(dm, &info));
508:   if (iis) PetscCall(PetscMalloc1(n, iis));
509:   if (ois) PetscCall(PetscMalloc1(n, ois));

511:   for (i = 0; i < n; i++) {
512:     PetscCall(DMDAGetLocalInfo(subdm[i], &subinfo));
513:     if (iis) {
514:       /* create the inner IS */
515:       lower.i = info.xs;
516:       lower.j = info.ys;
517:       lower.k = info.zs;
518:       upper.i = info.xs + info.xm;
519:       upper.j = info.ys + info.ym;
520:       upper.k = info.zs + info.zm;
521:       PetscCall(DMDACreatePatchIS(dm, &lower, &upper, &(*iis)[i], patchis_offproc));
522:     }

524:     if (ois) {
525:       /* create the outer IS */
526:       lower.i = subinfo.xs;
527:       lower.j = subinfo.ys;
528:       lower.k = subinfo.zs;
529:       upper.i = subinfo.xs + subinfo.xm;
530:       upper.j = subinfo.ys + subinfo.ym;
531:       upper.k = subinfo.zs + subinfo.zm;
532:       PetscCall(DMDACreatePatchIS(dm, &lower, &upper, &(*ois)[i], patchis_offproc));
533:     }
534:   }
535:   PetscFunctionReturn(PETSC_SUCCESS);
536: }

538: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm, PetscInt *len, char ***names, IS **iis, IS **ois, DM **subdm)
539: {
540:   DM      *sdm;
541:   PetscInt n, i;

543:   PetscFunctionBegin;
544:   PetscCall(DMDASubDomainDA_Private(dm, &n, &sdm));
545:   if (names) {
546:     PetscCall(PetscMalloc1(n, names));
547:     for (i = 0; i < n; i++) (*names)[i] = NULL;
548:   }
549:   PetscCall(DMDASubDomainIS_Private(dm, n, sdm, iis, ois));
550:   if (subdm) *subdm = sdm;
551:   else {
552:     for (i = 0; i < n; i++) PetscCall(DMDestroy(&sdm[i]));
553:   }
554:   if (len) *len = n;
555:   PetscFunctionReturn(PETSC_SUCCESS);
556: }