1 #include <petsc-private/dmdaimpl.h> 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 = PetscMalloc(sizeof(PetscInt)*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 + 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, DM *dddm) 132 { 133 DM da; 134 PetscErrorCode ierr; 135 DMDALocalInfo info; 136 PetscReal lmin[3],lmax[3]; 137 PetscInt xsize,ysize,zsize; 138 PetscInt xo,yo,zo; 139 PetscInt xol,yol,zol; 140 141 PetscFunctionBegin; 142 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 143 ierr = DMDAGetOverlap(dm,&xol,&yol,&zol);CHKERRQ(ierr); 144 145 ierr = DMDACreate(PETSC_COMM_SELF,&da);CHKERRQ(ierr); 146 ierr = DMSetOptionsPrefix(da,"sub_");CHKERRQ(ierr); 147 ierr = DMDASetDim(da, info.dim);CHKERRQ(ierr); 148 ierr = DMDASetDof(da, info.dof);CHKERRQ(ierr); 149 150 ierr = DMDASetStencilType(da,info.st);CHKERRQ(ierr); 151 ierr = DMDASetStencilWidth(da,info.sw);CHKERRQ(ierr); 152 153 /* local with overlap */ 154 xsize = info.xm; 155 ysize = info.ym; 156 zsize = info.zm; 157 xo = info.xs; 158 yo = info.ys; 159 zo = info.zs; 160 if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs != 0)) { 161 xsize += xol; 162 xo -= xol; 163 } 164 if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys != 0)) { 165 ysize += yol; 166 yo -= yol; 167 } 168 if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs != 0)) { 169 zsize += zol; 170 zo -= zol; 171 } 172 173 if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs+info.xm != info.mx)) xsize += xol; 174 if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys+info.ym != info.my)) ysize += yol; 175 if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs+info.zm != info.mz)) zsize += zol; 176 177 ierr = DMDASetSizes(da, xsize, ysize, zsize);CHKERRQ(ierr); 178 ierr = DMDASetNumProcs(da, 1, 1, 1);CHKERRQ(ierr); 179 ierr = DMDASetBoundaryType(da, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED);CHKERRQ(ierr); 180 181 /* set up as a block instead */ 182 ierr = DMSetUp(da);CHKERRQ(ierr); 183 184 /* todo - nonuniform coordinates */ 185 ierr = DMDAGetLocalBoundingBox(dm,lmin,lmax);CHKERRQ(ierr); 186 ierr = DMDASetUniformCoordinates(da,lmin[0],lmax[0],lmin[1],lmax[1],lmin[2],lmax[2]);CHKERRQ(ierr); 187 188 /* nonoverlapping region */ 189 ierr = DMDASetNonOverlappingRegion(da,info.xs,info.ys,info.zs,info.xm,info.ym,info.zm);CHKERRQ(ierr); 190 191 /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */ 192 ierr = DMDASetOffset(da,xo,yo,zo,info.mx,info.my,info.mz);CHKERRQ(ierr); 193 *dddm = da; 194 PetscFunctionReturn(0); 195 } 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "DMCreateDomainDecompositionScatters_DA" 199 /* 200 Fills the local vector problem on the subdomain from the global problem. 201 202 Right now this assumes one subdomain per processor. 203 204 */ 205 PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat) 206 { 207 PetscErrorCode ierr; 208 DMDALocalInfo info,subinfo; 209 DM subdm; 210 MatStencil upper,lower; 211 IS idis,isis,odis,osis,gdis; 212 Vec svec,dvec,slvec; 213 PetscInt xm,ym,zm,xs,ys,zs; 214 215 PetscFunctionBegin; 216 if (nsubdms != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have more than one subdomain per processor (yet)"); 217 218 /* allocate the arrays of scatters */ 219 if (iscat) {ierr = PetscMalloc(sizeof(VecScatter*),iscat);CHKERRQ(ierr);} 220 if (oscat) {ierr = PetscMalloc(sizeof(VecScatter*),oscat);CHKERRQ(ierr);} 221 if (lscat) {ierr = PetscMalloc(sizeof(VecScatter*),lscat);CHKERRQ(ierr);} 222 223 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 224 subdm = subdms[0]; 225 ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr); 226 227 /* create the global and subdomain index sets for the inner domain */ 228 /* TODO - make this actually support multiple subdomains -- subdomain needs to provide where it's nonoverlapping portion belongs */ 229 ierr = DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 230 lower.i = xs; 231 lower.j = ys; 232 lower.k = zs; 233 upper.i = xs+xm; 234 upper.j = ys+ym; 235 upper.k = zs+zm; 236 ierr = DMDACreatePatchIS(dm,&lower,&upper,&idis);CHKERRQ(ierr); 237 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&isis);CHKERRQ(ierr); 238 239 /* create the global and subdomain index sets for the outer subdomain */ 240 lower.i = subinfo.xs; 241 lower.j = subinfo.ys; 242 lower.k = subinfo.zs; 243 upper.i = subinfo.xs+subinfo.xm; 244 upper.j = subinfo.ys+subinfo.ym; 245 upper.k = subinfo.zs+subinfo.zm; 246 ierr = DMDACreatePatchIS(dm,&lower,&upper,&odis);CHKERRQ(ierr); 247 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&osis);CHKERRQ(ierr); 248 249 /* global and subdomain ISes for the local indices of the subdomain */ 250 /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */ 251 lower.i = subinfo.gxs; 252 lower.j = subinfo.gys; 253 lower.k = subinfo.gzs; 254 upper.i = subinfo.gxs+subinfo.gxm; 255 upper.j = subinfo.gys+subinfo.gym; 256 upper.k = subinfo.gzs+subinfo.gzm; 257 258 ierr = DMDACreatePatchIS(dm,&lower,&upper,&gdis);CHKERRQ(ierr); 259 260 /* form the scatter */ 261 ierr = DMGetGlobalVector(dm,&dvec);CHKERRQ(ierr); 262 ierr = DMGetGlobalVector(subdm,&svec);CHKERRQ(ierr); 263 ierr = DMGetLocalVector(subdm,&slvec);CHKERRQ(ierr); 264 265 if (iscat) {ierr = VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[0]);CHKERRQ(ierr);} 266 if (oscat) {ierr = VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[0]);CHKERRQ(ierr);} 267 if (lscat) {ierr = VecScatterCreate(dvec,gdis,slvec,NULL,&(*lscat)[0]);CHKERRQ(ierr);} 268 269 ierr = DMRestoreGlobalVector(dm,&dvec);CHKERRQ(ierr); 270 ierr = DMRestoreGlobalVector(subdm,&svec);CHKERRQ(ierr); 271 ierr = DMRestoreLocalVector(subdm,&slvec);CHKERRQ(ierr); 272 273 ierr = ISDestroy(&idis);CHKERRQ(ierr); 274 ierr = ISDestroy(&isis);CHKERRQ(ierr); 275 276 ierr = ISDestroy(&odis);CHKERRQ(ierr); 277 ierr = ISDestroy(&osis);CHKERRQ(ierr); 278 279 ierr = ISDestroy(&gdis);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 #undef __FUNCT__ 284 #define __FUNCT__ "DMDASubDomainIS_Private" 285 /* We have that the interior regions are going to be the same, but the ghost regions might not match up 286 287 ---------- 288 ---------- 289 --++++++o= 290 --++++++o= 291 --++++++o= 292 --++++++o= 293 --ooooooo= 294 --======== 295 296 Therefore, for each point in the overall, we must check if it's: 297 298 1. +: In the interior of the global dm; it lines up 299 2. o: In the overlap region -- for now the same as 1; no overlap 300 3. =: In the shared ghost region -- handled by DMCreateDomainDecompositionLocalScatter() 301 4. -: In the nonshared ghost region 302 */ 303 304 PetscErrorCode DMDASubDomainIS_Private(DM dm,DM subdm,IS *iis,IS *ois) 305 { 306 PetscErrorCode ierr; 307 DMDALocalInfo info,subinfo; 308 MatStencil lower,upper; 309 310 PetscFunctionBegin; 311 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 312 ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr); 313 314 /* create the inner IS */ 315 lower.i = info.xs; 316 lower.j = info.ys; 317 lower.k = info.zs; 318 upper.i = info.xs+info.xm; 319 upper.j = info.ys+info.ym; 320 upper.k = info.zs+info.zm; 321 322 ierr = DMDACreatePatchIS(dm,&lower,&upper,iis);CHKERRQ(ierr); 323 324 /* create the outer IS */ 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,ois);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 #undef __FUNCT__ 336 #define __FUNCT__ "DMCreateDomainDecomposition_DA" 337 PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm) 338 { 339 PetscErrorCode ierr; 340 IS iis0,ois0; 341 DM subdm0; 342 343 PetscFunctionBegin; 344 if (len) *len = 1; 345 346 if (iis) {ierr = PetscMalloc(sizeof(IS*),iis);CHKERRQ(ierr);} 347 if (ois) {ierr = PetscMalloc(sizeof(IS*),ois);CHKERRQ(ierr);} 348 if (subdm) {ierr = PetscMalloc(sizeof(DM*),subdm);CHKERRQ(ierr);} 349 if (names) {ierr = PetscMalloc(sizeof(char*),names);CHKERRQ(ierr);} 350 ierr = DMDASubDomainDA_Private(dm,&subdm0);CHKERRQ(ierr); 351 ierr = DMDASubDomainIS_Private(dm,subdm0,&iis0,&ois0);CHKERRQ(ierr); 352 if (iis) (*iis)[0] = iis0; 353 else { 354 ierr = ISDestroy(&iis0);CHKERRQ(ierr); 355 } 356 if (ois) (*ois)[0] = ois0; 357 else { 358 ierr = ISDestroy(&ois0);CHKERRQ(ierr); 359 } 360 if (subdm) (*subdm)[0] = subdm0; 361 if (names) (*names)[0] = 0; 362 PetscFunctionReturn(0); 363 } 364