1 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMDACreatePatchIS" 5 /*@ 6 DMDACreatePatchIS - Creates an index set corresponding to a patch of the DA. 7 8 Not Collective 9 10 Input Parameters: 11 + da - the DMDA 12 . lower - a matstencil with i, j and k corresponding to the lower corner of the patch 13 - upper - a matstencil with i, j and k corresponding to the upper corner of the patch 14 15 Output Parameters: 16 . is - the IS corresponding to the patch 17 18 Level: developer 19 20 .seealso: DMDACreateDomainDecomposition(), DMDACreateDomainDecompositionScatters() 21 @*/ 22 PetscErrorCode DMDACreatePatchIS(DM da,MatStencil *lower,MatStencil *upper,IS *is) 23 { 24 PetscInt ms=0,ns=0,ps=0; 25 PetscInt me=1,ne=1,pe=1; 26 PetscInt mr=0,nr=0,pr=0; 27 PetscInt ii,jj,kk; 28 PetscInt si,sj,sk; 29 PetscInt i,j,k,l,idx; 30 PetscInt base; 31 PetscInt xm=1,ym=1,zm=1; 32 const PetscInt *lx,*ly,*lz; 33 PetscInt ox,oy,oz; 34 PetscInt m,n,p,M,N,P,dof; 35 PetscInt nindices; 36 PetscInt *indices; 37 DM_DA *dd = (DM_DA*)da->data; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 /* need to get the sizes of the actual DM rather than the "global" space of a subdomain DM */ 42 M = dd->M;N = dd->N;P=dd->P; 43 m = dd->m;n = dd->n;p=dd->p; 44 dof = dd->w; 45 ierr = DMDAGetOffset(da,&ox,&oy,&oz,NULL,NULL,NULL);CHKERRQ(ierr); 46 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,&lz);CHKERRQ(ierr); 47 nindices = (upper->i - lower->i)*(upper->j - lower->j)*(upper->k - lower->k)*dof; 48 ierr = PetscMalloc1(nindices,&indices);CHKERRQ(ierr); 49 /* start at index 0 on processor 0 */ 50 mr = 0; 51 nr = 0; 52 pr = 0; 53 ms = 0; 54 ns = 0; 55 ps = 0; 56 if (lx) me = lx[0]; 57 if (ly) ne = ly[0]; 58 if (lz) pe = lz[0]; 59 idx = 0; 60 for (k=lower->k-oz;k<upper->k-oz;k++) { 61 for (j=lower->j-oy;j < upper->j-oy;j++) { 62 for (i=lower->i-ox;i < upper->i-ox;i++) { 63 /* "actual" indices rather than ones outside of the domain */ 64 ii = i; 65 jj = j; 66 kk = k; 67 if (ii < 0) ii = M + ii; 68 if (jj < 0) jj = N + jj; 69 if (kk < 0) kk = P + kk; 70 if (ii > M-1) ii = ii - M; 71 if (jj > N-1) jj = jj - N; 72 if (kk > P-1) kk = kk - P; 73 /* gone out of processor range on x axis */ 74 while(ii > me-1 || ii < ms) { 75 if (mr == m-1) { 76 ms = 0; 77 me = lx[0]; 78 mr = 0; 79 } else { 80 mr++; 81 ms = me; 82 me += lx[mr]; 83 } 84 } 85 /* gone out of processor range on y axis */ 86 while(jj > ne-1 || jj < ns) { 87 if (nr == n-1) { 88 ns = 0; 89 ne = ly[0]; 90 nr = 0; 91 } else { 92 nr++; 93 ns = ne; 94 ne += ly[nr]; 95 } 96 } 97 /* gone out of processor range on z axis */ 98 while(kk > pe-1 || kk < ps) { 99 if (pr == p-1) { 100 ps = 0; 101 pe = lz[0]; 102 pr = 0; 103 } else { 104 pr++; 105 ps = pe; 106 pe += lz[pr]; 107 } 108 } 109 /* compute the vector base on owning processor */ 110 xm = me - ms; 111 ym = ne - ns; 112 zm = pe - ps; 113 base = ms*ym*zm + ns*M*zm + ps*M*N; 114 /* compute the local coordinates on owning processor */ 115 si = ii - ms; 116 sj = jj - ns; 117 sk = kk - ps; 118 for (l=0;l<dof;l++) { 119 indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk); 120 idx++; 121 } 122 } 123 } 124 } 125 ISCreateGeneral(PETSC_COMM_SELF,idx,indices,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "DMDASubDomainDA_Private" 131 PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm) 132 { 133 DM *da; 134 PetscInt dim,size,i,j,k,idx; 135 PetscErrorCode ierr; 136 DMDALocalInfo info; 137 PetscInt xsize,ysize,zsize; 138 PetscInt xo,yo,zo; 139 PetscInt xs,ys,zs; 140 PetscInt xm=1,ym=1,zm=1; 141 PetscInt xol,yol,zol; 142 PetscInt m=1,n=1,p=1; 143 PetscInt M,N,P; 144 PetscInt pm,mtmp; 145 146 PetscFunctionBegin; 147 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 148 ierr = DMDAGetOverlap(dm,&xol,&yol,&zol);CHKERRQ(ierr); 149 ierr = DMDAGetNumLocalSubDomains(dm,&size);CHKERRQ(ierr); 150 ierr = PetscMalloc1(size,&da);CHKERRQ(ierr); 151 152 if (nlocal) *nlocal = size; 153 154 dim = info.dim; 155 156 M = info.xm; 157 N = info.ym; 158 P = info.zm; 159 160 if (dim == 1) { 161 m = size; 162 } else if (dim == 2) { 163 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N))); 164 while (m > 0) { 165 n = size/m; 166 if (m*n*p == size) break; 167 m--; 168 } 169 } else if (dim == 3) { 170 n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.))); if (!n) n = 1; 171 while (n > 0) { 172 pm = size/n; 173 if (n*pm == size) break; 174 n--; 175 } 176 if (!n) n = 1; 177 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n))); 178 if (!m) m = 1; 179 while (m > 0) { 180 p = size/(m*n); 181 if (m*n*p == size) break; 182 m--; 183 } 184 if (M > P && m < p) {mtmp = m; m = p; p = mtmp;} 185 } 186 187 zs = info.zs; 188 idx = 0; 189 for (k = 0; k < p; k++) { 190 ys = info.ys; 191 for (j = 0; j < n; j++) { 192 xs = info.xs; 193 for (i = 0; i < m; i++) { 194 if (dim == 1) { 195 xm = M/m + ((M % m) > i); 196 } else if (dim == 2) { 197 xm = M/m + ((M % m) > i); 198 ym = N/n + ((N % n) > j); 199 } else if (dim == 3) { 200 xm = M/m + ((M % m) > i); 201 ym = N/n + ((N % n) > j); 202 zm = P/p + ((P % p) > k); 203 } 204 205 xsize = xm; 206 ysize = ym; 207 zsize = zm; 208 xo = xs; 209 yo = ys; 210 zo = zs; 211 212 ierr = DMDACreate(PETSC_COMM_SELF,&(da[idx]));CHKERRQ(ierr); 213 ierr = DMSetOptionsPrefix(da[idx],"sub_");CHKERRQ(ierr); 214 ierr = DMSetDimension(da[idx], info.dim);CHKERRQ(ierr); 215 ierr = DMDASetDof(da[idx], info.dof);CHKERRQ(ierr); 216 217 ierr = DMDASetStencilType(da[idx],info.st);CHKERRQ(ierr); 218 ierr = DMDASetStencilWidth(da[idx],info.sw);CHKERRQ(ierr); 219 220 if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) { 221 xsize += xol; 222 xo -= xol; 223 } 224 if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) { 225 ysize += yol; 226 yo -= yol; 227 } 228 if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) { 229 zsize += zol; 230 zo -= zol; 231 } 232 233 if (info.bx == DM_BOUNDARY_PERIODIC || (xs+xm != info.mx)) xsize += xol; 234 if (info.by == DM_BOUNDARY_PERIODIC || (ys+ym != info.my)) ysize += yol; 235 if (info.bz == DM_BOUNDARY_PERIODIC || (zs+zm != info.mz)) zsize += zol; 236 237 if (info.bx != DM_BOUNDARY_PERIODIC) { 238 if (xo < 0) { 239 xsize += xo; 240 xo = 0; 241 } 242 if (xo+xsize > info.mx-1) { 243 xsize -= xo+xsize - info.mx; 244 } 245 } 246 if (info.by != DM_BOUNDARY_PERIODIC) { 247 if (yo < 0) { 248 ysize += yo; 249 yo = 0; 250 } 251 if (yo+ysize > info.my-1) { 252 ysize -= yo+ysize - info.my; 253 } 254 } 255 if (info.bz != DM_BOUNDARY_PERIODIC) { 256 if (zo < 0) { 257 zsize += zo; 258 zo = 0; 259 } 260 if (zo+zsize > info.mz-1) { 261 zsize -= zo+zsize - info.mz; 262 } 263 } 264 265 ierr = DMDASetSizes(da[idx], xsize, ysize, zsize);CHKERRQ(ierr); 266 ierr = DMDASetNumProcs(da[idx], 1, 1, 1);CHKERRQ(ierr); 267 ierr = DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);CHKERRQ(ierr); 268 269 /* set up as a block instead */ 270 ierr = DMSetUp(da[idx]);CHKERRQ(ierr); 271 272 /* nonoverlapping region */ 273 ierr = DMDASetNonOverlappingRegion(da[idx],xs,ys,zs,xm,ym,zm);CHKERRQ(ierr); 274 275 /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */ 276 ierr = DMDASetOffset(da[idx],xo,yo,zo,info.mx,info.my,info.mz);CHKERRQ(ierr); 277 xs += xm; 278 idx++; 279 } 280 ys += ym; 281 } 282 zs += zm; 283 } 284 *sdm = da; 285 PetscFunctionReturn(0); 286 } 287 288 #undef __FUNCT__ 289 #define __FUNCT__ "DMCreateDomainDecompositionScatters_DA" 290 /* 291 Fills the local vector problem on the subdomain from the global problem. 292 293 Right now this assumes one subdomain per processor. 294 295 */ 296 PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat) 297 { 298 PetscErrorCode ierr; 299 DMDALocalInfo info,subinfo; 300 DM subdm; 301 MatStencil upper,lower; 302 IS idis,isis,odis,osis,gdis; 303 Vec svec,dvec,slvec; 304 PetscInt xm,ym,zm,xs,ys,zs; 305 PetscInt i; 306 307 PetscFunctionBegin; 308 309 /* allocate the arrays of scatters */ 310 if (iscat) {ierr = PetscMalloc1(nsubdms,iscat);CHKERRQ(ierr);} 311 if (oscat) {ierr = PetscMalloc1(nsubdms,oscat);CHKERRQ(ierr);} 312 if (lscat) {ierr = PetscMalloc1(nsubdms,lscat);CHKERRQ(ierr);} 313 314 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 315 for (i = 0; i < nsubdms; i++) { 316 subdm = subdms[i]; 317 ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr); 318 ierr = DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 319 320 /* create the global and subdomain index sets for the inner domain */ 321 lower.i = xs; 322 lower.j = ys; 323 lower.k = zs; 324 upper.i = xs+xm; 325 upper.j = ys+ym; 326 upper.k = zs+zm; 327 ierr = DMDACreatePatchIS(dm,&lower,&upper,&idis);CHKERRQ(ierr); 328 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&isis);CHKERRQ(ierr); 329 330 /* create the global and subdomain index sets for the outer subdomain */ 331 lower.i = subinfo.xs; 332 lower.j = subinfo.ys; 333 lower.k = subinfo.zs; 334 upper.i = subinfo.xs+subinfo.xm; 335 upper.j = subinfo.ys+subinfo.ym; 336 upper.k = subinfo.zs+subinfo.zm; 337 ierr = DMDACreatePatchIS(dm,&lower,&upper,&odis);CHKERRQ(ierr); 338 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&osis);CHKERRQ(ierr); 339 340 /* global and subdomain ISes for the local indices of the subdomain */ 341 /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */ 342 lower.i = subinfo.gxs; 343 lower.j = subinfo.gys; 344 lower.k = subinfo.gzs; 345 upper.i = subinfo.gxs+subinfo.gxm; 346 upper.j = subinfo.gys+subinfo.gym; 347 upper.k = subinfo.gzs+subinfo.gzm; 348 349 ierr = DMDACreatePatchIS(dm,&lower,&upper,&gdis);CHKERRQ(ierr); 350 351 /* form the scatter */ 352 ierr = DMGetGlobalVector(dm,&dvec);CHKERRQ(ierr); 353 ierr = DMGetGlobalVector(subdm,&svec);CHKERRQ(ierr); 354 ierr = DMGetLocalVector(subdm,&slvec);CHKERRQ(ierr); 355 356 if (iscat) {ierr = VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[i]);CHKERRQ(ierr);} 357 if (oscat) {ierr = VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[i]);CHKERRQ(ierr);} 358 if (lscat) {ierr = VecScatterCreate(dvec,gdis,slvec,NULL,&(*lscat)[i]);CHKERRQ(ierr);} 359 360 ierr = DMRestoreGlobalVector(dm,&dvec);CHKERRQ(ierr); 361 ierr = DMRestoreGlobalVector(subdm,&svec);CHKERRQ(ierr); 362 ierr = DMRestoreLocalVector(subdm,&slvec);CHKERRQ(ierr); 363 364 ierr = ISDestroy(&idis);CHKERRQ(ierr); 365 ierr = ISDestroy(&isis);CHKERRQ(ierr); 366 367 ierr = ISDestroy(&odis);CHKERRQ(ierr); 368 ierr = ISDestroy(&osis);CHKERRQ(ierr); 369 370 ierr = ISDestroy(&gdis);CHKERRQ(ierr); 371 } 372 PetscFunctionReturn(0); 373 } 374 375 #undef __FUNCT__ 376 #define __FUNCT__ "DMDASubDomainIS_Private" 377 PetscErrorCode DMDASubDomainIS_Private(DM dm,PetscInt n,DM *subdm,IS **iis,IS **ois) 378 { 379 PetscErrorCode ierr; 380 PetscInt i; 381 DMDALocalInfo info,subinfo; 382 MatStencil lower,upper; 383 384 PetscFunctionBegin; 385 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 386 if (iis) {ierr = PetscMalloc1(n,iis);CHKERRQ(ierr);} 387 if (ois) {ierr = PetscMalloc1(n,ois);CHKERRQ(ierr);} 388 389 for (i = 0;i < n; i++) { 390 ierr = DMDAGetLocalInfo(subdm[i],&subinfo);CHKERRQ(ierr); 391 if (iis) { 392 /* create the inner IS */ 393 lower.i = info.xs; 394 lower.j = info.ys; 395 lower.k = info.zs; 396 upper.i = info.xs+info.xm; 397 upper.j = info.ys+info.ym; 398 upper.k = info.zs+info.zm; 399 ierr = DMDACreatePatchIS(dm,&lower,&upper,&(*iis)[i]);CHKERRQ(ierr); 400 } 401 402 if (ois) { 403 /* create the outer IS */ 404 lower.i = subinfo.xs; 405 lower.j = subinfo.ys; 406 lower.k = subinfo.zs; 407 upper.i = subinfo.xs+subinfo.xm; 408 upper.j = subinfo.ys+subinfo.ym; 409 upper.k = subinfo.zs+subinfo.zm; 410 ierr = DMDACreatePatchIS(dm,&lower,&upper,&(*ois)[i]);CHKERRQ(ierr); 411 } 412 } 413 PetscFunctionReturn(0); 414 } 415 416 #undef __FUNCT__ 417 #define __FUNCT__ "DMCreateDomainDecomposition_DA" 418 PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm) 419 { 420 PetscErrorCode ierr; 421 DM *sdm; 422 PetscInt n,i; 423 424 PetscFunctionBegin; 425 ierr = DMDASubDomainDA_Private(dm,&n,&sdm);CHKERRQ(ierr); 426 if (names) { 427 ierr = PetscMalloc1(n,names);CHKERRQ(ierr); 428 for (i=0;i<n;i++) (*names)[i] = 0; 429 } 430 ierr = DMDASubDomainIS_Private(dm,n,sdm,iis,ois);CHKERRQ(ierr); 431 if (subdm) *subdm = sdm; 432 else { 433 for (i=0;i<n;i++) { 434 ierr = DMDestroy(&sdm[i]);CHKERRQ(ierr); 435 } 436 } 437 if (len) *len = n; 438 PetscFunctionReturn(0); 439 } 440