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, k = 0, 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: PetscCheck(PetscLikely(nindices < 0), PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Lower and Upper stencils are identical! Please check inputs.");
82: if (PetscUnlikely(skip_i && skip_j && skip_k)) nindices = 0;
83: else nindices = -nindices;
85: PetscCall(PetscMalloc1(nindices * dof, &indices));
86: PetscCall(DMDAGetOffset(da, &ox, &oy, &oz, NULL, NULL, NULL));
88: if (!valid_k) {
89: upper->k = 0;
90: lower->k = 0;
91: }
92: if (!valid_j) {
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: PetscInt j;
124: if (skip_j) j = upper->j - oy;
125: else j = lower->j - oy;
126: do {
127: if (skip_i) i = upper->i - ox;
128: else i = lower->i - ox;
129: do {
130: /* "actual" indices rather than ones outside of the domain */
131: ii = i;
132: jj = j;
133: kk = k;
134: if (ii < 0) ii = M + ii;
135: if (jj < 0) jj = N + jj;
136: if (kk < 0) kk = P + kk;
137: if (ii > M - 1) ii = ii - M;
138: if (jj > N - 1) jj = jj - N;
139: if (kk > P - 1) kk = kk - P;
140: /* gone out of processor range on x axis */
141: while (ii > me - 1 || ii < ms) {
142: if (mr == m - 1) {
143: ms = 0;
144: me = lx[0];
145: mr = 0;
146: } else {
147: mr++;
148: ms = me;
149: me += lx[mr];
150: }
151: }
152: /* gone out of processor range on y axis */
153: while (jj > ne - 1 || jj < ns) {
154: if (nr == n - 1) {
155: ns = 0;
156: ne = ly[0];
157: nr = 0;
158: } else {
159: nr++;
160: ns = ne;
161: ne += ly[nr];
162: }
163: }
164: /* gone out of processor range on z axis */
165: while (kk > pe - 1 || kk < ps) {
166: if (pr == p - 1) {
167: ps = 0;
168: pe = lz[0];
169: pr = 0;
170: } else {
171: pr++;
172: ps = pe;
173: pe += lz[pr];
174: }
175: }
176: /* compute the vector base on owning processor */
177: xm = me - ms;
178: ym = ne - ns;
179: zm = pe - ps;
180: base = ms * ym * zm + ns * M * zm + ps * M * N;
181: /* compute the local coordinates on owning processor */
182: si = ii - ms;
183: sj = jj - ns;
184: sk = kk - ps;
185: for (l = 0; l < dof; l++) {
186: indices[idx] = l + dof * (base + si + xm * sj + xm * ym * sk);
187: idx++;
188: }
189: i++;
190: } while (i < upper->i - ox);
191: j++;
192: } while (j < upper->j - oy);
193: k++;
194: } while (k < upper->k - oz);
195: }
197: if (!offproc) {
198: PetscCall(DMDAGetCorners(da, &ms, &ns, &ps, &mw, &nw, &pw));
199: me = ms + mw;
200: if (N > 1) ne = ns + nw;
201: if (P > 1) pe = ps + pw;
202: /* Account for DM offsets */
203: ms = ms - ox;
204: me = me - ox;
205: ns = ns - oy;
206: ne = ne - oy;
207: ps = ps - oz;
208: pe = pe - oz;
210: /* compute the vector base on owning processor */
211: xm = me - ms;
212: ym = ne - ns;
213: zm = pe - ps;
214: base = ms * ym * zm + ns * M * zm + ps * M * N;
215: /*
216: if no indices are to be returned, create an empty is,
217: this prevents hanging in while loops
218: */
219: if (skip_i && skip_j && skip_k) goto createis;
220: /*
221: do..while loops to ensure the block gets entered once,
222: regardless of control condition being met, necessary for
223: cases when a subset of skip_i/j/k is true
224: */
225: if (skip_k) k = upper->k - oz;
226: else k = lower->k - oz;
227: do {
228: PetscInt j;
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(idx * sizeof(PetscInt), &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: }