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