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 Parameters:
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: M = dd->M;
50: N = dd->N;
51: P = dd->P;
52: m = dd->m;
53: n = dd->n;
54: p = dd->p;
55: dof = dd->w;
57: nindices = -1;
58: if (PetscLikely(upper->i - lower->i)) {
59: nindices = nindices * (upper->i - lower->i);
60: skip_i = PETSC_FALSE;
61: }
62: if (N > 1) {
63: valid_j = PETSC_TRUE;
64: if (PetscLikely(upper->j - lower->j)) {
65: nindices = nindices * (upper->j - lower->j);
66: skip_j = PETSC_FALSE;
67: }
68: }
69: if (P > 1) {
70: valid_k = PETSC_TRUE;
71: if (PetscLikely(upper->k - lower->k)) {
72: nindices = nindices * (upper->k - lower->k);
73: skip_k = PETSC_FALSE;
74: }
75: }
76: if (PetscLikely(nindices < 0)) {
77: if (PetscUnlikely(skip_i && skip_j && skip_k)) {
78: nindices = 0;
79: } else nindices = nindices * (-1);
80: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Lower and Upper stencils are identical! Please check inputs.");
82: PetscMalloc1(nindices * dof, &indices);
83: DMDAGetOffset(da, &ox, &oy, &oz, NULL, NULL, NULL);
85: if (!valid_k) {
86: k = 0;
87: upper->k = 0;
88: lower->k = 0;
89: }
90: if (!valid_j) {
91: j = 0;
92: upper->j = 0;
93: lower->j = 0;
94: }
96: if (offproc) {
97: DMDAGetOwnershipRanges(da, &lx, &ly, &lz);
98: /* start at index 0 on processor 0 */
99: mr = 0;
100: nr = 0;
101: pr = 0;
102: ms = 0;
103: ns = 0;
104: ps = 0;
105: if (lx) me = lx[0];
106: if (ly) ne = ly[0];
107: if (lz) pe = lz[0];
108: /*
109: If no indices are to be returned, create an empty is,
110: this prevents hanging in while loops
111: */
112: if (skip_i && skip_j && skip_k) goto createis;
113: /*
114: do..while loops to ensure the block gets entered once,
115: regardless of control condition being met, necessary for
116: cases when a subset of skip_i/j/k is true
117: */
118: if (skip_k) k = upper->k - oz;
119: else k = lower->k - oz;
120: do {
121: if (skip_j) j = upper->j - oy;
122: else j = lower->j - oy;
123: do {
124: if (skip_i) i = upper->i - ox;
125: else i = lower->i - ox;
126: do {
127: /* "actual" indices rather than ones outside of the domain */
128: ii = i;
129: jj = j;
130: kk = k;
131: if (ii < 0) ii = M + ii;
132: if (jj < 0) jj = N + jj;
133: if (kk < 0) kk = P + kk;
134: if (ii > M - 1) ii = ii - M;
135: if (jj > N - 1) jj = jj - N;
136: if (kk > P - 1) kk = kk - P;
137: /* gone out of processor range on x axis */
138: while (ii > me - 1 || ii < ms) {
139: if (mr == m - 1) {
140: ms = 0;
141: me = lx[0];
142: mr = 0;
143: } else {
144: mr++;
145: ms = me;
146: me += lx[mr];
147: }
148: }
149: /* gone out of processor range on y axis */
150: while (jj > ne - 1 || jj < ns) {
151: if (nr == n - 1) {
152: ns = 0;
153: ne = ly[0];
154: nr = 0;
155: } else {
156: nr++;
157: ns = ne;
158: ne += ly[nr];
159: }
160: }
161: /* gone out of processor range on z axis */
162: while (kk > pe - 1 || kk < ps) {
163: if (pr == p - 1) {
164: ps = 0;
165: pe = lz[0];
166: pr = 0;
167: } else {
168: pr++;
169: ps = pe;
170: pe += lz[pr];
171: }
172: }
173: /* compute the vector base on owning processor */
174: xm = me - ms;
175: ym = ne - ns;
176: zm = pe - ps;
177: base = ms * ym * zm + ns * M * zm + ps * M * N;
178: /* compute the local coordinates on owning processor */
179: si = ii - ms;
180: sj = jj - ns;
181: sk = kk - ps;
182: for (l = 0; l < dof; l++) {
183: indices[idx] = l + dof * (base + si + xm * sj + xm * ym * sk);
184: idx++;
185: }
186: i++;
187: } while (i < upper->i - ox);
188: j++;
189: } while (j < upper->j - oy);
190: k++;
191: } while (k < upper->k - oz);
192: }
194: if (!offproc) {
195: DMDAGetCorners(da, &ms, &ns, &ps, &mw, &nw, &pw);
196: me = ms + mw;
197: if (N > 1) ne = ns + nw;
198: if (P > 1) pe = ps + pw;
199: /* Account for DM offsets */
200: ms = ms - ox;
201: me = me - ox;
202: ns = ns - oy;
203: ne = ne - oy;
204: ps = ps - oz;
205: pe = pe - oz;
207: /* compute the vector base on owning processor */
208: xm = me - ms;
209: ym = ne - ns;
210: zm = pe - ps;
211: base = ms * ym * zm + ns * M * zm + ps * M * N;
212: /*
213: if no indices are to be returned, create an empty is,
214: this prevents hanging in while loops
215: */
216: if (skip_i && skip_j && skip_k) goto createis;
217: /*
218: do..while loops to ensure the block gets entered once,
219: regardless of control condition being met, necessary for
220: cases when a subset of skip_i/j/k is true
221: */
222: if (skip_k) k = upper->k - oz;
223: else k = lower->k - oz;
224: do {
225: if (skip_j) j = upper->j - oy;
226: else j = lower->j - oy;
227: do {
228: if (skip_i) i = upper->i - ox;
229: else i = lower->i - ox;
230: do {
231: if (k >= ps && k <= pe - 1) {
232: if (j >= ns && j <= ne - 1) {
233: if (i >= ms && i <= me - 1) {
234: /* compute the local coordinates on owning processor */
235: si = i - ms;
236: sj = j - ns;
237: sk = k - ps;
238: for (l = 0; l < dof; l++) {
239: indices[idx] = l + dof * (base + si + xm * sj + xm * ym * sk);
240: idx++;
241: }
242: }
243: }
244: }
245: i++;
246: } while (i < upper->i - ox);
247: j++;
248: } while (j < upper->j - oy);
249: k++;
250: } while (k < upper->k - oz);
252: PetscRealloc((size_t)(idx * sizeof(PetscInt)), (void *)&indices);
253: }
255: createis:
256: ISCreateGeneral(PetscObjectComm((PetscObject)da), idx, indices, PETSC_OWN_POINTER, is);
257: return 0;
258: }
260: PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm)
261: {
262: DM *da;
263: PetscInt dim, size, i, j, k, idx;
264: DMDALocalInfo info;
265: PetscInt xsize, ysize, zsize;
266: PetscInt xo, yo, zo;
267: PetscInt xs, ys, zs;
268: PetscInt xm = 1, ym = 1, zm = 1;
269: PetscInt xol, yol, zol;
270: PetscInt m = 1, n = 1, p = 1;
271: PetscInt M, N, P;
272: PetscInt pm, mtmp;
274: DMDAGetLocalInfo(dm, &info);
275: DMDAGetOverlap(dm, &xol, &yol, &zol);
276: DMDAGetNumLocalSubDomains(dm, &size);
277: PetscMalloc1(size, &da);
279: if (nlocal) *nlocal = size;
281: dim = info.dim;
283: M = info.xm;
284: N = info.ym;
285: P = info.zm;
287: if (dim == 1) {
288: m = size;
289: } else if (dim == 2) {
290: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N)));
291: while (m > 0) {
292: n = size / m;
293: if (m * n * p == size) break;
294: m--;
295: }
296: } else if (dim == 3) {
297: n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N * N) * ((PetscReal)size) / ((PetscReal)P * M), (PetscReal)(1. / 3.)));
298: if (!n) n = 1;
299: while (n > 0) {
300: pm = size / n;
301: if (n * pm == size) break;
302: n--;
303: }
304: if (!n) n = 1;
305: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)P * n)));
306: if (!m) m = 1;
307: while (m > 0) {
308: p = size / (m * n);
309: if (m * n * p == size) break;
310: m--;
311: }
312: if (M > P && m < p) {
313: mtmp = m;
314: m = p;
315: p = mtmp;
316: }
317: }
319: zs = info.zs;
320: idx = 0;
321: for (k = 0; k < p; k++) {
322: ys = info.ys;
323: for (j = 0; j < n; j++) {
324: xs = info.xs;
325: for (i = 0; i < m; i++) {
326: if (dim == 1) {
327: xm = M / m + ((M % m) > i);
328: } else if (dim == 2) {
329: xm = M / m + ((M % m) > i);
330: ym = N / n + ((N % n) > j);
331: } else if (dim == 3) {
332: xm = M / m + ((M % m) > i);
333: ym = N / n + ((N % n) > j);
334: zm = P / p + ((P % p) > k);
335: }
337: xsize = xm;
338: ysize = ym;
339: zsize = zm;
340: xo = xs;
341: yo = ys;
342: zo = zs;
344: DMDACreate(PETSC_COMM_SELF, &(da[idx]));
345: DMSetOptionsPrefix(da[idx], "sub_");
346: DMSetDimension(da[idx], info.dim);
347: DMDASetDof(da[idx], info.dof);
349: DMDASetStencilType(da[idx], info.st);
350: DMDASetStencilWidth(da[idx], info.sw);
352: if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) {
353: xsize += xol;
354: xo -= xol;
355: }
356: if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) {
357: ysize += yol;
358: yo -= yol;
359: }
360: if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) {
361: zsize += zol;
362: zo -= zol;
363: }
365: if (info.bx == DM_BOUNDARY_PERIODIC || (xs + xm != info.mx)) xsize += xol;
366: if (info.by == DM_BOUNDARY_PERIODIC || (ys + ym != info.my)) ysize += yol;
367: if (info.bz == DM_BOUNDARY_PERIODIC || (zs + zm != info.mz)) zsize += zol;
369: if (info.bx != DM_BOUNDARY_PERIODIC) {
370: if (xo < 0) {
371: xsize += xo;
372: xo = 0;
373: }
374: if (xo + xsize > info.mx - 1) xsize -= xo + xsize - info.mx;
375: }
376: if (info.by != DM_BOUNDARY_PERIODIC) {
377: if (yo < 0) {
378: ysize += yo;
379: yo = 0;
380: }
381: if (yo + ysize > info.my - 1) ysize -= yo + ysize - info.my;
382: }
383: if (info.bz != DM_BOUNDARY_PERIODIC) {
384: if (zo < 0) {
385: zsize += zo;
386: zo = 0;
387: }
388: if (zo + zsize > info.mz - 1) zsize -= zo + zsize - info.mz;
389: }
391: DMDASetSizes(da[idx], xsize, ysize, zsize);
392: DMDASetNumProcs(da[idx], 1, 1, 1);
393: DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);
395: /* set up as a block instead */
396: DMSetUp(da[idx]);
398: /* nonoverlapping region */
399: DMDASetNonOverlappingRegion(da[idx], xs, ys, zs, xm, ym, zm);
401: /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */
402: DMDASetOffset(da[idx], xo, yo, zo, info.mx, info.my, info.mz);
403: xs += xm;
404: idx++;
405: }
406: ys += ym;
407: }
408: zs += zm;
409: }
410: *sdm = da;
411: return 0;
412: }
414: /*
415: Fills the local vector problem on the subdomain from the global problem.
417: Right now this assumes one subdomain per processor.
419: */
420: PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm, PetscInt nsubdms, DM *subdms, VecScatter **iscat, VecScatter **oscat, VecScatter **lscat)
421: {
422: DMDALocalInfo info, subinfo;
423: DM subdm;
424: MatStencil upper, lower;
425: IS idis, isis, odis, osis, gdis;
426: Vec svec, dvec, slvec;
427: PetscInt xm, ym, zm, xs, ys, zs;
428: PetscInt i;
429: PetscBool patchis_offproc = PETSC_TRUE;
431: /* allocate the arrays of scatters */
432: if (iscat) PetscMalloc1(nsubdms, iscat);
433: if (oscat) PetscMalloc1(nsubdms, oscat);
434: if (lscat) PetscMalloc1(nsubdms, lscat);
436: DMDAGetLocalInfo(dm, &info);
437: for (i = 0; i < nsubdms; i++) {
438: subdm = subdms[i];
439: DMDAGetLocalInfo(subdm, &subinfo);
440: DMDAGetNonOverlappingRegion(subdm, &xs, &ys, &zs, &xm, &ym, &zm);
442: /* create the global and subdomain index sets for the inner domain */
443: lower.i = xs;
444: lower.j = ys;
445: lower.k = zs;
446: upper.i = xs + xm;
447: upper.j = ys + ym;
448: upper.k = zs + zm;
449: DMDACreatePatchIS(dm, &lower, &upper, &idis, patchis_offproc);
450: DMDACreatePatchIS(subdm, &lower, &upper, &isis, patchis_offproc);
452: /* create the global and subdomain index sets for the outer subdomain */
453: lower.i = subinfo.xs;
454: lower.j = subinfo.ys;
455: lower.k = subinfo.zs;
456: upper.i = subinfo.xs + subinfo.xm;
457: upper.j = subinfo.ys + subinfo.ym;
458: upper.k = subinfo.zs + subinfo.zm;
459: DMDACreatePatchIS(dm, &lower, &upper, &odis, patchis_offproc);
460: DMDACreatePatchIS(subdm, &lower, &upper, &osis, patchis_offproc);
462: /* global and subdomain ISes for the local indices of the subdomain */
463: /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */
464: lower.i = subinfo.gxs;
465: lower.j = subinfo.gys;
466: lower.k = subinfo.gzs;
467: upper.i = subinfo.gxs + subinfo.gxm;
468: upper.j = subinfo.gys + subinfo.gym;
469: upper.k = subinfo.gzs + subinfo.gzm;
470: DMDACreatePatchIS(dm, &lower, &upper, &gdis, patchis_offproc);
472: /* form the scatter */
473: DMGetGlobalVector(dm, &dvec);
474: DMGetGlobalVector(subdm, &svec);
475: DMGetLocalVector(subdm, &slvec);
477: if (iscat) VecScatterCreate(dvec, idis, svec, isis, &(*iscat)[i]);
478: if (oscat) VecScatterCreate(dvec, odis, svec, osis, &(*oscat)[i]);
479: if (lscat) VecScatterCreate(dvec, gdis, slvec, NULL, &(*lscat)[i]);
481: DMRestoreGlobalVector(dm, &dvec);
482: DMRestoreGlobalVector(subdm, &svec);
483: DMRestoreLocalVector(subdm, &slvec);
485: ISDestroy(&idis);
486: ISDestroy(&isis);
488: ISDestroy(&odis);
489: ISDestroy(&osis);
491: ISDestroy(&gdis);
492: }
493: return 0;
494: }
496: PetscErrorCode DMDASubDomainIS_Private(DM dm, PetscInt n, DM *subdm, IS **iis, IS **ois)
497: {
498: PetscInt i;
499: DMDALocalInfo info, subinfo;
500: MatStencil lower, upper;
501: PetscBool patchis_offproc = PETSC_TRUE;
503: DMDAGetLocalInfo(dm, &info);
504: if (iis) PetscMalloc1(n, iis);
505: if (ois) PetscMalloc1(n, ois);
507: for (i = 0; i < n; i++) {
508: DMDAGetLocalInfo(subdm[i], &subinfo);
509: if (iis) {
510: /* create the inner IS */
511: lower.i = info.xs;
512: lower.j = info.ys;
513: lower.k = info.zs;
514: upper.i = info.xs + info.xm;
515: upper.j = info.ys + info.ym;
516: upper.k = info.zs + info.zm;
517: DMDACreatePatchIS(dm, &lower, &upper, &(*iis)[i], patchis_offproc);
518: }
520: if (ois) {
521: /* create the outer IS */
522: lower.i = subinfo.xs;
523: lower.j = subinfo.ys;
524: lower.k = subinfo.zs;
525: upper.i = subinfo.xs + subinfo.xm;
526: upper.j = subinfo.ys + subinfo.ym;
527: upper.k = subinfo.zs + subinfo.zm;
528: DMDACreatePatchIS(dm, &lower, &upper, &(*ois)[i], patchis_offproc);
529: }
530: }
531: return 0;
532: }
534: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm, PetscInt *len, char ***names, IS **iis, IS **ois, DM **subdm)
535: {
536: DM *sdm;
537: PetscInt n, i;
539: DMDASubDomainDA_Private(dm, &n, &sdm);
540: if (names) {
541: PetscMalloc1(n, names);
542: for (i = 0; i < n; i++) (*names)[i] = NULL;
543: }
544: DMDASubDomainIS_Private(dm, n, sdm, iis, ois);
545: if (subdm) *subdm = sdm;
546: else {
547: for (i = 0; i < n; i++) DMDestroy(&sdm[i]);
548: }
549: if (len) *len = n;
550: return 0;
551: }