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