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