1 #include <petsc-private/dmdaimpl.h> 2 #include <petscsf.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMDACreatePatchIS" 6 7 PetscErrorCode DMDACreatePatchIS(DM da,MatStencil *lower,MatStencil *upper,IS *is) 8 { 9 PetscErrorCode ierr; 10 PetscInt i,j,k,idx; 11 PetscInt ii,jj,kk; 12 Vec v; 13 PetscInt n,pn,bs; 14 PetscMPIInt rank; 15 PetscSF sf,psf; 16 PetscLayout map; 17 MPI_Comm comm; 18 PetscInt *natidx,*globidx,*leafidx; 19 PetscInt *pnatidx,*pleafidx; 20 PetscInt base; 21 PetscInt ox,oy,oz; 22 DM_DA *dd; 23 24 PetscFunctionBegin; 25 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 26 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 27 dd = (DM_DA*)da->data; 28 29 /* construct the natural mapping */ 30 ierr = DMGetGlobalVector(da,&v);CHKERRQ(ierr); 31 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 32 ierr = VecGetOwnershipRange(v,&base,NULL);CHKERRQ(ierr); 33 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr); 34 ierr = DMRestoreGlobalVector(da,&v);CHKERRQ(ierr); 35 36 /* construct the layout */ 37 ierr = PetscLayoutCreate(comm,&map);CHKERRQ(ierr); 38 ierr = PetscLayoutSetBlockSize(map,1);CHKERRQ(ierr); 39 ierr = PetscLayoutSetLocalSize(map,n);CHKERRQ(ierr); 40 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 41 42 /* construct the list of natural indices on this process when PETSc ordering is considered */ 43 ierr = DMDAGetOffset(da,&ox,&oy,&oz,NULL,NULL,NULL);CHKERRQ(ierr); 44 ierr = PetscMalloc(sizeof(PetscInt)*n,&natidx);CHKERRQ(ierr); 45 ierr = PetscMalloc(sizeof(PetscInt)*n,&globidx);CHKERRQ(ierr); 46 ierr = PetscMalloc(sizeof(PetscInt)*n,&leafidx);CHKERRQ(ierr); 47 idx = 0; 48 for (k=dd->zs; k<dd->ze; k++) { 49 for (j=dd->ys; j<dd->ye; j++) { 50 for (i=dd->xs; i<dd->xe; i++) { 51 natidx[idx] = i + dd->w*(j*dd->M + k*dd->M*dd->N); 52 globidx[idx] = base + idx; 53 leafidx[idx] = 0; 54 idx++; 55 } 56 } 57 } 58 59 if (idx != n) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE, "for some reason the count is wrong."); 60 61 /* construct the SF going from the natural indices to the local set of PETSc indices */ 62 ierr = PetscSFCreate(comm,&sf);CHKERRQ(ierr); 63 ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr); 64 ierr = PetscSFSetGraphLayout(sf,map,n,NULL,PETSC_OWN_POINTER,natidx);CHKERRQ(ierr); 65 66 /* broadcast the global indices over to the corresponding natural indices */ 67 ierr = PetscSFGatherBegin(sf,MPIU_INT,globidx,leafidx);CHKERRQ(ierr); 68 ierr = PetscSFGatherEnd(sf,MPIU_INT,globidx,leafidx);CHKERRQ(ierr); 69 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 70 71 72 pn = dd->w*(upper->k - lower->k)*(upper->j - lower->j)*(upper->i - lower->i); 73 ierr = PetscMalloc(sizeof(PetscInt)*pn,&pnatidx);CHKERRQ(ierr); 74 ierr = PetscMalloc(sizeof(PetscInt)*pn,&pleafidx);CHKERRQ(ierr); 75 idx = 0; 76 for (k=lower->k-oz; k<upper->k-oz; k++) { 77 for (j=lower->j-oy; j<upper->j-oy; j++) { 78 for (i=dd->w*(lower->i-ox); i<dd->w*(upper->i-ox); i++) { 79 ii = i % (dd->w*dd->M); 80 jj = j % dd->N; 81 kk = k % dd->P; 82 if (ii < 0) ii = dd->w*dd->M + ii; 83 if (jj < 0) jj = dd->N + jj; 84 if (kk < 0) kk = dd->P + kk; 85 pnatidx[idx] = ii + dd->w*(jj*dd->M + kk*dd->M*dd->N); 86 idx++; 87 } 88 } 89 } 90 91 if (idx != pn) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE, "for some reason the count is wrong"); 92 93 ierr = PetscSFCreate(comm,&psf);CHKERRQ(ierr); 94 ierr = PetscSFSetFromOptions(psf);CHKERRQ(ierr); 95 ierr = PetscSFSetGraphLayout(psf,map,pn,NULL,PETSC_OWN_POINTER,pnatidx);CHKERRQ(ierr); 96 97 /* broadcast the global indices through to the patch */ 98 ierr = PetscSFBcastBegin(psf,MPIU_INT,leafidx,pleafidx);CHKERRQ(ierr); 99 ierr = PetscSFBcastEnd(psf,MPIU_INT,leafidx,pleafidx);CHKERRQ(ierr); 100 101 /* create the IS */ 102 ierr = ISCreateGeneral(comm,pn,pleafidx,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 103 104 ierr = PetscSFDestroy(&psf);CHKERRQ(ierr); 105 106 ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr); 107 108 ierr = PetscFree(globidx);CHKERRQ(ierr); 109 ierr = PetscFree(leafidx);CHKERRQ(ierr); 110 ierr = PetscFree(natidx);CHKERRQ(ierr); 111 ierr = PetscFree(pnatidx);CHKERRQ(ierr); 112 PetscFunctionReturn(0); 113 } 114 115 #undef __FUNCT__ 116 #define __FUNCT__ "DMDASubDomainDA_Private" 117 PetscErrorCode DMDASubDomainDA_Private(DM dm, DM *dddm) 118 { 119 DM da; 120 PetscErrorCode ierr; 121 DMDALocalInfo info; 122 PetscReal lmin[3],lmax[3]; 123 PetscInt xsize,ysize,zsize; 124 PetscInt xo,yo,zo; 125 PetscInt xol,yol,zol; 126 127 PetscFunctionBegin; 128 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 129 ierr = DMDAGetOverlap(dm,&xol,&yol,&zol);CHKERRQ(ierr); 130 131 ierr = DMDACreate(PETSC_COMM_SELF,&da);CHKERRQ(ierr); 132 ierr = DMSetOptionsPrefix(da,"sub_");CHKERRQ(ierr); 133 ierr = DMDASetDim(da, info.dim);CHKERRQ(ierr); 134 ierr = DMDASetDof(da, info.dof);CHKERRQ(ierr); 135 136 ierr = DMDASetStencilType(da,info.st);CHKERRQ(ierr); 137 ierr = DMDASetStencilWidth(da,info.sw);CHKERRQ(ierr); 138 139 /* local with overlap */ 140 xsize = info.xm; 141 ysize = info.ym; 142 zsize = info.zm; 143 xo = info.xs; 144 yo = info.ys; 145 zo = info.zs; 146 if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs != 0)) { 147 xsize += xol; 148 xo -= xol; 149 } 150 if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys != 0)) { 151 ysize += yol; 152 yo -= yol; 153 } 154 if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs != 0)) { 155 zsize += zol; 156 zo -= zol; 157 } 158 159 if (info.bx == DMDA_BOUNDARY_PERIODIC || (info.xs+info.xm != info.mx)) xsize += xol; 160 if (info.by == DMDA_BOUNDARY_PERIODIC || (info.ys+info.ym != info.my)) ysize += yol; 161 if (info.bz == DMDA_BOUNDARY_PERIODIC || (info.zs+info.zm != info.mz)) zsize += zol; 162 163 ierr = DMDASetSizes(da, xsize, ysize, zsize);CHKERRQ(ierr); 164 ierr = DMDASetNumProcs(da, 1, 1, 1);CHKERRQ(ierr); 165 ierr = DMDASetBoundaryType(da, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED);CHKERRQ(ierr); 166 167 /* set up as a block instead */ 168 ierr = DMSetUp(da);CHKERRQ(ierr); 169 170 /* todo - nonuniform coordinates */ 171 ierr = DMDAGetLocalBoundingBox(dm,lmin,lmax);CHKERRQ(ierr); 172 ierr = DMDASetUniformCoordinates(da,lmin[0],lmax[0],lmin[1],lmax[1],lmin[2],lmax[2]);CHKERRQ(ierr); 173 174 /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */ 175 ierr = DMDASetOffset(da,xo,yo,zo,info.mx,info.my,info.mz);CHKERRQ(ierr); 176 177 *dddm = da; 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "DMCreateDomainDecompositionScatters_DA" 183 /* 184 Fills the local vector problem on the subdomain from the global problem. 185 186 Right now this assumes one subdomain per processor. 187 188 */ 189 PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat) 190 { 191 PetscErrorCode ierr; 192 DMDALocalInfo info,subinfo; 193 DM subdm; 194 MatStencil upper,lower; 195 IS idis,isis,odis,osis,gdis; 196 Vec svec,dvec,slvec; 197 198 PetscFunctionBegin; 199 if (nsubdms != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have more than one subdomain per processor (yet)"); 200 201 /* allocate the arrays of scatters */ 202 if (iscat) {ierr = PetscMalloc(sizeof(VecScatter*),iscat);CHKERRQ(ierr);} 203 if (oscat) {ierr = PetscMalloc(sizeof(VecScatter*),oscat);CHKERRQ(ierr);} 204 if (lscat) {ierr = PetscMalloc(sizeof(VecScatter*),lscat);CHKERRQ(ierr);} 205 206 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 207 subdm = subdms[0]; 208 ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr); 209 210 /* create the global and subdomain index sets for the inner domain */ 211 /* TODO - make this actually support multiple subdomains -- subdomain needs to provide where it's nonoverlapping portion belongs */ 212 lower.i = info.xs; 213 lower.j = info.ys; 214 lower.k = info.zs; 215 upper.i = info.xs+info.xm; 216 upper.j = info.ys+info.ym; 217 upper.k = info.zs+info.zm; 218 ierr = DMDACreatePatchIS(dm,&lower,&upper,&idis);CHKERRQ(ierr); 219 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&isis);CHKERRQ(ierr); 220 221 /* create the global and subdomain index sets for the outer subdomain */ 222 lower.i = subinfo.xs; 223 lower.j = subinfo.ys; 224 lower.k = subinfo.zs; 225 upper.i = subinfo.xs+subinfo.xm; 226 upper.j = subinfo.ys+subinfo.ym; 227 upper.k = subinfo.zs+subinfo.zm; 228 ierr = DMDACreatePatchIS(dm,&lower,&upper,&odis);CHKERRQ(ierr); 229 ierr = DMDACreatePatchIS(subdm,&lower,&upper,&osis);CHKERRQ(ierr); 230 231 /* global and subdomain ISes for the local indices of the subdomain */ 232 /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */ 233 lower.i = subinfo.gxs; 234 lower.j = subinfo.gys; 235 lower.k = subinfo.gzs; 236 upper.i = subinfo.gxs+subinfo.gxm; 237 upper.j = subinfo.gys+subinfo.gym; 238 upper.k = subinfo.gzs+subinfo.gzm; 239 240 ierr = DMDACreatePatchIS(dm,&lower,&upper,&gdis);CHKERRQ(ierr); 241 242 /* form the scatter */ 243 ierr = DMGetGlobalVector(dm,&dvec);CHKERRQ(ierr); 244 ierr = DMGetGlobalVector(subdm,&svec);CHKERRQ(ierr); 245 ierr = DMGetLocalVector(subdm,&slvec);CHKERRQ(ierr); 246 247 if (iscat) {ierr = VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[0]);CHKERRQ(ierr);} 248 if (oscat) {ierr = VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[0]);CHKERRQ(ierr);} 249 if (lscat) {ierr = VecScatterCreate(dvec,gdis,slvec,NULL,&(*lscat)[0]);CHKERRQ(ierr);} 250 251 ierr = DMRestoreGlobalVector(dm,&dvec);CHKERRQ(ierr); 252 ierr = DMRestoreGlobalVector(subdm,&svec);CHKERRQ(ierr); 253 ierr = DMRestoreLocalVector(subdm,&slvec);CHKERRQ(ierr); 254 255 ierr = ISDestroy(&idis);CHKERRQ(ierr); 256 ierr = ISDestroy(&isis);CHKERRQ(ierr); 257 258 ierr = ISDestroy(&odis);CHKERRQ(ierr); 259 ierr = ISDestroy(&osis);CHKERRQ(ierr); 260 261 ierr = ISDestroy(&gdis);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 #undef __FUNCT__ 266 #define __FUNCT__ "DMDASubDomainIS_Private" 267 /* We have that the interior regions are going to be the same, but the ghost regions might not match up 268 269 ---------- 270 ---------- 271 --++++++o= 272 --++++++o= 273 --++++++o= 274 --++++++o= 275 --ooooooo= 276 --======== 277 278 Therefore, for each point in the overall, we must check if it's: 279 280 1. +: In the interior of the global dm; it lines up 281 2. o: In the overlap region -- for now the same as 1; no overlap 282 3. =: In the shared ghost region -- handled by DMCreateDomainDecompositionLocalScatter() 283 4. -: In the nonshared ghost region 284 */ 285 286 PetscErrorCode DMDASubDomainIS_Private(DM dm,DM subdm,IS *iis,IS *ois) 287 { 288 PetscErrorCode ierr; 289 DMDALocalInfo info,subinfo; 290 MatStencil lower,upper; 291 292 PetscFunctionBegin; 293 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 294 ierr = DMDAGetLocalInfo(subdm,&subinfo);CHKERRQ(ierr); 295 296 /* create the inner IS */ 297 lower.i = info.xs; 298 lower.j = info.ys; 299 lower.k = info.zs; 300 upper.i = info.xs+info.xm; 301 upper.j = info.ys+info.ym; 302 upper.k = info.zs+info.zm; 303 304 ierr = DMDACreatePatchIS(dm,&lower,&upper,iis);CHKERRQ(ierr); 305 306 /* create the outer IS */ 307 lower.i = subinfo.xs; 308 lower.j = subinfo.ys; 309 lower.k = subinfo.zs; 310 upper.i = subinfo.xs+subinfo.xm; 311 upper.j = subinfo.ys+subinfo.ym; 312 upper.k = subinfo.zs+subinfo.zm; 313 ierr = DMDACreatePatchIS(dm,&lower,&upper,ois);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 #undef __FUNCT__ 318 #define __FUNCT__ "DMCreateDomainDecomposition_DA" 319 PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm) 320 { 321 PetscErrorCode ierr; 322 IS iis0,ois0; 323 DM subdm0; 324 325 PetscFunctionBegin; 326 if (len) *len = 1; 327 328 if (iis) {ierr = PetscMalloc(sizeof(IS*),iis);CHKERRQ(ierr);} 329 if (ois) {ierr = PetscMalloc(sizeof(IS*),ois);CHKERRQ(ierr);} 330 if (subdm) {ierr = PetscMalloc(sizeof(DM*),subdm);CHKERRQ(ierr);} 331 if (names) {ierr = PetscMalloc(sizeof(char*),names);CHKERRQ(ierr);} 332 ierr = DMDASubDomainDA_Private(dm,&subdm0);CHKERRQ(ierr); 333 ierr = DMDASubDomainIS_Private(dm,subdm0,&iis0,&ois0);CHKERRQ(ierr); 334 if (iis) (*iis)[0] = iis0; 335 else { 336 ierr = ISDestroy(&iis0);CHKERRQ(ierr); 337 } 338 if (ois) (*ois)[0] = ois0; 339 else { 340 ierr = ISDestroy(&ois0);CHKERRQ(ierr); 341 } 342 if (subdm) (*subdm)[0] = subdm0; 343 if (names) (*names)[0] = 0; 344 PetscFunctionReturn(0); 345 } 346