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