1 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 2 3 /*@ 4 DMDACreatePatchIS - Creates an index set corresponding to a patch of the DA. 5 6 Collective 7 8 Input Parameters: 9 + da - the DMDA 10 . lower - a matstencil with i, j and k corresponding to the lower corner of the patch 11 . upper - a matstencil with i, j and k corresponding to the upper corner of the patch 12 - offproc - indicate whether the returned IS will contain off process indices 13 14 Output Parameters: 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's comm, if offproc is set to PETSC_TRUE, 21 the routine returns an IS with all the indices requested regardless of whether these indices 22 are present on the requesting rank or not. Thus, it is upon the caller to ensure that 23 the indices returned in this mode are appropriate. If offproc is set to PETSC_FALSE, 24 the IS only returns the subset of indices that are present on the requesting rank and there 25 is no duplication of indices. 26 27 .seealso: DMDACreateDomainDecomposition(), DMDACreateDomainDecompositionScatters() 28 @*/ 29 PetscErrorCode DMDACreatePatchIS(DM da,MatStencil *lower,MatStencil *upper,IS *is, PetscBool offproc) 30 { 31 PetscInt ms=0,ns=0,ps=0; 32 PetscInt mw=0,nw=0,pw=0; 33 PetscInt me=1,ne=1,pe=1; 34 PetscInt mr=0,nr=0,pr=0; 35 PetscInt ii,jj,kk; 36 PetscInt si,sj,sk; 37 PetscInt i,j,k,l,idx=0; 38 PetscInt base; 39 PetscInt xm=1,ym=1,zm=1; 40 PetscInt ox,oy,oz; 41 PetscInt m,n,p,M,N,P,dof; 42 const PetscInt *lx,*ly,*lz; 43 PetscInt nindices; 44 PetscInt *indices; 45 DM_DA *dd = (DM_DA*)da->data; 46 PetscBool skip_i=PETSC_TRUE, skip_j=PETSC_TRUE, skip_k=PETSC_TRUE; 47 PetscBool valid_j=PETSC_FALSE, valid_k=PETSC_FALSE; /* DMDA has at least 1 dimension */ 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 M = dd->M; N = dd->N; P = dd->P; 52 m = dd->m; n = dd->n; p = dd->p; 53 dof = dd->w; 54 55 nindices = -1; 56 if (PetscLikely(upper->i - lower->i)) { 57 nindices = nindices*(upper->i - lower->i); 58 skip_i=PETSC_FALSE; 59 } 60 if (N>1) { 61 valid_j = PETSC_TRUE; 62 if (PetscLikely(upper->j - lower->j)) { 63 nindices = nindices*(upper->j - lower->j); 64 skip_j=PETSC_FALSE; 65 } 66 } 67 if (P>1) { 68 valid_k = PETSC_TRUE; 69 if (PetscLikely(upper->k - lower->k)) { 70 nindices = nindices*(upper->k - lower->k); 71 skip_k=PETSC_FALSE; 72 } 73 } 74 if (PetscLikely(nindices<0)) { 75 if (PetscUnlikely(skip_i && skip_j && skip_k)) { 76 nindices = 0; 77 } else nindices = nindices*(-1); 78 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Lower and Upper stencils are identical! Please check inputs."); 79 80 ierr = PetscMalloc1(nindices*dof,&indices);CHKERRQ(ierr); 81 ierr = DMDAGetOffset(da,&ox,&oy,&oz,NULL,NULL,NULL);CHKERRQ(ierr); 82 83 if (!valid_k) { 84 k = 0; 85 upper->k=0; 86 lower->k=0; 87 } 88 if (!valid_j) { 89 j = 0; 90 upper->j=0; 91 lower->j=0; 92 } 93 94 if (offproc) { 95 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,&lz);CHKERRQ(ierr); 96 /* start at index 0 on processor 0 */ 97 mr = 0; 98 nr = 0; 99 pr = 0; 100 ms = 0; 101 ns = 0; 102 ps = 0; 103 if (lx) me = lx[0]; 104 if (ly) ne = ly[0]; 105 if (lz) pe = lz[0]; 106 /* 107 If no indices are to be returned, create an empty is, 108 this prevents hanging in while loops 109 */ 110 if (skip_i && skip_j && skip_k) goto createis; 111 /* 112 do..while loops to ensure the block gets entered once, 113 regardless of control condition being met, necessary for 114 cases when a subset of skip_i/j/k is true 115 */ 116 if (skip_k) k = upper->k-oz; else k = lower->k-oz; 117 do { 118 if (skip_j) j = upper->j-oy; else j = lower->j-oy; 119 do { 120 if (skip_i) i = upper->i-ox; else i = lower->i-ox; 121 do { 122 /* "actual" indices rather than ones outside of the domain */ 123 ii = i; 124 jj = j; 125 kk = k; 126 if (ii < 0) ii = M + ii; 127 if (jj < 0) jj = N + jj; 128 if (kk < 0) kk = P + kk; 129 if (ii > M-1) ii = ii - M; 130 if (jj > N-1) jj = jj - N; 131 if (kk > P-1) kk = kk - P; 132 /* gone out of processor range on x axis */ 133 while (ii > me-1 || ii < ms) { 134 if (mr == m-1) { 135 ms = 0; 136 me = lx[0]; 137 mr = 0; 138 } else { 139 mr++; 140 ms = me; 141 me += lx[mr]; 142 } 143 } 144 /* gone out of processor range on y axis */ 145 while (jj > ne-1 || jj < ns) { 146 if (nr == n-1) { 147 ns = 0; 148 ne = ly[0]; 149 nr = 0; 150 } else { 151 nr++; 152 ns = ne; 153 ne += ly[nr]; 154 } 155 } 156 /* gone out of processor range on z axis */ 157 while (kk > pe-1 || kk < ps) { 158 if (pr == p-1) { 159 ps = 0; 160 pe = lz[0]; 161 pr = 0; 162 } else { 163 pr++; 164 ps = pe; 165 pe += lz[pr]; 166 } 167 } 168 /* compute the vector base on owning processor */ 169 xm = me - ms; 170 ym = ne - ns; 171 zm = pe - ps; 172 base = ms*ym*zm + ns*M*zm + ps*M*N; 173 /* compute the local coordinates on owning processor */ 174 si = ii - ms; 175 sj = jj - ns; 176 sk = kk - ps; 177 for (l=0;l<dof;l++) { 178 indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk); 179 idx++; 180 } 181 i++; 182 } while (i<upper->i-ox); 183 j++; 184 } while (j<upper->j-oy); 185 k++; 186 } while (k<upper->k-oz); 187 } 188 189 if (!offproc) { 190 ierr = DMDAGetCorners(da, &ms, &ns, &ps, &mw, &nw, &pw);CHKERRQ(ierr); 191 me = ms + mw; 192 if (N>1) ne = ns + nw; 193 if (P>1) pe = ps + pw; 194 /* Account for DM offsets */ 195 ms = ms - ox; me = me - ox; 196 ns = ns - oy; ne = ne - oy; 197 ps = ps - oz; pe = pe - oz; 198 199 /* compute the vector base on owning processor */ 200 xm = me - ms; 201 ym = ne - ns; 202 zm = pe - ps; 203 base = ms*ym*zm + ns*M*zm + ps*M*N; 204 /* 205 if no indices are to be returned, create an empty is, 206 this prevents hanging in while loops 207 */ 208 if (skip_i && skip_j && skip_k) goto createis; 209 /* 210 do..while loops to ensure the block gets entered once, 211 regardless of control condition being met, necessary for 212 cases when a subset of skip_i/j/k is true 213 */ 214 if (skip_k) k = upper->k-oz; else k = lower->k-oz; 215 do { 216 if (skip_j) j = upper->j-oy; else j = lower->j-oy; 217 do { 218 if (skip_i) i = upper->i-ox; else i = lower->i-ox; 219 do { 220 if (k>=ps && k<=pe-1) { 221 if (j>=ns && j<=ne-1) { 222 if (i>=ms && i<=me-1) { 223 /* compute the local coordinates on owning processor */ 224 si = i - ms; 225 sj = j - ns; 226 sk = k - ps; 227 for (l=0; l<dof; l++) { 228 indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk); 229 idx++; 230 } 231 } 232 } 233 } 234 i++; 235 } while (i<upper->i-ox); 236 j++; 237 } while (j<upper->j-oy); 238 k++; 239 } while (k<upper->k-oz); 240 241 ierr = PetscRealloc((size_t)(idx*sizeof(PetscInt)), (void*)&indices);CHKERRQ(ierr); 242 } 243 244 createis: 245 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)da),idx,indices,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249 PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm) 250 { 251 DM *da; 252 PetscInt dim,size,i,j,k,idx; 253 PetscErrorCode ierr; 254 DMDALocalInfo info; 255 PetscInt xsize,ysize,zsize; 256 PetscInt xo,yo,zo; 257 PetscInt xs,ys,zs; 258 PetscInt xm=1,ym=1,zm=1; 259 PetscInt xol,yol,zol; 260 PetscInt m=1,n=1,p=1; 261 PetscInt M,N,P; 262 PetscInt pm,mtmp; 263 264 PetscFunctionBegin; 265 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 266 ierr = DMDAGetOverlap(dm,&xol,&yol,&zol);CHKERRQ(ierr); 267 ierr = DMDAGetNumLocalSubDomains(dm,&size);CHKERRQ(ierr); 268 ierr = PetscMalloc1(size,&da);CHKERRQ(ierr); 269 270 if (nlocal) *nlocal = size; 271 272 dim = info.dim; 273 274 M = info.xm; 275 N = info.ym; 276 P = info.zm; 277 278 if (dim == 1) { 279 m = size; 280 } else if (dim == 2) { 281 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N))); 282 while (m > 0) { 283 n = size/m; 284 if (m*n*p == size) break; 285 m--; 286 } 287 } else if (dim == 3) { 288 n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.))); if (!n) n = 1; 289 while (n > 0) { 290 pm = size/n; 291 if (n*pm == size) break; 292 n--; 293 } 294 if (!n) n = 1; 295 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n))); 296 if (!m) m = 1; 297 while (m > 0) { 298 p = size/(m*n); 299 if (m*n*p == size) break; 300 m--; 301 } 302 if (M > P && m < p) {mtmp = m; m = p; p = mtmp;} 303 } 304 305 zs = info.zs; 306 idx = 0; 307 for (k = 0; k < p; k++) { 308 ys = info.ys; 309 for (j = 0; j < n; j++) { 310 xs = info.xs; 311 for (i = 0; i < m; i++) { 312 if (dim == 1) { 313 xm = M/m + ((M % m) > i); 314 } else if (dim == 2) { 315 xm = M/m + ((M % m) > i); 316 ym = N/n + ((N % n) > j); 317 } else if (dim == 3) { 318 xm = M/m + ((M % m) > i); 319 ym = N/n + ((N % n) > j); 320 zm = P/p + ((P % p) > k); 321 } 322 323 xsize = xm; 324 ysize = ym; 325 zsize = zm; 326 xo = xs; 327 yo = ys; 328 zo = zs; 329 330 ierr = DMDACreate(PETSC_COMM_SELF,&(da[idx]));CHKERRQ(ierr); 331 ierr = DMSetOptionsPrefix(da[idx],"sub_");CHKERRQ(ierr); 332 ierr = DMSetDimension(da[idx], info.dim);CHKERRQ(ierr); 333 ierr = DMDASetDof(da[idx], info.dof);CHKERRQ(ierr); 334 335 ierr = DMDASetStencilType(da[idx],info.st);CHKERRQ(ierr); 336 ierr = DMDASetStencilWidth(da[idx],info.sw);CHKERRQ(ierr); 337 338 if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) { 339 xsize += xol; 340 xo -= xol; 341 } 342 if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) { 343 ysize += yol; 344 yo -= yol; 345 } 346 if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) { 347 zsize += zol; 348 zo -= zol; 349 } 350 351 if (info.bx == DM_BOUNDARY_PERIODIC || (xs+xm != info.mx)) xsize += xol; 352 if (info.by == DM_BOUNDARY_PERIODIC || (ys+ym != info.my)) ysize += yol; 353 if (info.bz == DM_BOUNDARY_PERIODIC || (zs+zm != info.mz)) zsize += zol; 354 355 if (info.bx != DM_BOUNDARY_PERIODIC) { 356 if (xo < 0) { 357 xsize += xo; 358 xo = 0; 359 } 360 if (xo+xsize > info.mx-1) { 361 xsize -= xo+xsize - info.mx; 362 } 363 } 364 if (info.by != DM_BOUNDARY_PERIODIC) { 365 if (yo < 0) { 366 ysize += yo; 367 yo = 0; 368 } 369 if (yo+ysize > info.my-1) { 370 ysize -= yo+ysize - info.my; 371 } 372 } 373 if (info.bz != DM_BOUNDARY_PERIODIC) { 374 if (zo < 0) { 375 zsize += zo; 376 zo = 0; 377 } 378 if (zo+zsize > info.mz-1) { 379 zsize -= zo+zsize - info.mz; 380 } 381 } 382 383 ierr = DMDASetSizes(da[idx], xsize, ysize, zsize);CHKERRQ(ierr); 384 ierr = DMDASetNumProcs(da[idx], 1, 1, 1);CHKERRQ(ierr); 385 ierr = DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);CHKERRQ(ierr); 386 387 /* set up as a block instead */ 388 ierr = DMSetUp(da[idx]);CHKERRQ(ierr); 389 390 /* nonoverlapping region */ 391 ierr = DMDASetNonOverlappingRegion(da[idx],xs,ys,zs,xm,ym,zm);CHKERRQ(ierr); 392 393 /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */ 394 ierr = DMDASetOffset(da[idx],xo,yo,zo,info.mx,info.my,info.mz);CHKERRQ(ierr); 395 xs += xm; 396 idx++; 397 } 398 ys += ym; 399 } 400 zs += zm; 401 } 402 *sdm = da; 403 PetscFunctionReturn(0); 404 } 405 406 /* 407 Fills the local vector problem on the subdomain from the global problem. 408 409 Right now this assumes one subdomain per processor. 410 411 */ 412 PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat) 413 { 414 PetscErrorCode ierr; 415 DMDALocalInfo info,subinfo; 416 DM subdm; 417 MatStencil upper,lower; 418 IS idis,isis,odis,osis,gdis; 419 Vec svec,dvec,slvec; 420 PetscInt xm,ym,zm,xs,ys,zs; 421 PetscInt i; 422 PetscBool patchis_offproc = PETSC_TRUE; 423 424 PetscFunctionBegin; 425 /* allocate the arrays of scatters */ 426 if (iscat) {ierr = PetscMalloc1(nsubdms,iscat);CHKERRQ(ierr);} 427 if (oscat) {ierr = PetscMalloc1(nsubdms,oscat);CHKERRQ(ierr);} 428 if (lscat) {ierr = PetscMalloc1(nsubdms,lscat);CHKERRQ(ierr);} 429 430 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 431 for (i = 0; i < nsubdms; i++) { 432 subdm = subdms[i]; 433 ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr); 434 ierr = DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 435 436 /* create the global and subdomain index sets for the inner domain */ 437 lower.i = xs; 438 lower.j = ys; 439 lower.k = zs; 440 upper.i = xs+xm; 441 upper.j = ys+ym; 442 upper.k = zs+zm; 443 ierr = DMDACreatePatchIS(dm,&lower,&upper,&idis,patchis_offproc);CHKERRQ(ierr); 444 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&isis,patchis_offproc);CHKERRQ(ierr); 445 446 /* create the global and subdomain index sets for the outer subdomain */ 447 lower.i = subinfo.xs; 448 lower.j = subinfo.ys; 449 lower.k = subinfo.zs; 450 upper.i = subinfo.xs+subinfo.xm; 451 upper.j = subinfo.ys+subinfo.ym; 452 upper.k = subinfo.zs+subinfo.zm; 453 ierr = DMDACreatePatchIS(dm,&lower,&upper,&odis,patchis_offproc);CHKERRQ(ierr); 454 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&osis,patchis_offproc);CHKERRQ(ierr); 455 456 /* global and subdomain ISes for the local indices of the subdomain */ 457 /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */ 458 lower.i = subinfo.gxs; 459 lower.j = subinfo.gys; 460 lower.k = subinfo.gzs; 461 upper.i = subinfo.gxs+subinfo.gxm; 462 upper.j = subinfo.gys+subinfo.gym; 463 upper.k = subinfo.gzs+subinfo.gzm; 464 ierr = DMDACreatePatchIS(dm,&lower,&upper,&gdis,patchis_offproc);CHKERRQ(ierr); 465 466 /* form the scatter */ 467 ierr = DMGetGlobalVector(dm,&dvec);CHKERRQ(ierr); 468 ierr = DMGetGlobalVector(subdm,&svec);CHKERRQ(ierr); 469 ierr = DMGetLocalVector(subdm,&slvec);CHKERRQ(ierr); 470 471 if (iscat) {ierr = VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[i]);CHKERRQ(ierr);} 472 if (oscat) {ierr = VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[i]);CHKERRQ(ierr);} 473 if (lscat) {ierr = VecScatterCreate(dvec,gdis,slvec,NULL,&(*lscat)[i]);CHKERRQ(ierr);} 474 475 ierr = DMRestoreGlobalVector(dm,&dvec);CHKERRQ(ierr); 476 ierr = DMRestoreGlobalVector(subdm,&svec);CHKERRQ(ierr); 477 ierr = DMRestoreLocalVector(subdm,&slvec);CHKERRQ(ierr); 478 479 ierr = ISDestroy(&idis);CHKERRQ(ierr); 480 ierr = ISDestroy(&isis);CHKERRQ(ierr); 481 482 ierr = ISDestroy(&odis);CHKERRQ(ierr); 483 ierr = ISDestroy(&osis);CHKERRQ(ierr); 484 485 ierr = ISDestroy(&gdis);CHKERRQ(ierr); 486 } 487 PetscFunctionReturn(0); 488 } 489 490 PetscErrorCode DMDASubDomainIS_Private(DM dm,PetscInt n,DM *subdm,IS **iis,IS **ois) 491 { 492 PetscErrorCode ierr; 493 PetscInt i; 494 DMDALocalInfo info,subinfo; 495 MatStencil lower,upper; 496 PetscBool patchis_offproc = PETSC_TRUE; 497 498 PetscFunctionBegin; 499 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 500 if (iis) {ierr = PetscMalloc1(n,iis);CHKERRQ(ierr);} 501 if (ois) {ierr = PetscMalloc1(n,ois);CHKERRQ(ierr);} 502 503 for (i = 0;i < n; i++) { 504 ierr = DMDAGetLocalInfo(subdm[i],&subinfo);CHKERRQ(ierr); 505 if (iis) { 506 /* create the inner IS */ 507 lower.i = info.xs; 508 lower.j = info.ys; 509 lower.k = info.zs; 510 upper.i = info.xs+info.xm; 511 upper.j = info.ys+info.ym; 512 upper.k = info.zs+info.zm; 513 ierr = DMDACreatePatchIS(dm,&lower,&upper,&(*iis)[i],patchis_offproc);CHKERRQ(ierr); 514 } 515 516 if (ois) { 517 /* create the outer IS */ 518 lower.i = subinfo.xs; 519 lower.j = subinfo.ys; 520 lower.k = subinfo.zs; 521 upper.i = subinfo.xs+subinfo.xm; 522 upper.j = subinfo.ys+subinfo.ym; 523 upper.k = subinfo.zs+subinfo.zm; 524 ierr = DMDACreatePatchIS(dm,&lower,&upper,&(*ois)[i],patchis_offproc);CHKERRQ(ierr); 525 } 526 } 527 PetscFunctionReturn(0); 528 } 529 530 PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm) 531 { 532 PetscErrorCode ierr; 533 DM *sdm; 534 PetscInt n,i; 535 536 PetscFunctionBegin; 537 ierr = DMDASubDomainDA_Private(dm,&n,&sdm);CHKERRQ(ierr); 538 if (names) { 539 ierr = PetscMalloc1(n,names);CHKERRQ(ierr); 540 for (i=0;i<n;i++) (*names)[i] = NULL; 541 } 542 ierr = DMDASubDomainIS_Private(dm,n,sdm,iis,ois);CHKERRQ(ierr); 543 if (subdm) *subdm = sdm; 544 else { 545 for (i=0;i<n;i++) { 546 ierr = DMDestroy(&sdm[i]);CHKERRQ(ierr); 547 } 548 } 549 if (len) *len = n; 550 PetscFunctionReturn(0); 551 } 552