1 #define PETSCDM_DLL 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 "petscda.h" I*/ 8 9 const char *DAPeriodicTypes[] = {"NONPERIODIC","XPERIODIC","YPERIODIC","XYPERIODIC", 10 "XYZPERIODIC","XZPERIODIC","YZPERIODIC","ZPERIODIC","XYZGHOSTED","DAPeriodicType","DA_",0}; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "DAView_1d" 14 PetscErrorCode DAView_1d(DA da,PetscViewer viewer) 15 { 16 PetscErrorCode ierr; 17 PetscMPIInt rank; 18 PetscBool iascii,isdraw; 19 DM_DA *dd = (DM_DA*)da->data; 20 21 PetscFunctionBegin; 22 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 23 24 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 25 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 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 DALocalInfo info; 32 ierr = DAGetLocalInfo(da,&info);CHKERRQ(ierr); 33 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); 34 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D\n",info.xs,info.xs+info.xm);CHKERRQ(ierr); 35 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 36 } 37 } else if (isdraw) { 38 PetscDraw draw; 39 double ymin = -1,ymax = 1,xmin = -1,xmax = dd->M,x; 40 PetscInt base; 41 char node[10]; 42 PetscBool isnull; 43 44 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 45 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 46 47 ierr = PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);CHKERRQ(ierr); 48 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 49 50 /* first processor draws all node lines */ 51 if (!rank) { 52 PetscInt xmin_tmp; 53 ymin = 0.0; ymax = 0.3; 54 55 /* ADIC doesn't like doubles in a for loop */ 56 for (xmin_tmp =0; xmin_tmp < dd->M; xmin_tmp++) { 57 ierr = PetscDrawLine(draw,(double)xmin_tmp,ymin,(double)xmin_tmp,ymax,PETSC_DRAW_BLACK);CHKERRQ(ierr); 58 } 59 60 xmin = 0.0; xmax = dd->M - 1; 61 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);CHKERRQ(ierr); 62 ierr = PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_BLACK);CHKERRQ(ierr); 63 } 64 65 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 66 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 67 68 /* draw my box */ 69 ymin = 0; ymax = 0.3; xmin = dd->xs / dd->w; xmax = (dd->xe / dd->w) - 1; 70 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);CHKERRQ(ierr); 71 ierr = PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 72 ierr = PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 73 ierr = PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 74 75 /* Put in index numbers */ 76 base = dd->base / dd->w; 77 for (x=xmin; x<=xmax; x++) { 78 sprintf(node,"%d",(int)base++); 79 ierr = PetscDrawString(draw,x,ymin,PETSC_DRAW_RED,node);CHKERRQ(ierr); 80 } 81 82 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 83 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 84 } else { 85 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for DA 1d",((PetscObject)viewer)->type_name); 86 } 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "DAView_Private" 92 /* 93 Processes command line options to determine if/how a DA 94 is to be viewed. Called by DACreateXX() 95 */ 96 PetscErrorCode DAView_Private(DA da) 97 { 98 PetscErrorCode ierr; 99 PetscBool flg1 = PETSC_FALSE; 100 PetscViewer view; 101 102 PetscFunctionBegin; 103 ierr = PetscOptionsBegin(((PetscObject)da)->comm,((PetscObject)da)->prefix,"DA viewing options","DA");CHKERRQ(ierr); 104 ierr = PetscOptionsTruth("-da_view","Print information about the DA's distribution","DAView",PETSC_FALSE,&flg1,PETSC_NULL);CHKERRQ(ierr); 105 if (flg1) { 106 ierr = PetscViewerASCIIGetStdout(((PetscObject)da)->comm,&view);CHKERRQ(ierr); 107 ierr = DAView(da,view);CHKERRQ(ierr); 108 } 109 flg1 = PETSC_FALSE; 110 ierr = PetscOptionsTruth("-da_view_draw","Draw how the DA is distributed","DAView",PETSC_FALSE,&flg1,PETSC_NULL);CHKERRQ(ierr); 111 if (flg1) {ierr = DAView(da,PETSC_VIEWER_DRAW_(((PetscObject)da)->comm));CHKERRQ(ierr);} 112 ierr = PetscOptionsEnd();CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 EXTERN_C_BEGIN 117 #undef __FUNCT__ 118 #define __FUNCT__ "DACreate_1D" 119 PetscErrorCode PETSCDM_DLLEXPORT DACreate_1D(DA da) 120 { 121 DM_DA *dd = (DM_DA*)da->data; 122 const PetscInt dim = dd->dim; 123 const PetscInt M = dd->M; 124 const PetscInt dof = dd->w; 125 const PetscInt s = dd->s; 126 const PetscInt sDist = s*dof; /* absolute stencil distance */ 127 const PetscInt *lx = dd->lx; 128 const DAPeriodicType wrap = dd->wrap; 129 MPI_Comm comm; 130 Vec local, global; 131 VecScatter ltog, gtol; 132 IS to, from; 133 PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE; 134 PetscMPIInt rank, size; 135 PetscInt i,*idx,nn,left,xs,xe,x,Xs,Xe,start,end,m; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 if (dim != PETSC_DECIDE && dim != 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Dimension should be 1: %D",dim); 140 if (dof < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof); 141 if (s < 0) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s); 142 143 dd->dim = 1; 144 ierr = PetscMalloc(dof*sizeof(char*),&dd->fieldname);CHKERRQ(ierr); 145 ierr = PetscMemzero(dd->fieldname,dof*sizeof(char*));CHKERRQ(ierr); 146 ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr); 147 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 148 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 149 150 dd->m = size; 151 m = dd->m; 152 153 if (s > 0) { 154 /* if not communicating data then should be ok to have nothing on some processes */ 155 if (M < m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"More processes than data points! %D %D",m,M); 156 if ((M-1) < s) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Array is too small for stencil! %D %D",M-1,s); 157 } 158 159 /* 160 Determine locally owned region 161 xs is the first local node number, x is the number of local nodes 162 */ 163 if (!lx) { 164 ierr = PetscOptionsGetTruth(PETSC_NULL,"-da_partition_blockcomm",&flg1,PETSC_NULL);CHKERRQ(ierr); 165 ierr = PetscOptionsGetTruth(PETSC_NULL,"-da_partition_nodes_at_end",&flg2,PETSC_NULL);CHKERRQ(ierr); 166 if (flg1) { /* Block Comm type Distribution */ 167 xs = rank*M/m; 168 x = (rank + 1)*M/m - xs; 169 } else if (flg2) { /* The odd nodes are evenly distributed across last nodes */ 170 x = (M + rank)/m; 171 if (M/m == x) { xs = rank*x; } 172 else { xs = rank*(x-1) + (M+rank)%(x*m); } 173 } else { /* The odd nodes are evenly distributed across the first k nodes */ 174 /* Regular PETSc Distribution */ 175 x = M/m + ((M % m) > rank); 176 if (rank >= (M % m)) {xs = (rank * (PetscInt)(M/m) + M % m);} 177 else {xs = rank * (PetscInt)(M/m) + rank;} 178 } 179 } else { 180 x = lx[rank]; 181 xs = 0; 182 for (i=0; i<rank; i++) { 183 xs += lx[i]; 184 } 185 /* verify that data user provided is consistent */ 186 left = xs; 187 for (i=rank; i<size; i++) { 188 left += lx[i]; 189 } 190 if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M %D %D",left,M); 191 } 192 193 /* From now on x,xs,xe,Xs,Xe are the exact location in the array */ 194 x *= dof; 195 xs *= dof; 196 xe = xs + x; 197 198 /* determine ghost region */ 199 if (wrap == DA_XPERIODIC || wrap == DA_XYZGHOSTED) { 200 Xs = xs - sDist; 201 Xe = xe + sDist; 202 } else { 203 if ((xs-sDist) >= 0) Xs = xs-sDist; else Xs = 0; 204 if ((xe+sDist) <= M*dof) Xe = xe+sDist; else Xe = M*dof; 205 } 206 207 /* allocate the base parallel and sequential vectors */ 208 dd->Nlocal = x; 209 ierr = VecCreateMPIWithArray(comm,dd->Nlocal,PETSC_DECIDE,0,&global);CHKERRQ(ierr); 210 ierr = VecSetBlockSize(global,dof);CHKERRQ(ierr); 211 dd->nlocal = (Xe-Xs); 212 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dd->nlocal,0,&local);CHKERRQ(ierr); 213 ierr = VecSetBlockSize(local,dof);CHKERRQ(ierr); 214 215 /* Create Local to Global Vector Scatter Context */ 216 /* local to global inserts non-ghost point region into global */ 217 ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr); 218 ierr = ISCreateStride(comm,x,start,1,&to);CHKERRQ(ierr); 219 ierr = ISCreateStride(comm,x,xs-Xs,1,&from);CHKERRQ(ierr); 220 ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr); 221 ierr = PetscLogObjectParent(da,ltog);CHKERRQ(ierr); 222 ierr = ISDestroy(from);CHKERRQ(ierr); 223 ierr = ISDestroy(to);CHKERRQ(ierr); 224 225 /* Create Global to Local Vector Scatter Context */ 226 /* global to local must retrieve ghost points */ 227 if (wrap == DA_XYZGHOSTED) { 228 if (size == 1) { 229 ierr = ISCreateStride(comm,(xe-xs),sDist,1,&to);CHKERRQ(ierr); 230 } else if (!rank) { 231 ierr = ISCreateStride(comm,(Xe-xs),sDist,1,&to);CHKERRQ(ierr); 232 } else if (rank == size-1) { 233 ierr = ISCreateStride(comm,(xe-Xs),0,1,&to);CHKERRQ(ierr); 234 } else { 235 ierr = ISCreateStride(comm,(Xe-Xs),0,1,&to);CHKERRQ(ierr); 236 } 237 } else { 238 ierr = ISCreateStride(comm,(Xe-Xs),0,1,&to);CHKERRQ(ierr); 239 } 240 241 ierr = PetscMalloc((x+2*sDist)*sizeof(PetscInt),&idx);CHKERRQ(ierr); 242 ierr = PetscLogObjectMemory(da,(x+2*sDist)*sizeof(PetscInt));CHKERRQ(ierr); 243 244 nn = 0; 245 if (wrap == DA_XPERIODIC) { /* Handle all cases with wrap first */ 246 247 for (i=0; i<sDist; i++) { /* Left ghost points */ 248 if ((xs-sDist+i)>=0) { idx[nn++] = xs-sDist+i;} 249 else { idx[nn++] = M*dof+(xs-sDist+i);} 250 } 251 252 for (i=0; i<x; i++) { idx [nn++] = xs + i;} /* Non-ghost points */ 253 254 for (i=0; i<sDist; 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 (wrap == DA_XYZGHOSTED) { 259 260 if (sDist <= xs) {for (i=0; i<sDist; i++) {idx[nn++] = xs - sDist + i;}} 261 262 for (i=0; i<x; i++) { idx [nn++] = xs + i;} 263 264 if ((xe+sDist)<=M*dof) {for (i=0; i<sDist; i++) {idx[nn++]=xe+i;}} 265 266 } else { /* Now do all cases with no wrapping */ 267 268 if (sDist <= xs) {for (i=0; i<sDist; i++) {idx[nn++] = xs - sDist + i;}} 269 else {for (i=0; i<xs; i++) {idx[nn++] = i;}} 270 271 for (i=0; i<x; i++) { idx [nn++] = xs + i;} 272 273 if ((xe+sDist)<=M*dof) {for (i=0; i<sDist; i++) {idx[nn++]=xe+i;}} 274 else {for (i=xe; i<(M*dof); i++) {idx[nn++]=i;}} 275 } 276 277 ierr = ISCreateGeneral(comm,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 278 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 279 ierr = PetscLogObjectParent(da,to);CHKERRQ(ierr); 280 ierr = PetscLogObjectParent(da,from);CHKERRQ(ierr); 281 ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr); 282 ierr = ISDestroy(to);CHKERRQ(ierr); 283 ierr = ISDestroy(from);CHKERRQ(ierr); 284 ierr = VecDestroy(local);CHKERRQ(ierr); 285 ierr = VecDestroy(global);CHKERRQ(ierr); 286 287 dd->xs = xs; dd->xe = xe; dd->ys = 0; dd->ye = 1; dd->zs = 0; dd->ze = 1; 288 dd->Xs = Xs; dd->Xe = Xe; dd->Ys = 0; dd->Ye = 1; dd->Zs = 0; dd->Ze = 1; 289 290 dd->gtol = gtol; 291 dd->ltog = ltog; 292 dd->base = xs; 293 da->ops->view = DAView_1d; 294 295 /* 296 Set the local to global ordering in the global vector, this allows use 297 of VecSetValuesLocal(). 298 */ 299 if (wrap == DA_XYZGHOSTED) { 300 PetscInt *tmpidx; 301 if (size == 1) { 302 ierr = PetscMalloc((nn+2*sDist)*sizeof(PetscInt),&tmpidx);CHKERRQ(ierr); 303 for (i=0; i<sDist; i++) tmpidx[i] = -1; 304 ierr = PetscMemcpy(tmpidx+sDist,idx,nn*sizeof(PetscInt));CHKERRQ(ierr); 305 for (i=nn+sDist; i<nn+2*sDist; i++) tmpidx[i] = -1; 306 ierr = PetscFree(idx);CHKERRQ(ierr); 307 idx = tmpidx; 308 nn += 2*sDist; 309 } else if (!rank) { /* must preprend -1 marker for ghost location that have no global value */ 310 ierr = PetscMalloc((nn+sDist)*sizeof(PetscInt),&tmpidx);CHKERRQ(ierr); 311 for (i=0; i<sDist; i++) tmpidx[i] = -1; 312 ierr = PetscMemcpy(tmpidx+sDist,idx,nn*sizeof(PetscInt));CHKERRQ(ierr); 313 ierr = PetscFree(idx);CHKERRQ(ierr); 314 idx = tmpidx; 315 nn += sDist; 316 } else if (rank == size-1) { /* must postpend -1 marker for ghost location that have no global value */ 317 ierr = PetscMalloc((nn+sDist)*sizeof(PetscInt),&tmpidx);CHKERRQ(ierr); 318 ierr = PetscMemcpy(tmpidx,idx,nn*sizeof(PetscInt));CHKERRQ(ierr); 319 for (i=nn; i<nn+sDist; i++) tmpidx[i] = -1; 320 ierr = PetscFree(idx);CHKERRQ(ierr); 321 idx = tmpidx; 322 nn += sDist; 323 } 324 } 325 ierr = ISLocalToGlobalMappingCreate(comm,nn,idx,PETSC_OWN_POINTER,&dd->ltogmap);CHKERRQ(ierr); 326 ierr = ISLocalToGlobalMappingBlock(dd->ltogmap,dd->w,&dd->ltogmapb);CHKERRQ(ierr); 327 ierr = PetscLogObjectParent(da,dd->ltogmap);CHKERRQ(ierr); 328 329 dd->idx = idx; 330 dd->Nl = nn; 331 332 PetscFunctionReturn(0); 333 } 334 EXTERN_C_END 335 336 #undef __FUNCT__ 337 #define __FUNCT__ "DACreate1d" 338 /*@C 339 DACreate1d - Creates an object that will manage the communication of one-dimensional 340 regular array data that is distributed across some processors. 341 342 Collective on MPI_Comm 343 344 Input Parameters: 345 + comm - MPI communicator 346 . wrap - type of periodicity should the array have, if any. Use 347 either DA_NONPERIODIC or DA_XPERIODIC 348 . M - global dimension of the array (use -M to indicate that it may be set to a different value 349 from the command line with -da_grid_x <M>) 350 . dof - number of degrees of freedom per node 351 . s - stencil width 352 - lx - array containing number of nodes in the X direction on each processor, 353 or PETSC_NULL. If non-null, must be of length as m. 354 355 Output Parameter: 356 . da - the resulting distributed array object 357 358 Options Database Key: 359 + -da_view - Calls DAView() at the conclusion of DACreate1d() 360 . -da_grid_x <nx> - number of grid points in x direction; can set if M < 0 361 - -da_refine_x - refinement factor 362 363 Level: beginner 364 365 Notes: 366 The array data itself is NOT stored in the DA, it is stored in Vec objects; 367 The appropriate vector objects can be obtained with calls to DACreateGlobalVector() 368 and DACreateLocalVector() and calls to VecDuplicate() if more are needed. 369 370 371 .keywords: distributed array, create, one-dimensional 372 373 .seealso: DADestroy(), DAView(), DACreate2d(), DACreate3d(), DAGlobalToLocalBegin(), DASetRefinementFactor(), 374 DAGlobalToLocalEnd(), DALocalToGlobal(), DALocalToLocalBegin(), DALocalToLocalEnd(), DAGetRefinementFactor(), 375 DAGetInfo(), DACreateGlobalVector(), DACreateLocalVector(), DACreateNaturalVector(), DALoad(), DAView(), DAGetOwnershipRanges() 376 377 @*/ 378 PetscErrorCode PETSCDM_DLLEXPORT DACreate1d(MPI_Comm comm, DAPeriodicType wrap, PetscInt M, PetscInt dof, PetscInt s, const PetscInt lx[], DA *da) 379 { 380 PetscErrorCode ierr; 381 PetscMPIInt size; 382 383 PetscFunctionBegin; 384 ierr = DACreate(comm, da);CHKERRQ(ierr); 385 ierr = DASetDim(*da, 1);CHKERRQ(ierr); 386 ierr = DASetSizes(*da, M, 1, 1);CHKERRQ(ierr); 387 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 388 ierr = DASetNumProcs(*da, size, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr); 389 ierr = DASetPeriodicity(*da, wrap);CHKERRQ(ierr); 390 ierr = DASetDof(*da, dof);CHKERRQ(ierr); 391 ierr = DASetStencilWidth(*da, s);CHKERRQ(ierr); 392 ierr = DASetOwnershipRanges(*da, lx, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 393 /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */ 394 ierr = DASetFromOptions(*da);CHKERRQ(ierr); 395 ierr = DASetType(*da, DA1D);CHKERRQ(ierr); 396 PetscFunctionReturn(0); 397 } 398