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