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