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