1 2 /* 3 Code for manipulating distributed regular 1d arrays in parallel. 4 This file was created by Peter Mell 6/30/95 5 */ 6 7 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 8 9 const char *const DMDABoundaryTypes[] = {"BOUNDARY_NONE","BOUNDARY_GHOSTED","BOUNDARY_PERIODIC","DMDA_",0}; 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "DMView_DA_1d" 13 PetscErrorCode DMView_DA_1d(DM da,PetscViewer viewer) 14 { 15 PetscErrorCode ierr; 16 PetscMPIInt rank; 17 PetscBool iascii,isdraw,isbinary; 18 DM_DA *dd = (DM_DA*)da->data; 19 #if defined(PETSC_HAVE_MATLAB_ENGINE) 20 PetscBool ismatlab; 21 #endif 22 23 PetscFunctionBegin; 24 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 25 26 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 27 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 28 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 29 #if defined(PETSC_HAVE_MATLAB_ENGINE) 30 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);CHKERRQ(ierr); 31 #endif 32 if (iascii) { 33 PetscViewerFormat format; 34 35 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 36 if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) { 37 DMDALocalInfo info; 38 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 39 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 40 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D m %D w %D s %D\n",rank,dd->M,dd->m,dd->w,dd->s);CHKERRQ(ierr); 41 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D\n",info.xs,info.xs+info.xm);CHKERRQ(ierr); 42 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 43 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 44 } else { 45 ierr = DMView_DA_VTK(da, viewer);CHKERRQ(ierr); 46 } 47 } else if (isdraw) { 48 PetscDraw draw; 49 double ymin = -1,ymax = 1,xmin = -1,xmax = dd->M,x; 50 PetscInt base; 51 char node[10]; 52 PetscBool isnull; 53 54 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 55 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 56 57 ierr = PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);CHKERRQ(ierr); 58 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 59 60 /* first processor draws all node lines */ 61 if (!rank) { 62 PetscInt xmin_tmp; 63 ymin = 0.0; ymax = 0.3; 64 65 for (xmin_tmp=0; xmin_tmp < dd->M; xmin_tmp++) { 66 ierr = PetscDrawLine(draw,(double)xmin_tmp,ymin,(double)xmin_tmp,ymax,PETSC_DRAW_BLACK);CHKERRQ(ierr); 67 } 68 69 xmin = 0.0; xmax = dd->M - 1; 70 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);CHKERRQ(ierr); 71 ierr = PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_BLACK);CHKERRQ(ierr); 72 } 73 74 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 75 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 76 77 /* draw my box */ 78 ymin = 0; ymax = 0.3; xmin = dd->xs / dd->w; xmax = (dd->xe / dd->w) - 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 84 /* Put in index numbers */ 85 base = dd->base / dd->w; 86 for (x=xmin; x<=xmax; x++) { 87 sprintf(node,"%d",(int)base++); 88 ierr = PetscDrawString(draw,x,ymin,PETSC_DRAW_RED,node);CHKERRQ(ierr); 89 } 90 91 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 92 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 93 } else if (isbinary){ 94 ierr = DMView_DA_Binary(da,viewer);CHKERRQ(ierr); 95 #if defined(PETSC_HAVE_MATLAB_ENGINE) 96 } else if (ismatlab) { 97 ierr = DMView_DA_Matlab(da,viewer);CHKERRQ(ierr); 98 #endif 99 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for DMDA 1d",((PetscObject)viewer)->type_name); 100 PetscFunctionReturn(0); 101 } 102 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "DMSetUp_DA_1D" 106 PetscErrorCode DMSetUp_DA_1D(DM da) 107 { 108 DM_DA *dd = (DM_DA*)da->data; 109 const PetscInt M = dd->M; 110 const PetscInt dof = dd->w; 111 const PetscInt s = dd->s; 112 const PetscInt o = dd->overlap; 113 const PetscInt sDist = s*dof; /* absolute stencil distance */ 114 const PetscInt oDist = o*dof; 115 const PetscInt *lx = dd->lx; 116 DMDABoundaryType bx = dd->bx; 117 MPI_Comm comm; 118 Vec local, global; 119 VecScatter ltog, gtol; 120 IS to, from; 121 PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE; 122 PetscMPIInt rank, size; 123 PetscInt i,j,*idx,nn,left,xs,xe,x,Xs,Xe,start,end,m,IXs,IXe; 124 PetscErrorCode ierr; 125 126 PetscFunctionBegin; 127 ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr); 128 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 129 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 130 131 dd->m = size; 132 m = dd->m; 133 134 if (s > 0) { 135 /* if not communicating data then should be ok to have nothing on some processes */ 136 if (M < m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"More processes than data points! %D %D",m,M); 137 if ((M-1) < s) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Array is too small for stencil! %D %D",M-1,s); 138 } 139 140 /* 141 Determine locally owned region 142 xs is the first local node number, x is the number of local nodes 143 */ 144 if (!lx) { 145 ierr = PetscMalloc(m*sizeof(PetscInt), &dd->lx);CHKERRQ(ierr); 146 ierr = PetscOptionsGetBool(PETSC_NULL,"-da_partition_blockcomm",&flg1,PETSC_NULL);CHKERRQ(ierr); 147 ierr = PetscOptionsGetBool(PETSC_NULL,"-da_partition_nodes_at_end",&flg2,PETSC_NULL);CHKERRQ(ierr); 148 if (flg1) { /* Block Comm type Distribution */ 149 xs = rank*M/m; 150 x = (rank + 1)*M/m - xs; 151 } else if (flg2) { /* The odd nodes are evenly distributed across last nodes */ 152 x = (M + rank)/m; 153 if (M/m == x) { xs = rank*x; } 154 else { xs = rank*(x-1) + (M+rank)%(x*m); } 155 } else { /* The odd nodes are evenly distributed across the first k nodes */ 156 /* Regular PETSc Distribution */ 157 x = M/m + ((M % m) > rank); 158 if (rank >= (M % m)) {xs = (rank * (PetscInt)(M/m) + M % m);} 159 else {xs = rank * (PetscInt)(M/m) + rank;} 160 } 161 ierr = MPI_Allgather(&xs,1,MPIU_INT,dd->lx,1,MPIU_INT,comm);CHKERRQ(ierr); 162 for (i=0; i<m-1; i++) dd->lx[i] = dd->lx[i+1] - dd->lx[i]; 163 dd->lx[m-1] = M - dd->lx[m-1]; 164 } else { 165 x = lx[rank]; 166 xs = 0; 167 for (i=0; i<rank; i++) { 168 xs += lx[i]; 169 } 170 /* verify that data user provided is consistent */ 171 left = xs; 172 for (i=rank; i<size; i++) { 173 left += lx[i]; 174 } 175 if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M %D %D",left,M); 176 } 177 178 /* 179 check if the scatter requires more than one process neighbor or wraps around 180 the domain more than once 181 */ 182 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); 183 184 /* From now on x,xs,xe,Xs,Xe are the exact location in the array */ 185 x *= dof; 186 xs *= dof; 187 xe = xs + x; 188 189 /* determine ghost region (Xs) and region scattered into (IXs) */ 190 if (xs-sDist-oDist > 0) { 191 Xs = xs - sDist - oDist; 192 IXs = xs - sDist - oDist; 193 } else { 194 if (bx) { 195 Xs = xs - sDist; 196 } else { 197 Xs = 0; 198 } 199 IXs = 0; 200 } 201 if (xe+sDist+oDist <= M*dof) { 202 Xe = xe + sDist + oDist; 203 IXe = xe + sDist + oDist; 204 } else { 205 if (bx) { 206 Xe = xe + sDist; 207 } else { 208 Xe = M*dof; 209 } 210 IXe = M*dof; 211 } 212 213 if (bx == DMDA_BOUNDARY_PERIODIC || bx == DMDA_BOUNDARY_MIRROR) { 214 Xs = xs - sDist - oDist; 215 Xe = xe + sDist + oDist; 216 IXs = xs - sDist - oDist; 217 IXe = xe + sDist + oDist; 218 } 219 220 /* allocate the base parallel and sequential vectors */ 221 dd->Nlocal = x; 222 ierr = VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);CHKERRQ(ierr); 223 dd->nlocal = (Xe-Xs); 224 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);CHKERRQ(ierr); 225 226 /* Create Local to Global Vector Scatter Context */ 227 /* local to global inserts non-ghost point region into global */ 228 ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr); 229 ierr = ISCreateStride(comm,x,start,1,&to);CHKERRQ(ierr); 230 ierr = ISCreateStride(comm,x,xs-Xs,1,&from);CHKERRQ(ierr); 231 ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr); 232 ierr = PetscLogObjectParent(da,ltog);CHKERRQ(ierr); 233 ierr = ISDestroy(&from);CHKERRQ(ierr); 234 ierr = ISDestroy(&to);CHKERRQ(ierr); 235 236 /* Create Global to Local Vector Scatter Context */ 237 /* global to local must retrieve ghost points */ 238 ierr = ISCreateStride(comm,(IXe-IXs),IXs-Xs,1,&to);CHKERRQ(ierr); 239 240 ierr = PetscMalloc((x+2*(sDist+oDist))*sizeof(PetscInt),&idx);CHKERRQ(ierr); 241 ierr = PetscLogObjectMemory(da,(x+2*(sDist+oDist))*sizeof(PetscInt));CHKERRQ(ierr); 242 243 for (i=0; i<IXs-Xs; i++) {idx[i] = -1; } /* prepend with -1s if needed for ghosted case*/ 244 245 nn = IXs-Xs; 246 if (bx == DMDA_BOUNDARY_PERIODIC) { /* Handle all cases with periodic first */ 247 for (i=0; i<sDist+oDist; i++) { /* Left ghost points */ 248 if ((xs-sDist-oDist+i)>=0) {idx[nn++] = xs-sDist-oDist+i;} 249 else {idx[nn++] = M*dof+(xs-sDist-oDist+i);} 250 } 251 252 for (i=0; i<x; i++) { idx [nn++] = xs + i;} /* Non-ghost points */ 253 254 for (i=0; i<sDist+oDist; i++) { /* Right ghost points */ 255 if ((xe+i)<M*dof) {idx [nn++] = xe+i; } 256 else {idx [nn++] = (xe+i) - M*dof;} 257 } 258 } else if (bx == DMDA_BOUNDARY_MIRROR) { /* Handle all cases with periodic first */ 259 for (i=0; i<(sDist+oDist)/dof; i++) { /* Left ghost points */ 260 for (j=0; j<dof; j++) { 261 if ((xs-sDist-oDist+i*dof + j)>=0) {idx[nn++] = xs-sDist-oDist+i*dof +j;} 262 else {idx[nn++] = sDist+oDist - dof*(i) + j;} 263 } 264 } 265 266 for (i=0; i<x; i++) { idx [nn++] = xs + i;} /* Non-ghost points */ 267 268 for (i=0; i<(sDist+oDist)/dof; i++) { /* Right ghost points */ 269 for (j=0; j<dof; j++) { 270 if ((xe+i)<M*dof) {idx[nn++] = xe+i*dof+j; } 271 else {idx[nn++] = M*dof - dof*(i + 2) + j ;} 272 } 273 } 274 } else { /* Now do all cases with no periodicity */ 275 if (0 <= xs-sDist-oDist) {for (i=0; i<sDist+oDist; i++) {idx[nn++] = xs - sDist - oDist + i;}} 276 else {for (i=0; i<xs; i++) {idx[nn++] = i;}} 277 278 for (i=0; i<x; i++) { idx [nn++] = xs + i;} 279 280 if ((xe+sDist+oDist)<=M*dof) {for (i=0; i<sDist+oDist; i++) {idx[nn++]=xe+i;}} 281 else {for (i=xe; i<(M*dof); i++) {idx[nn++]=i;}} 282 } 283 284 ierr = ISCreateGeneral(comm,nn-IXs+Xs,&idx[IXs-Xs],PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 285 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 286 ierr = PetscLogObjectParent(da,to);CHKERRQ(ierr); 287 ierr = PetscLogObjectParent(da,from);CHKERRQ(ierr); 288 ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr); 289 ierr = ISDestroy(&to);CHKERRQ(ierr); 290 ierr = ISDestroy(&from);CHKERRQ(ierr); 291 ierr = VecDestroy(&local);CHKERRQ(ierr); 292 ierr = VecDestroy(&global);CHKERRQ(ierr); 293 294 dd->xs = xs; dd->xe = xe; dd->ys = 0; dd->ye = 1; dd->zs = 0; dd->ze = 1; 295 dd->Xs = Xs; dd->Xe = Xe; dd->Ys = 0; dd->Ye = 1; dd->Zs = 0; dd->Ze = 1; 296 297 dd->gtol = gtol; 298 dd->ltog = ltog; 299 dd->base = xs; 300 da->ops->view = DMView_DA_1d; 301 302 /* 303 Set the local to global ordering in the global vector, this allows use 304 of VecSetValuesLocal(). 305 */ 306 for (i=0; i<Xe-IXe; i++) {idx[nn++] = -1; } /* pad with -1s if needed for ghosted case*/ 307 308 ierr = ISLocalToGlobalMappingCreate(comm,nn,idx,PETSC_COPY_VALUES,&da->ltogmap);CHKERRQ(ierr); 309 ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr); 310 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 311 312 dd->idx = idx; 313 dd->Nl = nn; 314 315 PetscFunctionReturn(0); 316 } 317 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "DMDACreate1d" 321 /*@C 322 DMDACreate1d - Creates an object that will manage the communication of one-dimensional 323 regular array data that is distributed across some processors. 324 325 Collective on MPI_Comm 326 327 Input Parameters: 328 + comm - MPI communicator 329 . bx - type of ghost cells at the boundary the array should have, if any. Use 330 DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, or DMDA_BOUNDARY_PERIODIC. 331 . M - global dimension of the array (use -M to indicate that it may be set to a different value 332 from the command line with -da_grid_x <M>) 333 . dof - number of degrees of freedom per node 334 . s - stencil width 335 - lx - array containing number of nodes in the X direction on each processor, 336 or PETSC_NULL. If non-null, must be of length as the number of processes in the MPI_Comm. 337 338 Output Parameter: 339 . da - the resulting distributed array object 340 341 Options Database Key: 342 + -dm_view - Calls DMView() at the conclusion of DMDACreate1d() 343 . -da_grid_x <nx> - number of grid points in x direction; can set if M < 0 344 . -da_refine_x <rx> - refinement factor 345 - -da_refine <n> - refine the DMDA n times before creating it, if M < 0 346 347 Level: beginner 348 349 Notes: 350 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 351 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 352 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 353 354 .keywords: distributed array, create, one-dimensional 355 356 .seealso: DMDestroy(), DMView(), DMDACreate2d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDASetRefinementFactor(), 357 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDAGetRefinementFactor(), 358 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges() 359 360 @*/ 361 PetscErrorCode DMDACreate1d(MPI_Comm comm, DMDABoundaryType bx, PetscInt M, PetscInt dof, PetscInt s, const PetscInt lx[], DM *da) 362 { 363 PetscErrorCode ierr; 364 PetscMPIInt size; 365 366 PetscFunctionBegin; 367 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 368 ierr = DMDASetDim(*da, 1);CHKERRQ(ierr); 369 ierr = DMDASetSizes(*da, M, 1, 1);CHKERRQ(ierr); 370 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 371 ierr = DMDASetNumProcs(*da, size, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr); 372 ierr = DMDASetBoundaryType(*da, bx, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE);CHKERRQ(ierr); 373 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 374 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 375 ierr = DMDASetOwnershipRanges(*da, lx, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 376 /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */ 377 ierr = DMSetFromOptions(*da);CHKERRQ(ierr); 378 ierr = DMSetUp(*da);CHKERRQ(ierr); 379 ierr = DMViewFromOptions(*da,"-dm_view");CHKERRQ(ierr); 380 PetscFunctionReturn(0); 381 } 382