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