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