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