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