1 #include <petsc-private/daimpl.h> 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMDASubDomainDA_Private" 5 PetscErrorCode DMDASubDomainDA_Private(DM dm, DM *dddm) { 6 DM da; 7 DM_DA *dd; 8 PetscErrorCode ierr; 9 DMDALocalInfo info; 10 PetscReal lmin[3],lmax[3]; 11 PetscInt i,xsize,ysize,zsize; 12 PetscInt xo,yo,zo; 13 14 PetscFunctionBegin; 15 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 16 ierr = DMDACreate(PETSC_COMM_SELF,&da);CHKERRQ(ierr); 17 ierr = DMSetOptionsPrefix(da,"sub_");CHKERRQ(ierr); 18 19 ierr = DMDASetDim(da, info.dim);CHKERRQ(ierr); 20 ierr = DMDASetDof(da, info.dof);CHKERRQ(ierr); 21 22 ierr = DMDASetStencilType(da,info.st);CHKERRQ(ierr); 23 ierr = DMDASetStencilWidth(da,info.sw);CHKERRQ(ierr); 24 25 dd = (DM_DA *)dm->data; 26 27 /* local with overlap */ 28 xsize = info.xm; 29 ysize = info.ym; 30 zsize = info.zm; 31 xo = info.xs; 32 yo = info.ys; 33 zo = info.zs; 34 if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs != 0)) { 35 xsize += dd->overlap; 36 xo -= dd->overlap; 37 } 38 if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys != 0)) { 39 ysize += dd->overlap; 40 yo -= dd->overlap; 41 } 42 if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs != 0)) { 43 zsize += dd->overlap; 44 zo -= dd->overlap; 45 } 46 47 if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs+info.xm != info.mx)) { 48 xsize += dd->overlap; 49 } 50 if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys+info.ym != info.my)) { 51 ysize += dd->overlap; 52 } 53 if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs+info.zm != info.mz)) { 54 zsize += dd->overlap; 55 } 56 57 ierr = DMDASetSizes(da, xsize, ysize, zsize);CHKERRQ(ierr); 58 ierr = DMDASetNumProcs(da, 1, 1, 1);CHKERRQ(ierr); 59 ierr = DMDASetBoundaryType(da, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED);CHKERRQ(ierr); 60 61 /* set up as a block instead */ 62 ierr = DMSetUp(da);CHKERRQ(ierr); 63 ierr = DMDASetOffset(da,xo,yo,zo);CHKERRQ(ierr); 64 65 66 /* todo - nonuniform coordinates */ 67 ierr = DMDAGetLocalBoundingBox(dm,lmin,lmax);CHKERRQ(ierr); 68 for (i=info.dim; i<3; i++) {lmin[i] = 0; lmax[i] = 0;} 69 ierr = DMDASetUniformCoordinates(da,lmin[0],lmax[0],lmin[1],lmax[1],lmin[2],lmax[2]);CHKERRQ(ierr); 70 71 dd = (DM_DA *)da->data; 72 dd->Mo = info.mx; 73 dd->No = info.my; 74 dd->Po = info.mz; 75 76 *dddm = da; 77 78 PetscFunctionReturn(0); 79 } 80 81 #undef __FUNCT__ 82 #define __FUNCT__ "DMCreateDomainDecompositionScatters_DA" 83 /* 84 Fills the local vector problem on the subdomain from the global problem. 85 86 */ 87 PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat) { 88 PetscErrorCode ierr; 89 DMDALocalInfo dinfo,sinfo; 90 IS isis,idis,osis,odis,gsis,gdis; 91 PetscInt *ididx,*isidx,*odidx,*osidx,*gdidx,*gsidx,*idx_global,n_global,*idx_sub,n_sub; 92 PetscInt l,i,j,k,d,n_i,n_o,n_g,sl,dl,di,dj,dk,si,sj,sk; 93 Vec dvec,svec,slvec; 94 DM subdm; 95 96 PetscFunctionBegin; 97 98 /* allocate the arrays of scatters */ 99 if (iscat) {ierr = PetscMalloc(sizeof(VecScatter *),iscat);CHKERRQ(ierr);} 100 if (oscat) {ierr = PetscMalloc(sizeof(VecScatter *),oscat);CHKERRQ(ierr);} 101 if (lscat) {ierr = PetscMalloc(sizeof(VecScatter *),lscat);CHKERRQ(ierr);} 102 103 ierr = DMDAGetLocalInfo(dm,&dinfo);CHKERRQ(ierr); 104 ierr = DMDAGetGlobalIndices(dm,&n_global,&idx_global);CHKERRQ(ierr); 105 for (l = 0;l < nsubdms;l++) { 106 n_i = 0; 107 n_o = 0; 108 n_g = 0; 109 subdm = subdms[l]; 110 ierr = DMDAGetLocalInfo(subdm,&sinfo);CHKERRQ(ierr); 111 ierr = DMDAGetGlobalIndices(subdm,&n_sub,&idx_sub);CHKERRQ(ierr); 112 /* count the three region sizes */ 113 for (k=sinfo.gzs;k<sinfo.gzs+sinfo.gzm;k++) { 114 for (j=sinfo.gys;j<sinfo.gys+sinfo.gym;j++) { 115 for (i=sinfo.gxs;i<sinfo.gxs+sinfo.gxm;i++) { 116 for (d=0;d<sinfo.dof;d++) { 117 if (k >= sinfo.zs && j >= sinfo.ys && i >= sinfo.xs && 118 k < sinfo.zs+sinfo.zm && j < sinfo.ys+sinfo.ym && i < sinfo.xs+sinfo.xm) { 119 120 /* interior - subinterior overlap */ 121 if (k >= dinfo.zs && j >= dinfo.ys && i >= dinfo.xs && 122 k < dinfo.zs+dinfo.zm && j < dinfo.ys+dinfo.ym && i < dinfo.xs+dinfo.xm) { 123 n_i++; 124 } 125 /* ghost - subinterior overlap */ 126 if (k >= dinfo.gzs && j >= dinfo.gys && i >= dinfo.gxs && 127 k < dinfo.gzs+dinfo.gzm && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) { 128 n_o++; 129 } 130 } 131 132 /* ghost - subghost overlap */ 133 if (k >= dinfo.gzs && j >= dinfo.gys && i >= dinfo.gxs && 134 k < dinfo.gzs+dinfo.gzm && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) { 135 n_g++; 136 } 137 } 138 } 139 } 140 } 141 142 if (n_g == 0) SETERRQ(((PetscObject)subdm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Processor-local domain and subdomain do not intersect!"); 143 144 /* local and subdomain local index set indices */ 145 ierr = PetscMalloc(sizeof(PetscInt)*n_i,&ididx);CHKERRQ(ierr); 146 ierr = PetscMalloc(sizeof(PetscInt)*n_i,&isidx);CHKERRQ(ierr); 147 148 ierr = PetscMalloc(sizeof(PetscInt)*n_o,&odidx);CHKERRQ(ierr); 149 ierr = PetscMalloc(sizeof(PetscInt)*n_o,&osidx);CHKERRQ(ierr); 150 151 ierr = PetscMalloc(sizeof(PetscInt)*n_g,&gdidx);CHKERRQ(ierr); 152 ierr = PetscMalloc(sizeof(PetscInt)*n_g,&gsidx);CHKERRQ(ierr); 153 154 n_i = 0; n_o = 0;n_g = 0; 155 for (k=sinfo.gzs;k<sinfo.gzs+sinfo.gzm;k++) { 156 for (j=sinfo.gys;j<sinfo.gys+sinfo.gym;j++) { 157 for (i=sinfo.gxs;i<sinfo.gxs+sinfo.gxm;i++) { 158 for (d=0;d<sinfo.dof;d++) { 159 si = i - sinfo.gxs; 160 sj = j - sinfo.gys; 161 sk = k - sinfo.gzs; 162 sl = d + sinfo.dof*(si + sinfo.gxm*(sj + sinfo.gym*sk)); 163 di = i - dinfo.gxs; 164 dj = j - dinfo.gys; 165 dk = k - dinfo.gzs; 166 dl = d + dinfo.dof*(di + dinfo.gxm*(dj + dinfo.gym*dk)); 167 168 if (k >= sinfo.zs && j >= sinfo.ys && i >= sinfo.xs && 169 k < sinfo.zs+sinfo.zm && j < sinfo.ys+sinfo.ym && i < sinfo.xs+sinfo.xm) { 170 171 /* interior - subinterior overlap */ 172 if (k >= dinfo.zs && j >= dinfo.ys && i >= dinfo.xs && 173 k < dinfo.zs+dinfo.zm && j < dinfo.ys+dinfo.ym && i < dinfo.xs+dinfo.xm) { 174 ididx[n_i] = idx_global[dl]; 175 isidx[n_i] = idx_sub[sl]; 176 n_i++; 177 } 178 /* ghost - subinterior overlap */ 179 if (k >= dinfo.gzs && j >= dinfo.gys && i >= dinfo.gxs && 180 k < dinfo.gzs+dinfo.gzm && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) { 181 odidx[n_o] = idx_global[dl]; 182 osidx[n_o] = idx_sub[sl]; 183 n_o++; 184 } 185 } 186 187 /* ghost - subghost overlap */ 188 if (k >= dinfo.gzs && j >= dinfo.gys && i >= dinfo.gxs && 189 k < dinfo.gzs+dinfo.gzm && j < dinfo.gys+dinfo.gym && i < dinfo.gxs+dinfo.gxm) { 190 gdidx[n_g] = idx_global[dl]; 191 gsidx[n_g] = sl; 192 n_g++; 193 } 194 } 195 } 196 } 197 } 198 199 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_i,ididx,PETSC_OWN_POINTER,&idis);CHKERRQ(ierr); 200 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_i,isidx,PETSC_OWN_POINTER,&isis);CHKERRQ(ierr); 201 202 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_o,odidx,PETSC_OWN_POINTER,&odis);CHKERRQ(ierr); 203 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_o,osidx,PETSC_OWN_POINTER,&osis);CHKERRQ(ierr); 204 205 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_g,gdidx,PETSC_OWN_POINTER,&gdis);CHKERRQ(ierr); 206 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_g,gsidx,PETSC_OWN_POINTER,&gsis);CHKERRQ(ierr); 207 208 /* form the scatter */ 209 ierr = DMGetGlobalVector(dm,&dvec);CHKERRQ(ierr); 210 ierr = DMGetGlobalVector(subdm,&svec);CHKERRQ(ierr); 211 ierr = DMGetLocalVector(subdm,&slvec);CHKERRQ(ierr); 212 213 if (iscat) {ierr = VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[l]);CHKERRQ(ierr);} 214 if (oscat) {ierr = VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[l]);CHKERRQ(ierr);} 215 if (lscat) {ierr = VecScatterCreate(dvec,gdis,slvec,gsis,&(*lscat)[l]);CHKERRQ(ierr);} 216 217 ierr = DMRestoreGlobalVector(dm,&dvec);CHKERRQ(ierr); 218 ierr = DMRestoreGlobalVector(subdm,&svec);CHKERRQ(ierr); 219 ierr = DMRestoreLocalVector(subdm,&slvec);CHKERRQ(ierr); 220 221 ierr = ISDestroy(&idis);CHKERRQ(ierr); 222 ierr = ISDestroy(&isis);CHKERRQ(ierr); 223 224 ierr = ISDestroy(&odis);CHKERRQ(ierr); 225 ierr = ISDestroy(&osis);CHKERRQ(ierr); 226 227 ierr = ISDestroy(&gdis);CHKERRQ(ierr); 228 ierr = ISDestroy(&gsis);CHKERRQ(ierr); 229 } 230 PetscFunctionReturn(0); 231 232 } 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "DMDASubDomainIS_Private" 236 /* We have that the interior regions are going to be the same, but the ghost regions might not match up 237 238 ---------- 239 ---------- 240 --++++++o= 241 --++++++o= 242 --++++++o= 243 --++++++o= 244 --ooooooo= 245 --======== 246 247 Therefore, for each point in the overall, we must check if it's: 248 249 1. +: In the interior of the global dm; it lines up 250 2. o: In the overlap region -- for now the same as 1; no overlap 251 3. =: In the shared ghost region -- handled by DMCreateDomainDecompositionLocalScatter() 252 4. -: In the nonshared ghost region 253 */ 254 255 PetscErrorCode DMDASubDomainIS_Private(DM dm,DM subdm,IS *iis,IS *ois) { 256 PetscErrorCode ierr; 257 DMDALocalInfo info,subinfo; 258 PetscInt *iiidx,*oiidx,*gidx,gindx; 259 PetscInt i,j,k,d,n,nsub,nover,llindx,lindx,li,lj,lk,gi,gj,gk; 260 261 PetscFunctionBegin; 262 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 263 ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr); 264 ierr = DMDAGetGlobalIndices(dm,&n,&gidx);CHKERRQ(ierr); 265 /* todo -- overlap */ 266 nsub = info.xm*info.ym*info.zm*info.dof; 267 nover = subinfo.xm*subinfo.ym*subinfo.zm*subinfo.dof; 268 /* iis is going to have size of the local problem's global part but have a lot of fill-in */ 269 ierr = PetscMalloc(sizeof(PetscInt)*nsub,&iiidx);CHKERRQ(ierr); 270 /* ois is going to have size of the local problem's global part */ 271 ierr = PetscMalloc(sizeof(PetscInt)*nover,&oiidx);CHKERRQ(ierr); 272 /* loop over the ghost region of the subdm and fill in the indices */ 273 for (k=subinfo.gzs;k<subinfo.gzs+subinfo.gzm;k++) { 274 for (j=subinfo.gys;j<subinfo.gys+subinfo.gym;j++) { 275 for (i=subinfo.gxs;i<subinfo.gxs+subinfo.gxm;i++) { 276 for (d=0;d<subinfo.dof;d++) { 277 li = i - subinfo.xs; 278 lj = j - subinfo.ys; 279 lk = k - subinfo.zs; 280 lindx = d + subinfo.dof*(li + subinfo.xm*(lj + subinfo.ym*lk)); 281 li = i - info.xs; 282 lj = j - info.ys; 283 lk = k - info.zs; 284 llindx = d + info.dof*(li + info.xm*(lj + info.ym*lk)); 285 gi = i - info.gxs; 286 gj = j - info.gys; 287 gk = k - info.gzs; 288 gindx = d + info.dof*(gi + info.gxm*(gj + info.gym*gk)); 289 290 /* check if the current point is inside the interior region */ 291 if (k >= info.zs && j >= info.ys && i >= info.xs && 292 k < info.zs+info.zm && j < info.ys+info.ym && i < info.xs+info.xm) { 293 iiidx[llindx] = gidx[gindx]; 294 oiidx[lindx] = gidx[gindx]; 295 /* overlap region */ 296 } else if (k >= subinfo.zs && j >= subinfo.ys && i >= subinfo.xs && 297 k < subinfo.zs+subinfo.zm && j < subinfo.ys+subinfo.ym && i < subinfo.xs+subinfo.xm) { 298 oiidx[lindx] = gidx[gindx]; 299 } 300 } 301 } 302 } 303 } 304 305 /* create the index sets */ 306 ierr = ISCreateGeneral(PETSC_COMM_SELF,nsub,iiidx,PETSC_OWN_POINTER,iis);CHKERRQ(ierr); 307 ierr = ISCreateGeneral(PETSC_COMM_SELF,nover,oiidx,PETSC_OWN_POINTER,ois);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "DMCreateDomainDecomposition_DA" 313 PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm) { 314 PetscErrorCode ierr; 315 IS iis0,ois0; 316 DM subdm0; 317 PetscFunctionBegin; 318 if (len)*len = 1; 319 320 if (iis) {ierr = PetscMalloc(sizeof(IS *),iis);CHKERRQ(ierr);} 321 if (ois) {ierr = PetscMalloc(sizeof(IS *),ois);CHKERRQ(ierr);} 322 if (subdm) {ierr = PetscMalloc(sizeof(DM *),subdm);CHKERRQ(ierr);} 323 if (names) {ierr = PetscMalloc(sizeof(char *),names);CHKERRQ(ierr);} 324 ierr = DMDASubDomainDA_Private(dm,&subdm0);CHKERRQ(ierr); 325 ierr = DMDASubDomainIS_Private(dm,subdm0,&iis0,&ois0);CHKERRQ(ierr); 326 if (iis) { 327 (*iis)[0] = iis0; 328 } else { 329 ierr = ISDestroy(&iis0);CHKERRQ(ierr); 330 } 331 if (ois) { 332 (*ois)[0] = ois0; 333 } else { 334 ierr = ISDestroy(&ois0);CHKERRQ(ierr); 335 } 336 if (subdm) (*subdm)[0] = subdm0; 337 if (names) (*names)[0] = 0; 338 PetscFunctionReturn(0); 339 } 340