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 @*/ 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 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 */ 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 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 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