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