1 2 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMView_DA_2d" 6 PetscErrorCode DMView_DA_2d(DM da,PetscViewer viewer) 7 { 8 PetscErrorCode ierr; 9 PetscMPIInt rank; 10 PetscBool iascii,isdraw,isbinary; 11 DM_DA *dd = (DM_DA*)da->data; 12 #if defined(PETSC_HAVE_MATLAB_ENGINE) 13 PetscBool ismatlab; 14 #endif 15 16 PetscFunctionBegin; 17 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 18 19 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 20 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 21 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 22 #if defined(PETSC_HAVE_MATLAB_ENGINE) 23 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);CHKERRQ(ierr); 24 #endif 25 if (iascii) { 26 PetscViewerFormat format; 27 28 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 29 if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) { 30 DMDALocalInfo info; 31 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 32 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 33 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D m %D n %D w %D s %D\n",rank,dd->M,dd->N,dd->m,dd->n,dd->w,dd->s);CHKERRQ(ierr); 34 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\n",info.xs,info.xs+info.xm,info.ys,info.ys+info.ym);CHKERRQ(ierr); 35 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 36 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 37 } else { 38 ierr = DMView_DA_VTK(da,viewer);CHKERRQ(ierr); 39 } 40 } else if (isdraw) { 41 PetscDraw draw; 42 double ymin = -1*dd->s-1,ymax = dd->N+dd->s; 43 double xmin = -1*dd->s-1,xmax = dd->M+dd->s; 44 double x,y; 45 PetscInt base,*idx; 46 char node[10]; 47 PetscBool isnull; 48 49 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 50 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 51 if (!da->coordinates) { 52 ierr = PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);CHKERRQ(ierr); 53 } 54 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 55 56 /* first processor draw all node lines */ 57 if (!rank) { 58 ymin = 0.0; ymax = dd->N - 1; 59 for (xmin=0; xmin<dd->M; xmin++) { 60 ierr = PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);CHKERRQ(ierr); 61 } 62 xmin = 0.0; xmax = dd->M - 1; 63 for (ymin=0; ymin<dd->N; ymin++) { 64 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);CHKERRQ(ierr); 65 } 66 } 67 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 68 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 69 70 /* draw my box */ 71 ymin = dd->ys; ymax = dd->ye - 1; xmin = dd->xs/dd->w; 72 xmax =(dd->xe-1)/dd->w; 73 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);CHKERRQ(ierr); 74 ierr = PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 75 ierr = PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 76 ierr = PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 77 78 /* put in numbers */ 79 base = (dd->base)/dd->w; 80 for (y=ymin; y<=ymax; y++) { 81 for (x=xmin; x<=xmax; x++) { 82 sprintf(node,"%d",(int)base++); 83 ierr = PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);CHKERRQ(ierr); 84 } 85 } 86 87 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 88 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 89 /* overlay ghost numbers, useful for error checking */ 90 /* put in numbers */ 91 92 base = 0; idx = dd->idx; 93 ymin = dd->Ys; ymax = dd->Ye; xmin = dd->Xs; xmax = dd->Xe; 94 for (y=ymin; y<ymax; y++) { 95 for (x=xmin; x<xmax; x++) { 96 if ((base % dd->w) == 0) { 97 sprintf(node,"%d",(int)(idx[base]/dd->w)); 98 ierr = PetscDrawString(draw,x/dd->w,y,PETSC_DRAW_BLUE,node);CHKERRQ(ierr); 99 } 100 base++; 101 } 102 } 103 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 104 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 105 } else if (isbinary) { 106 ierr = DMView_DA_Binary(da,viewer);CHKERRQ(ierr); 107 #if defined(PETSC_HAVE_MATLAB_ENGINE) 108 } else if (ismatlab) { 109 ierr = DMView_DA_Matlab(da,viewer);CHKERRQ(ierr); 110 #endif 111 } 112 PetscFunctionReturn(0); 113 } 114 115 /* 116 M is number of grid points 117 m is number of processors 118 119 */ 120 #undef __FUNCT__ 121 #define __FUNCT__ "DMDASplitComm2d" 122 PetscErrorCode DMDASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm) 123 { 124 PetscErrorCode ierr; 125 PetscInt m,n = 0,x = 0,y = 0; 126 PetscMPIInt size,csize,rank; 127 128 PetscFunctionBegin; 129 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 130 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 131 132 csize = 4*size; 133 do { 134 if (csize % 4) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot split communicator of size %d tried %d %D %D",size,csize,x,y); 135 csize = csize/4; 136 137 m = (PetscInt)(0.5 + sqrt(((double)M)*((double)csize)/((double)N))); 138 if (!m) m = 1; 139 while (m > 0) { 140 n = csize/m; 141 if (m*n == csize) break; 142 m--; 143 } 144 if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;} 145 146 x = M/m + ((M % m) > ((csize-1) % m)); 147 y = (N + (csize-1)/m)/n; 148 } while ((x < 4 || y < 4) && csize > 1); 149 if (size != csize) { 150 MPI_Group entire_group,sub_group; 151 PetscMPIInt i,*groupies; 152 153 ierr = MPI_Comm_group(comm,&entire_group);CHKERRQ(ierr); 154 ierr = PetscMalloc(csize*sizeof(PetscInt),&groupies);CHKERRQ(ierr); 155 for (i=0; i<csize; i++) { 156 groupies[i] = (rank/csize)*csize + i; 157 } 158 ierr = MPI_Group_incl(entire_group,csize,groupies,&sub_group);CHKERRQ(ierr); 159 ierr = PetscFree(groupies);CHKERRQ(ierr); 160 ierr = MPI_Comm_create(comm,sub_group,outcomm);CHKERRQ(ierr); 161 ierr = MPI_Group_free(&entire_group);CHKERRQ(ierr); 162 ierr = MPI_Group_free(&sub_group);CHKERRQ(ierr); 163 ierr = PetscInfo1(0,"DMDASplitComm2d:Creating redundant coarse problems of size %d\n",csize);CHKERRQ(ierr); 164 } else { 165 *outcomm = comm; 166 } 167 PetscFunctionReturn(0); 168 } 169 170 #if defined(new) 171 #undef __FUNCT__ 172 #define __FUNCT__ "DMDAGetDiagonal_MFFD" 173 /* 174 DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local 175 function lives on a DMDA 176 177 y ~= (F(u + ha) - F(u))/h, 178 where F = nonlinear function, as set by SNESSetFunction() 179 u = current iterate 180 h = difference interval 181 */ 182 PetscErrorCode DMDAGetDiagonal_MFFD(DM da,Vec U,Vec a) 183 { 184 PetscScalar h,*aa,*ww,v; 185 PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON; 186 PetscErrorCode ierr; 187 PetscInt gI,nI; 188 MatStencil stencil; 189 DMDALocalInfo info; 190 191 PetscFunctionBegin; 192 ierr = (*ctx->func)(0,U,a,ctx->funcctx);CHKERRQ(ierr); 193 ierr = (*ctx->funcisetbase)(U,ctx->funcctx);CHKERRQ(ierr); 194 195 ierr = VecGetArray(U,&ww);CHKERRQ(ierr); 196 ierr = VecGetArray(a,&aa);CHKERRQ(ierr); 197 198 nI = 0; 199 h = ww[gI]; 200 if (h == 0.0) h = 1.0; 201 if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin; 202 else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin; 203 h *= epsilon; 204 205 ww[gI] += h; 206 ierr = (*ctx->funci)(i,w,&v,ctx->funcctx);CHKERRQ(ierr); 207 aa[nI] = (v - aa[nI])/h; 208 ww[gI] -= h; 209 nI++; 210 } 211 ierr = VecRestoreArray(U,&ww);CHKERRQ(ierr); 212 ierr = VecRestoreArray(a,&aa);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 #endif 216 217 #undef __FUNCT__ 218 #define __FUNCT__ "DMSetUp_DA_2D" 219 PetscErrorCode DMSetUp_DA_2D(DM da) 220 { 221 DM_DA *dd = (DM_DA*)da->data; 222 const PetscInt M = dd->M; 223 const PetscInt N = dd->N; 224 PetscInt m = dd->m; 225 PetscInt n = dd->n; 226 const PetscInt dof = dd->w; 227 const PetscInt s = dd->s; 228 DMDABoundaryType bx = dd->bx; 229 DMDABoundaryType by = dd->by; 230 DMDAStencilType stencil_type = dd->stencil_type; 231 PetscInt *lx = dd->lx; 232 PetscInt *ly = dd->ly; 233 MPI_Comm comm; 234 PetscMPIInt rank,size; 235 PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end,IXs,IXe,IYs,IYe; 236 PetscInt up,down,left,right,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn,*idx_cpy; 237 const PetscInt *idx_full; 238 PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count; 239 PetscInt s_x,s_y; /* s proportionalized to w */ 240 PetscInt sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0; 241 Vec local,global; 242 VecScatter ltog,gtol; 243 IS to,from,ltogis; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 if (stencil_type == DMDA_STENCIL_BOX && (bx == DMDA_BOUNDARY_MIRROR || by == DMDA_BOUNDARY_MIRROR)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Mirror boundary and box stencil"); 248 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 249 #if !defined(PETSC_USE_64BIT_INDICES) 250 if (((Petsc64bitInt) M)*((Petsc64bitInt) N)*((Petsc64bitInt) dof) > (Petsc64bitInt) PETSC_MPI_INT_MAX) SETERRQ3(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D (dof) is too large for 32 bit indices",M,N,dof); 251 #endif 252 253 if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof); 254 if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s); 255 256 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 257 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 258 259 if (m != PETSC_DECIDE) { 260 if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m); 261 else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size); 262 } 263 if (n != PETSC_DECIDE) { 264 if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n); 265 else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size); 266 } 267 268 if (m == PETSC_DECIDE || n == PETSC_DECIDE) { 269 if (n != PETSC_DECIDE) { 270 m = size/n; 271 } else if (m != PETSC_DECIDE) { 272 n = size/m; 273 } else { 274 /* try for squarish distribution */ 275 m = (PetscInt)(0.5 + sqrt(((double)M)*((double)size)/((double)N))); 276 if (!m) m = 1; 277 while (m > 0) { 278 n = size/m; 279 if (m*n == size) break; 280 m--; 281 } 282 if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;} 283 } 284 if (m*n != size) SETERRQ(comm,PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n "); 285 } else if (m*n != size) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition"); 286 287 if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m); 288 if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n); 289 290 /* 291 Determine locally owned region 292 xs is the first local node number, x is the number of local nodes 293 */ 294 if (!lx) { 295 ierr = PetscMalloc(m*sizeof(PetscInt), &dd->lx);CHKERRQ(ierr); 296 lx = dd->lx; 297 for (i=0; i<m; i++) { 298 lx[i] = M/m + ((M % m) > i); 299 } 300 } 301 x = lx[rank % m]; 302 xs = 0; 303 for (i=0; i<(rank % m); i++) { 304 xs += lx[i]; 305 } 306 #if defined(PETSC_USE_DEBUG) 307 left = xs; 308 for (i=(rank % m); i<m; i++) { 309 left += lx[i]; 310 } 311 if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M); 312 #endif 313 314 /* 315 Determine locally owned region 316 ys is the first local node number, y is the number of local nodes 317 */ 318 if (!ly) { 319 ierr = PetscMalloc(n*sizeof(PetscInt), &dd->ly);CHKERRQ(ierr); 320 ly = dd->ly; 321 for (i=0; i<n; i++) { 322 ly[i] = N/n + ((N % n) > i); 323 } 324 } 325 y = ly[rank/m]; 326 ys = 0; 327 for (i=0; i<(rank/m); i++) { 328 ys += ly[i]; 329 } 330 #if defined(PETSC_USE_DEBUG) 331 left = ys; 332 for (i=(rank/m); i<n; i++) { 333 left += ly[i]; 334 } 335 if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N); 336 #endif 337 338 /* 339 check if the scatter requires more than one process neighbor or wraps around 340 the domain more than once 341 */ 342 if ((x < s) && ((m > 1) || (bx == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s); 343 if ((y < s) && ((n > 1) || (by == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s); 344 xe = xs + x; 345 ye = ys + y; 346 347 /* determine ghost region (Xs) and region scattered into (IXs) */ 348 if (xs-s > 0) { 349 Xs = xs - s; IXs = xs - s; 350 } else { 351 if (bx) { 352 Xs = xs - s; 353 } else { 354 Xs = 0; 355 } 356 IXs = 0; 357 } 358 if (xe+s <= M) { 359 Xe = xe + s; IXe = xe + s; 360 } else { 361 if (bx) { 362 Xs = xs - s; Xe = xe + s; 363 } else { 364 Xe = M; 365 } 366 IXe = M; 367 } 368 369 if (bx == DMDA_BOUNDARY_PERIODIC || bx == DMDA_BOUNDARY_MIRROR) { 370 IXs = xs - s; 371 IXe = xe + s; 372 Xs = xs - s; 373 Xe = xe + s; 374 } 375 376 if (ys-s > 0) { 377 Ys = ys - s; IYs = ys - s; 378 } else { 379 if (by) { 380 Ys = ys - s; 381 } else { 382 Ys = 0; 383 } 384 IYs = 0; 385 } 386 if (ye+s <= N) { 387 Ye = ye + s; IYe = ye + s; 388 } else { 389 if (by) { 390 Ye = ye + s; 391 } else { 392 Ye = N; 393 } 394 IYe = N; 395 } 396 397 if (by == DMDA_BOUNDARY_PERIODIC || by == DMDA_BOUNDARY_MIRROR) { 398 IYs = ys - s; 399 IYe = ye + s; 400 Ys = ys - s; 401 Ye = ye + s; 402 } 403 404 /* stencil length in each direction */ 405 s_x = s; 406 s_y = s; 407 408 /* determine starting point of each processor */ 409 nn = x*y; 410 ierr = PetscMalloc2(size+1,PetscInt,&bases,size,PetscInt,&ldims);CHKERRQ(ierr); 411 ierr = MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);CHKERRQ(ierr); 412 bases[0] = 0; 413 for (i=1; i<=size; i++) { 414 bases[i] = ldims[i-1]; 415 } 416 for (i=1; i<=size; i++) { 417 bases[i] += bases[i-1]; 418 } 419 base = bases[rank]*dof; 420 421 /* allocate the base parallel and sequential vectors */ 422 dd->Nlocal = x*y*dof; 423 ierr = VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);CHKERRQ(ierr); 424 dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof; 425 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);CHKERRQ(ierr); 426 427 /* generate appropriate vector scatters */ 428 /* local to global inserts non-ghost point region into global */ 429 ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr); 430 ierr = ISCreateStride(comm,x*y*dof,start,1,&to);CHKERRQ(ierr); 431 432 ierr = PetscMalloc(x*y*sizeof(PetscInt),&idx);CHKERRQ(ierr); 433 left = xs - Xs; right = left + x; 434 down = ys - Ys; up = down + y; 435 count = 0; 436 for (i=down; i<up; i++) { 437 for (j=left; j<right; j++) { 438 idx[count++] = i*(Xe-Xs) + j; 439 } 440 } 441 442 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr); 443 ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr); 444 ierr = PetscLogObjectParent(dd,ltog);CHKERRQ(ierr); 445 ierr = ISDestroy(&from);CHKERRQ(ierr); 446 ierr = ISDestroy(&to);CHKERRQ(ierr); 447 448 /* global to local must include ghost points within the domain, 449 but not ghost points outside the domain that aren't periodic */ 450 if (stencil_type == DMDA_STENCIL_BOX) { 451 count = (IXe-IXs)*(IYe-IYs); 452 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 453 454 left = IXs - Xs; right = left + (IXe-IXs); 455 down = IYs - Ys; up = down + (IYe-IYs); 456 count = 0; 457 for (i=down; i<up; i++) { 458 for (j=left; j<right; j++) { 459 idx[count++] = j + i*(Xe-Xs); 460 } 461 } 462 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 463 464 } else { 465 /* must drop into cross shape region */ 466 /* ---------| 467 | top | 468 |--- ---| up 469 | middle | 470 | | 471 ---- ---- down 472 | bottom | 473 ----------- 474 Xs xs xe Xe */ 475 count = (ys-IYs)*x + y*(IXe-IXs) + (IYe-ye)*x; 476 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 477 478 left = xs - Xs; right = left + x; 479 down = ys - Ys; up = down + y; 480 count = 0; 481 /* bottom */ 482 for (i=(IYs-Ys); i<down; i++) { 483 for (j=left; j<right; j++) { 484 idx[count++] = j + i*(Xe-Xs); 485 } 486 } 487 /* middle */ 488 for (i=down; i<up; i++) { 489 for (j=(IXs-Xs); j<(IXe-Xs); j++) { 490 idx[count++] = j + i*(Xe-Xs); 491 } 492 } 493 /* top */ 494 for (i=up; i<up+IYe-ye; i++) { 495 for (j=left; j<right; j++) { 496 idx[count++] = j + i*(Xe-Xs); 497 } 498 } 499 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 500 } 501 502 503 /* determine who lies on each side of us stored in n6 n7 n8 504 n3 n5 505 n0 n1 n2 506 */ 507 508 /* Assume the Non-Periodic Case */ 509 n1 = rank - m; 510 if (rank % m) { 511 n0 = n1 - 1; 512 } else { 513 n0 = -1; 514 } 515 if ((rank+1) % m) { 516 n2 = n1 + 1; 517 n5 = rank + 1; 518 n8 = rank + m + 1; if (n8 >= m*n) n8 = -1; 519 } else { 520 n2 = -1; n5 = -1; n8 = -1; 521 } 522 if (rank % m) { 523 n3 = rank - 1; 524 n6 = n3 + m; if (n6 >= m*n) n6 = -1; 525 } else { 526 n3 = -1; n6 = -1; 527 } 528 n7 = rank + m; if (n7 >= m*n) n7 = -1; 529 530 if (bx == DMDA_BOUNDARY_PERIODIC && by == DMDA_BOUNDARY_PERIODIC) { 531 /* Modify for Periodic Cases */ 532 /* Handle all four corners */ 533 if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1; 534 if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0; 535 if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m; 536 if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1; 537 538 /* Handle Top and Bottom Sides */ 539 if (n1 < 0) n1 = rank + m * (n-1); 540 if (n7 < 0) n7 = rank - m * (n-1); 541 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 542 if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1; 543 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 544 if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1; 545 546 /* Handle Left and Right Sides */ 547 if (n3 < 0) n3 = rank + (m-1); 548 if (n5 < 0) n5 = rank - (m-1); 549 if ((n1 >= 0) && (n0 < 0)) n0 = rank-1; 550 if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1; 551 if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1; 552 if ((n7 >= 0) && (n8 < 0)) n8 = rank+1; 553 } else if (by == DMDA_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */ 554 if (n1 < 0) n1 = rank + m * (n-1); 555 if (n7 < 0) n7 = rank - m * (n-1); 556 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 557 if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1; 558 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 559 if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1; 560 } else if (bx == DMDA_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */ 561 if (n3 < 0) n3 = rank + (m-1); 562 if (n5 < 0) n5 = rank - (m-1); 563 if ((n1 >= 0) && (n0 < 0)) n0 = rank-1; 564 if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1; 565 if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1; 566 if ((n7 >= 0) && (n8 < 0)) n8 = rank+1; 567 } 568 569 ierr = PetscMalloc(9*sizeof(PetscInt),&dd->neighbors);CHKERRQ(ierr); 570 dd->neighbors[0] = n0; 571 dd->neighbors[1] = n1; 572 dd->neighbors[2] = n2; 573 dd->neighbors[3] = n3; 574 dd->neighbors[4] = rank; 575 dd->neighbors[5] = n5; 576 dd->neighbors[6] = n6; 577 dd->neighbors[7] = n7; 578 dd->neighbors[8] = n8; 579 580 if (stencil_type == DMDA_STENCIL_STAR) { 581 /* save corner processor numbers */ 582 sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8; 583 n0 = n2 = n6 = n8 = -1; 584 } 585 586 ierr = PetscMalloc((Xe-Xs)*(Ye-Ys)*sizeof(PetscInt),&idx);CHKERRQ(ierr); 587 ierr = PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*sizeof(PetscInt));CHKERRQ(ierr); 588 589 nn = 0; 590 xbase = bases[rank]; 591 for (i=1; i<=s_y; i++) { 592 if (n0 >= 0) { /* left below */ 593 x_t = lx[n0 % m]; 594 y_t = ly[(n0/m)]; 595 s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x; 596 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 597 } 598 599 if (n1 >= 0) { /* directly below */ 600 x_t = x; 601 y_t = ly[(n1/m)]; 602 s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t; 603 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 604 } else if (by == DMDA_BOUNDARY_MIRROR) { 605 for (j=0; j<x; j++) { idx[nn++] = bases[rank] + x*(s_y - i + 1) + j;} 606 } 607 608 if (n2 >= 0) { /* right below */ 609 x_t = lx[n2 % m]; 610 y_t = ly[(n2/m)]; 611 s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t; 612 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 613 } 614 } 615 616 for (i=0; i<y; i++) { 617 if (n3 >= 0) { /* directly left */ 618 x_t = lx[n3 % m]; 619 /* y_t = y; */ 620 s_t = bases[n3] + (i+1)*x_t - s_x; 621 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 622 } else if (bx == DMDA_BOUNDARY_MIRROR) { 623 for (j=0; j<s_x; j++) { idx[nn++] = bases[rank] + x*i + s_x - j;} 624 } 625 626 for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */ 627 628 if (n5 >= 0) { /* directly right */ 629 x_t = lx[n5 % m]; 630 /* y_t = y; */ 631 s_t = bases[n5] + (i)*x_t; 632 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 633 } else if (bx == DMDA_BOUNDARY_MIRROR) { 634 for (j=0; j<s_x; j++) { idx[nn++] = bases[rank] + x*(i + 1) - 2 - j;} 635 } 636 } 637 638 for (i=1; i<=s_y; i++) { 639 if (n6 >= 0) { /* left above */ 640 x_t = lx[n6 % m]; 641 /* y_t = ly[(n6/m)]; */ 642 s_t = bases[n6] + (i)*x_t - s_x; 643 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 644 } 645 646 if (n7 >= 0) { /* directly above */ 647 x_t = x; 648 /* y_t = ly[(n7/m)]; */ 649 s_t = bases[n7] + (i-1)*x_t; 650 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 651 } else if (by == DMDA_BOUNDARY_MIRROR) { 652 for (j=0; j<x; j++) { idx[nn++] = bases[rank] + x*(y - i - 1) + j;} 653 } 654 655 if (n8 >= 0) { /* right above */ 656 x_t = lx[n8 % m]; 657 /* y_t = ly[(n8/m)]; */ 658 s_t = bases[n8] + (i-1)*x_t; 659 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 660 } 661 } 662 663 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 664 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 665 ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr); 666 ierr = ISDestroy(&to);CHKERRQ(ierr); 667 ierr = ISDestroy(&from);CHKERRQ(ierr); 668 669 if (stencil_type == DMDA_STENCIL_STAR) { 670 n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8; 671 } 672 673 if (((stencil_type == DMDA_STENCIL_STAR) || 674 (bx && bx != DMDA_BOUNDARY_PERIODIC) || 675 (by && by != DMDA_BOUNDARY_PERIODIC))) { 676 /* 677 Recompute the local to global mappings, this time keeping the 678 information about the cross corner processor numbers and any ghosted 679 but not periodic indices. 680 */ 681 nn = 0; 682 xbase = bases[rank]; 683 for (i=1; i<=s_y; i++) { 684 if (n0 >= 0) { /* left below */ 685 x_t = lx[n0 % m]; 686 y_t = ly[(n0/m)]; 687 s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x; 688 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 689 } else if (xs-Xs > 0 && ys-Ys > 0) { 690 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 691 } 692 if (n1 >= 0) { /* directly below */ 693 x_t = x; 694 y_t = ly[(n1/m)]; 695 s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t; 696 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 697 } else if (ys-Ys > 0) { 698 if (by == DMDA_BOUNDARY_MIRROR) { 699 for (j=0; j<x; j++) { idx[nn++] = bases[rank] + x*(s_y - i + 1) + j;} 700 } else { 701 for (j=0; j<x; j++) { idx[nn++] = -1;} 702 } 703 } 704 if (n2 >= 0) { /* right below */ 705 x_t = lx[n2 % m]; 706 y_t = ly[(n2/m)]; 707 s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t; 708 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 709 } else if (Xe-xe> 0 && ys-Ys > 0) { 710 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 711 } 712 } 713 714 for (i=0; i<y; i++) { 715 if (n3 >= 0) { /* directly left */ 716 x_t = lx[n3 % m]; 717 /* y_t = y; */ 718 s_t = bases[n3] + (i+1)*x_t - s_x; 719 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 720 } else if (xs-Xs > 0) { 721 if (bx == DMDA_BOUNDARY_MIRROR) { 722 for (j=0; j<s_x; j++) { idx[nn++] = bases[rank] + x*i + s_x - j;} 723 } else { 724 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 725 } 726 } 727 728 for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */ 729 730 if (n5 >= 0) { /* directly right */ 731 x_t = lx[n5 % m]; 732 /* y_t = y; */ 733 s_t = bases[n5] + (i)*x_t; 734 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 735 } else if (Xe-xe > 0) { 736 if (bx == DMDA_BOUNDARY_MIRROR) { 737 for (j=0; j<s_x; j++) { idx[nn++] = bases[rank] + x*(i + 1) - 2 - j;} 738 } else { 739 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 740 } 741 } 742 } 743 744 for (i=1; i<=s_y; i++) { 745 if (n6 >= 0) { /* left above */ 746 x_t = lx[n6 % m]; 747 /* y_t = ly[(n6/m)]; */ 748 s_t = bases[n6] + (i)*x_t - s_x; 749 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 750 } else if (xs-Xs > 0 && Ye-ye > 0) { 751 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 752 } 753 if (n7 >= 0) { /* directly above */ 754 x_t = x; 755 /* y_t = ly[(n7/m)]; */ 756 s_t = bases[n7] + (i-1)*x_t; 757 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 758 } else if (Ye-ye > 0) { 759 if (by == DMDA_BOUNDARY_MIRROR) { 760 for (j=0; j<x; j++) { idx[nn++] = bases[rank] + x*(y - i - 1) + j;} 761 } else { 762 for (j=0; j<x; j++) { idx[nn++] = -1;} 763 } 764 } 765 if (n8 >= 0) { /* right above */ 766 x_t = lx[n8 % m]; 767 /* y_t = ly[(n8/m)]; */ 768 s_t = bases[n8] + (i-1)*x_t; 769 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 770 } else if (Xe-xe > 0 && Ye-ye > 0) { 771 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 772 } 773 } 774 } 775 /* 776 Set the local to global ordering in the global vector, this allows use 777 of VecSetValuesLocal(). 778 */ 779 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);CHKERRQ(ierr); 780 ierr = PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);CHKERRQ(ierr); 781 ierr = PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 782 ierr = ISGetIndices(ltogis, &idx_full);CHKERRQ(ierr); 783 ierr = PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 784 ierr = ISRestoreIndices(ltogis, &idx_full);CHKERRQ(ierr); 785 ierr = ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);CHKERRQ(ierr); 786 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 787 ierr = ISDestroy(<ogis);CHKERRQ(ierr); 788 ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr); 789 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 790 791 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 792 dd->m = m; dd->n = n; 793 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 794 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1; 795 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1; 796 797 ierr = VecDestroy(&local);CHKERRQ(ierr); 798 ierr = VecDestroy(&global);CHKERRQ(ierr); 799 800 dd->gtol = gtol; 801 dd->ltog = ltog; 802 dd->idx = idx_cpy; 803 dd->Nl = nn*dof; 804 dd->base = base; 805 da->ops->view = DMView_DA_2d; 806 dd->ltol = PETSC_NULL; 807 dd->ao = PETSC_NULL; 808 PetscFunctionReturn(0); 809 } 810 811 #undef __FUNCT__ 812 #define __FUNCT__ "DMDACreate2d" 813 /*@C 814 DMDACreate2d - Creates an object that will manage the communication of two-dimensional 815 regular array data that is distributed across some processors. 816 817 Collective on MPI_Comm 818 819 Input Parameters: 820 + comm - MPI communicator 821 . bx,by - type of ghost nodes the array have. 822 Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC. 823 . stencil_type - stencil type. Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR. 824 . M,N - global dimension in each direction of the array (use -M and or -N to indicate that it may be set to a different value 825 from the command line with -da_grid_x <M> -da_grid_y <N>) 826 . m,n - corresponding number of processors in each dimension 827 (or PETSC_DECIDE to have calculated) 828 . dof - number of degrees of freedom per node 829 . s - stencil width 830 - lx, ly - arrays containing the number of nodes in each cell along 831 the x and y coordinates, or PETSC_NULL. If non-null, these 832 must be of length as m and n, and the corresponding 833 m and n cannot be PETSC_DECIDE. The sum of the lx[] entries 834 must be M, and the sum of the ly[] entries must be N. 835 836 Output Parameter: 837 . da - the resulting distributed array object 838 839 Options Database Key: 840 + -dm_view - Calls DMView() at the conclusion of DMDACreate2d() 841 . -da_grid_x <nx> - number of grid points in x direction, if M < 0 842 . -da_grid_y <ny> - number of grid points in y direction, if N < 0 843 . -da_processors_x <nx> - number of processors in x direction 844 . -da_processors_y <ny> - number of processors in y direction 845 . -da_refine_x <rx> - refinement ratio in x direction 846 . -da_refine_y <ry> - refinement ratio in y direction 847 - -da_refine <n> - refine the DMDA n times before creating, if M or N < 0 848 849 850 Level: beginner 851 852 Notes: 853 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 854 standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 855 the standard 9-pt stencil. 856 857 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 858 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 859 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 860 861 .keywords: distributed array, create, two-dimensional 862 863 .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(), 864 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(), 865 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges() 866 867 @*/ 868 869 PetscErrorCode DMDACreate2d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDAStencilType stencil_type, 870 PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da) 871 { 872 PetscErrorCode ierr; 873 874 PetscFunctionBegin; 875 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 876 ierr = DMDASetDim(*da, 2);CHKERRQ(ierr); 877 ierr = DMDASetSizes(*da, M, N, 1);CHKERRQ(ierr); 878 ierr = DMDASetNumProcs(*da, m, n, PETSC_DECIDE);CHKERRQ(ierr); 879 ierr = DMDASetBoundaryType(*da, bx, by, DMDA_BOUNDARY_NONE);CHKERRQ(ierr); 880 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 881 ierr = DMDASetStencilType(*da, stencil_type);CHKERRQ(ierr); 882 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 883 ierr = DMDASetOwnershipRanges(*da, lx, ly, PETSC_NULL);CHKERRQ(ierr); 884 /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */ 885 ierr = DMSetFromOptions(*da);CHKERRQ(ierr); 886 ierr = DMSetUp(*da);CHKERRQ(ierr); 887 ierr = DMViewFromOptions(*da,"-dm_view");CHKERRQ(ierr); 888 PetscFunctionReturn(0); 889 } 890