Actual source code: dadd.c

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

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

  6:   Collective

  8:   Input Parameters:
  9: + da      - the `DMDA`
 10: . lower   - a `MatStencil` with i, j and k entries corresponding to the lower corner of the patch
 11: . upper   - a `MatStencil` with i, j and k entries 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` communicator.

 22:   If `offproc` is set to `PETSC_TRUE`,
 23:   the routine returns an `IS` with all the indices requested regardless of whether these indices
 24:   are present on the requesting MPI process or not. Thus, it is upon the caller to ensure that
 25:   the indices returned in this mode are appropriate.

 27:   If `offproc` is set to `PETSC_FALSE`,
 28:   the `IS` only returns the subset of indices that are present on the requesting MPI process and there
 29:   is no duplication of indices between multiple MPI processes.

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

 53:   PetscFunctionBegin;
 54:   M   = dd->M;
 55:   N   = dd->N;
 56:   P   = dd->P;
 57:   m   = dd->m;
 58:   n   = dd->n;
 59:   p   = dd->p;
 60:   dof = dd->w;

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

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

 90:   if (!valid_k) {
 91:     k        = 0;
 92:     upper->k = 0;
 93:     lower->k = 0;
 94:   }
 95:   if (!valid_j) {
 96:     j        = 0;
 97:     upper->j = 0;
 98:     lower->j = 0;
 99:   }

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

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

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

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

260: createis:
261:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)da), idx, indices, PETSC_OWN_POINTER, is));
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

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

279:   PetscFunctionBegin;
280:   PetscCall(DMDAGetLocalInfo(dm, &info));
281:   PetscCall(DMDAGetOverlap(dm, &xol, &yol, &zol));
282:   PetscCall(DMDAGetNumLocalSubDomains(dm, &size));
283:   PetscCall(PetscMalloc1(size, &da));

285:   if (nlocal) *nlocal = size;

287:   dim = info.dim;

289:   M = info.xm;
290:   N = info.ym;
291:   P = info.zm;

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

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

343:         xsize = xm;
344:         ysize = ym;
345:         zsize = zm;
346:         xo    = xs;
347:         yo    = ys;
348:         zo    = zs;

350:         PetscCall(DMDACreate(PETSC_COMM_SELF, &(da[idx])));
351:         PetscCall(DMSetOptionsPrefix(da[idx], "sub_"));
352:         PetscCall(DMSetDimension(da[idx], info.dim));
353:         PetscCall(DMDASetDof(da[idx], info.dof));

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

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

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

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

397:         PetscCall(DMDASetSizes(da[idx], xsize, ysize, zsize));
398:         PetscCall(DMDASetNumProcs(da[idx], 1, 1, 1));
399:         PetscCall(DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED));

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

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

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

420: /*
421:    Fills the local vector problem on the subdomain from the global problem.

423:    Right now this assumes one subdomain per processor.

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

437:   PetscFunctionBegin;
438:   /* allocate the arrays of scatters */
439:   if (iscat) PetscCall(PetscMalloc1(nsubdms, iscat));
440:   if (oscat) PetscCall(PetscMalloc1(nsubdms, oscat));
441:   if (lscat) PetscCall(PetscMalloc1(nsubdms, lscat));

443:   PetscCall(DMDAGetLocalInfo(dm, &info));
444:   for (i = 0; i < nsubdms; i++) {
445:     subdm = subdms[i];
446:     PetscCall(DMDAGetLocalInfo(subdm, &subinfo));
447:     PetscCall(DMDAGetNonOverlappingRegion(subdm, &xs, &ys, &zs, &xm, &ym, &zm));

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

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

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

479:     /* form the scatter */
480:     PetscCall(DMGetGlobalVector(dm, &dvec));
481:     PetscCall(DMGetGlobalVector(subdm, &svec));
482:     PetscCall(DMGetLocalVector(subdm, &slvec));

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

488:     PetscCall(DMRestoreGlobalVector(dm, &dvec));
489:     PetscCall(DMRestoreGlobalVector(subdm, &svec));
490:     PetscCall(DMRestoreLocalVector(subdm, &slvec));

492:     PetscCall(ISDestroy(&idis));
493:     PetscCall(ISDestroy(&isis));

495:     PetscCall(ISDestroy(&odis));
496:     PetscCall(ISDestroy(&osis));

498:     PetscCall(ISDestroy(&gdis));
499:   }
500:   PetscFunctionReturn(PETSC_SUCCESS);
501: }

503: static PetscErrorCode DMDASubDomainIS_Private(DM dm, PetscInt n, DM *subdm, IS **iis, IS **ois)
504: {
505:   PetscInt      i;
506:   DMDALocalInfo info, subinfo;
507:   MatStencil    lower, upper;
508:   PetscBool     patchis_offproc = PETSC_TRUE;

510:   PetscFunctionBegin;
511:   PetscCall(DMDAGetLocalInfo(dm, &info));
512:   if (iis) PetscCall(PetscMalloc1(n, iis));
513:   if (ois) PetscCall(PetscMalloc1(n, ois));

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

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

542: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm, PetscInt *len, char ***names, IS **iis, IS **ois, DM **subdm)
543: {
544:   DM      *sdm = NULL;
545:   PetscInt n   = 0, i;

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