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