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 if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof); 224 if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s); 225 226 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 227 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 228 229 dd->p = 1; 230 if (m != PETSC_DECIDE) { 231 if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m); 232 else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size); 233 } 234 if (n != PETSC_DECIDE) { 235 if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n); 236 else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size); 237 } 238 239 if (m == PETSC_DECIDE || n == PETSC_DECIDE) { 240 if (n != PETSC_DECIDE) { 241 m = size/n; 242 } else if (m != PETSC_DECIDE) { 243 n = size/m; 244 } else { 245 /* try for squarish distribution */ 246 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N))); 247 if (!m) m = 1; 248 while (m > 0) { 249 n = size/m; 250 if (m*n == size) break; 251 m--; 252 } 253 if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;} 254 } 255 if (m*n != size) SETERRQ(comm,PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n "); 256 } else if (m*n != size) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition"); 257 258 if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m); 259 if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n); 260 261 /* 262 Determine locally owned region 263 xs is the first local node number, x is the number of local nodes 264 */ 265 if (!lx) { 266 ierr = PetscMalloc1(m, &dd->lx);CHKERRQ(ierr); 267 lx = dd->lx; 268 for (i=0; i<m; i++) { 269 lx[i] = M/m + ((M % m) > i); 270 } 271 } 272 x = lx[rank % m]; 273 xs = 0; 274 for (i=0; i<(rank % m); i++) { 275 xs += lx[i]; 276 } 277 #if defined(PETSC_USE_DEBUG) 278 left = xs; 279 for (i=(rank % m); i<m; i++) { 280 left += lx[i]; 281 } 282 if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M); 283 #endif 284 285 /* 286 Determine locally owned region 287 ys is the first local node number, y is the number of local nodes 288 */ 289 if (!ly) { 290 ierr = PetscMalloc1(n, &dd->ly);CHKERRQ(ierr); 291 ly = dd->ly; 292 for (i=0; i<n; i++) { 293 ly[i] = N/n + ((N % n) > i); 294 } 295 } 296 y = ly[rank/m]; 297 ys = 0; 298 for (i=0; i<(rank/m); i++) { 299 ys += ly[i]; 300 } 301 #if defined(PETSC_USE_DEBUG) 302 left = ys; 303 for (i=(rank/m); i<n; i++) { 304 left += ly[i]; 305 } 306 if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N); 307 #endif 308 309 /* 310 check if the scatter requires more than one process neighbor or wraps around 311 the domain more than once 312 */ 313 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); 314 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); 315 xe = xs + x; 316 ye = ys + y; 317 318 /* determine ghost region (Xs) and region scattered into (IXs) */ 319 if (xs-s > 0) { 320 Xs = xs - s; IXs = xs - s; 321 } else { 322 if (bx) { 323 Xs = xs - s; 324 } else { 325 Xs = 0; 326 } 327 IXs = 0; 328 } 329 if (xe+s <= M) { 330 Xe = xe + s; IXe = xe + s; 331 } else { 332 if (bx) { 333 Xs = xs - s; Xe = xe + s; 334 } else { 335 Xe = M; 336 } 337 IXe = M; 338 } 339 340 if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) { 341 IXs = xs - s; 342 IXe = xe + s; 343 Xs = xs - s; 344 Xe = xe + s; 345 } 346 347 if (ys-s > 0) { 348 Ys = ys - s; IYs = ys - s; 349 } else { 350 if (by) { 351 Ys = ys - s; 352 } else { 353 Ys = 0; 354 } 355 IYs = 0; 356 } 357 if (ye+s <= N) { 358 Ye = ye + s; IYe = ye + s; 359 } else { 360 if (by) { 361 Ye = ye + s; 362 } else { 363 Ye = N; 364 } 365 IYe = N; 366 } 367 368 if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) { 369 IYs = ys - s; 370 IYe = ye + s; 371 Ys = ys - s; 372 Ye = ye + s; 373 } 374 375 /* stencil length in each direction */ 376 s_x = s; 377 s_y = s; 378 379 /* determine starting point of each processor */ 380 nn = x*y; 381 ierr = PetscMalloc2(size+1,&bases,size,&ldims);CHKERRQ(ierr); 382 ierr = MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);CHKERRQ(ierr); 383 bases[0] = 0; 384 for (i=1; i<=size; i++) { 385 bases[i] = ldims[i-1]; 386 } 387 for (i=1; i<=size; i++) { 388 bases[i] += bases[i-1]; 389 } 390 base = bases[rank]*dof; 391 392 /* allocate the base parallel and sequential vectors */ 393 dd->Nlocal = x*y*dof; 394 ierr = VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,NULL,&global);CHKERRQ(ierr); 395 dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof; 396 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,NULL,&local);CHKERRQ(ierr); 397 398 /* generate global to local vector scatter and local to global mapping*/ 399 400 /* global to local must include ghost points within the domain, 401 but not ghost points outside the domain that aren't periodic */ 402 ierr = PetscMalloc1((IXe-IXs)*(IYe-IYs),&idx);CHKERRQ(ierr); 403 if (stencil_type == DMDA_STENCIL_BOX) { 404 left = IXs - Xs; right = left + (IXe-IXs); 405 down = IYs - Ys; up = down + (IYe-IYs); 406 count = 0; 407 for (i=down; i<up; i++) { 408 for (j=left; j<right; j++) { 409 idx[count++] = j + i*(Xe-Xs); 410 } 411 } 412 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 413 414 } else { 415 /* must drop into cross shape region */ 416 /* ---------| 417 | top | 418 |--- ---| up 419 | middle | 420 | | 421 ---- ---- down 422 | bottom | 423 ----------- 424 Xs xs xe Xe */ 425 left = xs - Xs; right = left + x; 426 down = ys - Ys; up = down + y; 427 count = 0; 428 /* bottom */ 429 for (i=(IYs-Ys); i<down; i++) { 430 for (j=left; j<right; j++) { 431 idx[count++] = j + i*(Xe-Xs); 432 } 433 } 434 /* middle */ 435 for (i=down; i<up; i++) { 436 for (j=(IXs-Xs); j<(IXe-Xs); j++) { 437 idx[count++] = j + i*(Xe-Xs); 438 } 439 } 440 /* top */ 441 for (i=up; i<up+IYe-ye; i++) { 442 for (j=left; j<right; j++) { 443 idx[count++] = j + i*(Xe-Xs); 444 } 445 } 446 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 447 } 448 449 450 /* determine who lies on each side of us stored in n6 n7 n8 451 n3 n5 452 n0 n1 n2 453 */ 454 455 /* Assume the Non-Periodic Case */ 456 n1 = rank - m; 457 if (rank % m) { 458 n0 = n1 - 1; 459 } else { 460 n0 = -1; 461 } 462 if ((rank+1) % m) { 463 n2 = n1 + 1; 464 n5 = rank + 1; 465 n8 = rank + m + 1; if (n8 >= m*n) n8 = -1; 466 } else { 467 n2 = -1; n5 = -1; n8 = -1; 468 } 469 if (rank % m) { 470 n3 = rank - 1; 471 n6 = n3 + m; if (n6 >= m*n) n6 = -1; 472 } else { 473 n3 = -1; n6 = -1; 474 } 475 n7 = rank + m; if (n7 >= m*n) n7 = -1; 476 477 if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) { 478 /* Modify for Periodic Cases */ 479 /* Handle all four corners */ 480 if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1; 481 if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0; 482 if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m; 483 if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1; 484 485 /* Handle Top and Bottom Sides */ 486 if (n1 < 0) n1 = rank + m * (n-1); 487 if (n7 < 0) n7 = rank - m * (n-1); 488 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 489 if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1; 490 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 491 if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1; 492 493 /* Handle Left and Right Sides */ 494 if (n3 < 0) n3 = rank + (m-1); 495 if (n5 < 0) n5 = rank - (m-1); 496 if ((n1 >= 0) && (n0 < 0)) n0 = rank-1; 497 if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1; 498 if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1; 499 if ((n7 >= 0) && (n8 < 0)) n8 = rank+1; 500 } else if (by == DM_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */ 501 if (n1 < 0) n1 = rank + m * (n-1); 502 if (n7 < 0) n7 = rank - m * (n-1); 503 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 504 if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1; 505 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 506 if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1; 507 } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */ 508 if (n3 < 0) n3 = rank + (m-1); 509 if (n5 < 0) n5 = rank - (m-1); 510 if ((n1 >= 0) && (n0 < 0)) n0 = rank-1; 511 if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1; 512 if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1; 513 if ((n7 >= 0) && (n8 < 0)) n8 = rank+1; 514 } 515 516 ierr = PetscMalloc1(9,&dd->neighbors);CHKERRQ(ierr); 517 518 dd->neighbors[0] = n0; 519 dd->neighbors[1] = n1; 520 dd->neighbors[2] = n2; 521 dd->neighbors[3] = n3; 522 dd->neighbors[4] = rank; 523 dd->neighbors[5] = n5; 524 dd->neighbors[6] = n6; 525 dd->neighbors[7] = n7; 526 dd->neighbors[8] = n8; 527 528 if (stencil_type == DMDA_STENCIL_STAR) { 529 /* save corner processor numbers */ 530 sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8; 531 n0 = n2 = n6 = n8 = -1; 532 } 533 534 ierr = PetscMalloc1((Xe-Xs)*(Ye-Ys),&idx);CHKERRQ(ierr); 535 536 nn = 0; 537 xbase = bases[rank]; 538 for (i=1; i<=s_y; i++) { 539 if (n0 >= 0) { /* left below */ 540 x_t = lx[n0 % m]; 541 y_t = ly[(n0/m)]; 542 s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x; 543 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 544 } 545 546 if (n1 >= 0) { /* directly below */ 547 x_t = x; 548 y_t = ly[(n1/m)]; 549 s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t; 550 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 551 } else if (by == DM_BOUNDARY_MIRROR) { 552 for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(s_y - i + 1) + j; 553 } 554 555 if (n2 >= 0) { /* right below */ 556 x_t = lx[n2 % m]; 557 y_t = ly[(n2/m)]; 558 s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t; 559 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 560 } 561 } 562 563 for (i=0; i<y; i++) { 564 if (n3 >= 0) { /* directly left */ 565 x_t = lx[n3 % m]; 566 /* y_t = y; */ 567 s_t = bases[n3] + (i+1)*x_t - s_x; 568 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 569 } else if (bx == DM_BOUNDARY_MIRROR) { 570 for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*i + s_x - j; 571 } 572 573 for (j=0; j<x; j++) idx[nn++] = xbase++; /* interior */ 574 575 if (n5 >= 0) { /* directly right */ 576 x_t = lx[n5 % m]; 577 /* y_t = y; */ 578 s_t = bases[n5] + (i)*x_t; 579 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 580 } else if (bx == DM_BOUNDARY_MIRROR) { 581 for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*(i + 1) - 2 - j; 582 } 583 } 584 585 for (i=1; i<=s_y; i++) { 586 if (n6 >= 0) { /* left above */ 587 x_t = lx[n6 % m]; 588 /* y_t = ly[(n6/m)]; */ 589 s_t = bases[n6] + (i)*x_t - s_x; 590 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 591 } 592 593 if (n7 >= 0) { /* directly above */ 594 x_t = x; 595 /* y_t = ly[(n7/m)]; */ 596 s_t = bases[n7] + (i-1)*x_t; 597 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 598 } else if (by == DM_BOUNDARY_MIRROR) { 599 for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(y - i - 1) + j; 600 } 601 602 if (n8 >= 0) { /* right above */ 603 x_t = lx[n8 % m]; 604 /* y_t = ly[(n8/m)]; */ 605 s_t = bases[n8] + (i-1)*x_t; 606 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 607 } 608 } 609 610 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);CHKERRQ(ierr); 611 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 612 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);CHKERRQ(ierr); 613 ierr = ISDestroy(&to);CHKERRQ(ierr); 614 ierr = ISDestroy(&from);CHKERRQ(ierr); 615 616 if (stencil_type == DMDA_STENCIL_STAR) { 617 n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8; 618 } 619 620 if (((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC))) { 621 /* 622 Recompute the local to global mappings, this time keeping the 623 information about the cross corner processor numbers and any ghosted 624 but not periodic indices. 625 */ 626 nn = 0; 627 xbase = bases[rank]; 628 for (i=1; i<=s_y; i++) { 629 if (n0 >= 0) { /* left below */ 630 x_t = lx[n0 % m]; 631 y_t = ly[(n0/m)]; 632 s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x; 633 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 634 } else if (xs-Xs > 0 && ys-Ys > 0) { 635 for (j=0; j<s_x; j++) idx[nn++] = -1; 636 } 637 if (n1 >= 0) { /* directly below */ 638 x_t = x; 639 y_t = ly[(n1/m)]; 640 s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t; 641 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 642 } else if (ys-Ys > 0) { 643 if (by == DM_BOUNDARY_MIRROR) { 644 for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(s_y - i + 1) + j; 645 } else { 646 for (j=0; j<x; j++) idx[nn++] = -1; 647 } 648 } 649 if (n2 >= 0) { /* right below */ 650 x_t = lx[n2 % m]; 651 y_t = ly[(n2/m)]; 652 s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t; 653 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 654 } else if (Xe-xe> 0 && ys-Ys > 0) { 655 for (j=0; j<s_x; j++) idx[nn++] = -1; 656 } 657 } 658 659 for (i=0; i<y; i++) { 660 if (n3 >= 0) { /* directly left */ 661 x_t = lx[n3 % m]; 662 /* y_t = y; */ 663 s_t = bases[n3] + (i+1)*x_t - s_x; 664 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 665 } else if (xs-Xs > 0) { 666 if (bx == DM_BOUNDARY_MIRROR) { 667 for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*i + s_x - j; 668 } else { 669 for (j=0; j<s_x; j++) idx[nn++] = -1; 670 } 671 } 672 673 for (j=0; j<x; j++) idx[nn++] = xbase++; /* interior */ 674 675 if (n5 >= 0) { /* directly right */ 676 x_t = lx[n5 % m]; 677 /* y_t = y; */ 678 s_t = bases[n5] + (i)*x_t; 679 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 680 } else if (Xe-xe > 0) { 681 if (bx == DM_BOUNDARY_MIRROR) { 682 for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*(i + 1) - 2 - j; 683 } else { 684 for (j=0; j<s_x; j++) idx[nn++] = -1; 685 } 686 } 687 } 688 689 for (i=1; i<=s_y; i++) { 690 if (n6 >= 0) { /* left above */ 691 x_t = lx[n6 % m]; 692 /* y_t = ly[(n6/m)]; */ 693 s_t = bases[n6] + (i)*x_t - s_x; 694 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 695 } else if (xs-Xs > 0 && Ye-ye > 0) { 696 for (j=0; j<s_x; j++) idx[nn++] = -1; 697 } 698 if (n7 >= 0) { /* directly above */ 699 x_t = x; 700 /* y_t = ly[(n7/m)]; */ 701 s_t = bases[n7] + (i-1)*x_t; 702 for (j=0; j<x_t; j++) idx[nn++] = s_t++; 703 } else if (Ye-ye > 0) { 704 if (by == DM_BOUNDARY_MIRROR) { 705 for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(y - i - 1) + j; 706 } else { 707 for (j=0; j<x; j++) idx[nn++] = -1; 708 } 709 } 710 if (n8 >= 0) { /* right above */ 711 x_t = lx[n8 % m]; 712 /* y_t = ly[(n8/m)]; */ 713 s_t = bases[n8] + (i-1)*x_t; 714 for (j=0; j<s_x; j++) idx[nn++] = s_t++; 715 } else if (Xe-xe > 0 && Ye-ye > 0) { 716 for (j=0; j<s_x; j++) idx[nn++] = -1; 717 } 718 } 719 } 720 /* 721 Set the local to global ordering in the global vector, this allows use 722 of VecSetValuesLocal(). 723 */ 724 ierr = ISLocalToGlobalMappingCreate(comm,dof,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);CHKERRQ(ierr); 725 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)da->ltogmap);CHKERRQ(ierr); 726 727 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 728 dd->m = m; dd->n = n; 729 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 730 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1; 731 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1; 732 733 ierr = VecDestroy(&local);CHKERRQ(ierr); 734 ierr = VecDestroy(&global);CHKERRQ(ierr); 735 736 dd->gtol = gtol; 737 dd->base = base; 738 da->ops->view = DMView_DA_2d; 739 dd->ltol = NULL; 740 dd->ao = NULL; 741 PetscFunctionReturn(0); 742 } 743 744 /*@C 745 DMDACreate2d - Creates an object that will manage the communication of two-dimensional 746 regular array data that is distributed across some processors. 747 748 Collective 749 750 Input Parameters: 751 + comm - MPI communicator 752 . bx,by - type of ghost nodes the array have. 753 Use one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC. 754 . stencil_type - stencil type. Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR. 755 . M,N - global dimension in each direction of the array 756 . m,n - corresponding number of processors in each dimension 757 (or PETSC_DECIDE to have calculated) 758 . dof - number of degrees of freedom per node 759 . s - stencil width 760 - lx, ly - arrays containing the number of nodes in each cell along 761 the x and y coordinates, or NULL. If non-null, these 762 must be of length as m and n, and the corresponding 763 m and n cannot be PETSC_DECIDE. The sum of the lx[] entries 764 must be M, and the sum of the ly[] entries must be N. 765 766 Output Parameter: 767 . da - the resulting distributed array object 768 769 Options Database Key: 770 + -dm_view - Calls DMView() at the conclusion of DMDACreate2d() 771 . -da_grid_x <nx> - number of grid points in x direction 772 . -da_grid_y <ny> - number of grid points in y direction 773 . -da_processors_x <nx> - number of processors in x direction 774 . -da_processors_y <ny> - number of processors in y direction 775 . -da_refine_x <rx> - refinement ratio in x direction 776 . -da_refine_y <ry> - refinement ratio in y direction 777 - -da_refine <n> - refine the DMDA n times before creating 778 779 780 Level: beginner 781 782 Notes: 783 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 784 standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 785 the standard 9-pt stencil. 786 787 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 788 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 789 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 790 791 You must call DMSetUp() after this call before using this DM. 792 793 If you wish to use the options database to change values in the DMDA call DMSetFromOptions() after this call 794 but before DMSetUp(). 795 796 .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(), 797 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDASetRefinementFactor(), 798 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges(), 799 DMStagCreate2d() 800 801 @*/ 802 803 PetscErrorCode DMDACreate2d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMDAStencilType stencil_type, 804 PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da) 805 { 806 PetscErrorCode ierr; 807 808 PetscFunctionBegin; 809 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 810 ierr = DMSetDimension(*da, 2);CHKERRQ(ierr); 811 ierr = DMDASetSizes(*da, M, N, 1);CHKERRQ(ierr); 812 ierr = DMDASetNumProcs(*da, m, n, PETSC_DECIDE);CHKERRQ(ierr); 813 ierr = DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE);CHKERRQ(ierr); 814 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 815 ierr = DMDASetStencilType(*da, stencil_type);CHKERRQ(ierr); 816 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 817 ierr = DMDASetOwnershipRanges(*da, lx, ly, NULL);CHKERRQ(ierr); 818 PetscFunctionReturn(0); 819 } 820