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, 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 */ 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 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."); 86 87 PetscCall(PetscMalloc1(nindices * dof, &indices)); 88 PetscCall(DMDAGetOffset(da, &ox, &oy, &oz, NULL, NULL, NULL)); 89 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 } 100 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 } 198 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; 211 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); 256 257 PetscCall(PetscRealloc((size_t)(idx * sizeof(PetscInt)), (void *)&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