1 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
2
3 /*@
4 DMDACreatePatchIS - Creates an index set corresponding to a logically rectangular patch of the `DMDA`.
5
6 Collective
7
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
13
14 Output Parameter:
15 . is - the `IS` corresponding to the patch
16
17 Level: developer
18
19 Notes:
20 This routine always returns an `IS` on the `DMDA` communicator.
21
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.
26
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.
30
31 .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMCreateDomainDecompositionScatters()`
32 @*/
DMDACreatePatchIS(DM da,MatStencil * lower,MatStencil * upper,IS * is,PetscBool offproc)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 */
52
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;
61
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;
84
85 PetscCall(PetscMalloc1(nindices * dof, &indices));
86 PetscCall(DMDAGetOffset(da, &ox, &oy, &oz, NULL, NULL, NULL));
87
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 }
96
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;
123
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 }
196
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;
209
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;
229
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);
256
257 PetscCall(PetscRealloc(idx * sizeof(PetscInt), &indices));
258 }
259
260 createis:
261 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)da), idx, indices, PETSC_OWN_POINTER, is));
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
DMDASubDomainDA_Private(DM dm,PetscInt * nlocal,DM ** sdm)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;
278
279 PetscFunctionBegin;
280 PetscCall(DMDAGetLocalInfo(dm, &info));
281 PetscCall(DMDAGetOverlap(dm, &xol, &yol, &zol));
282 PetscCall(DMDAGetNumLocalSubDomains(dm, &size));
283 PetscCall(PetscMalloc1(size, &da));
284
285 if (nlocal) *nlocal = size;
286
287 dim = info.dim;
288
289 M = info.xm;
290 N = info.ym;
291 P = info.zm;
292
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 }
324
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 }
342
343 xsize = xm;
344 ysize = ym;
345 zsize = zm;
346 xo = xs;
347 yo = ys;
348 zo = zs;
349
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));
354
355 PetscCall(DMDASetStencilType(da[idx], info.st));
356 PetscCall(DMDASetStencilWidth(da[idx], info.sw));
357
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 }
370
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;
374
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 }
396
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));
400
401 /* set up as a block instead */
402 PetscCall(DMSetUp(da[idx]));
403
404 /* nonoverlapping region */
405 PetscCall(DMDASetNonOverlappingRegion(da[idx], xs, ys, zs, xm, ym, zm));
406
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 }
419
420 /*
421 Fills the local vector problem on the subdomain from the global problem.
422
423 Right now this assumes one subdomain per processor.
424
425 */
DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM * subdms,VecScatter ** iscat,VecScatter ** oscat,VecScatter ** lscat)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;
436
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));
442
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));
448
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));
458
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));
468
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));
478
479 /* form the scatter */
480 PetscCall(DMGetGlobalVector(dm, &dvec));
481 PetscCall(DMGetGlobalVector(subdm, &svec));
482 PetscCall(DMGetLocalVector(subdm, &slvec));
483
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]));
487
488 PetscCall(DMRestoreGlobalVector(dm, &dvec));
489 PetscCall(DMRestoreGlobalVector(subdm, &svec));
490 PetscCall(DMRestoreLocalVector(subdm, &slvec));
491
492 PetscCall(ISDestroy(&idis));
493 PetscCall(ISDestroy(&isis));
494
495 PetscCall(ISDestroy(&odis));
496 PetscCall(ISDestroy(&osis));
497
498 PetscCall(ISDestroy(&gdis));
499 }
500 PetscFunctionReturn(PETSC_SUCCESS);
501 }
502
DMDASubDomainIS_Private(DM dm,PetscInt n,DM * subdm,IS ** iis,IS ** ois)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;
509
510 PetscFunctionBegin;
511 PetscCall(DMDAGetLocalInfo(dm, &info));
512 if (iis) PetscCall(PetscMalloc1(n, iis));
513 if (ois) PetscCall(PetscMalloc1(n, ois));
514
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 }
527
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 }
541
DMCreateDomainDecomposition_DA(DM dm,PetscInt * len,char *** names,IS ** iis,IS ** ois,DM ** subdm)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;
546
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 }
561