1 2 /* 3 Plots vectors obtained with DMDACreate2d() 4 */ 5 6 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 7 #include <petsc-private/vecimpl.h> 8 9 /* 10 The data that is passed into the graphics callback 11 */ 12 typedef struct { 13 PetscInt m,n,step,k; 14 PetscReal min,max,scale; 15 PetscScalar *xy,*v; 16 PetscBool showgrid; 17 const char *name0,*name1; 18 } ZoomCtx; 19 20 /* 21 This does the drawing for one particular field 22 in one particular set of coordinates. It is a callback 23 called from PetscDrawZoom() 24 */ 25 #undef __FUNCT__ 26 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom" 27 PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx) 28 { 29 ZoomCtx *zctx = (ZoomCtx*)ctx; 30 PetscErrorCode ierr; 31 PetscInt m,n,i,j,k,step,id,c1,c2,c3,c4; 32 PetscReal s,min,max,x1,x2,x3,x4,y_1,y2,y3,y4,xmin = PETSC_MAX_REAL,xmax = PETSC_MIN_REAL,ymin = PETSC_MAX_REAL,ymax = PETSC_MIN_REAL; 33 PetscReal xminf,xmaxf,yminf,ymaxf,w; 34 PetscScalar *v,*xy; 35 char value[16]; 36 size_t len; 37 38 PetscFunctionBegin; 39 m = zctx->m; 40 n = zctx->n; 41 step = zctx->step; 42 k = zctx->k; 43 v = zctx->v; 44 xy = zctx->xy; 45 s = zctx->scale; 46 min = zctx->min; 47 max = zctx->max; 48 49 /* PetscDraw the contour plot patch */ 50 for (j=0; j<n-1; j++) { 51 for (i=0; i<m-1; i++) { 52 id = i+j*m; 53 x1 = PetscRealPart(xy[2*id]); 54 y_1 = PetscRealPart(xy[2*id+1]); 55 c1 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min)); 56 xmin = PetscMin(xmin,x1); 57 ymin = PetscMin(ymin,y_1); 58 xmax = PetscMax(xmax,x1); 59 ymax = PetscMax(ymax,y_1); 60 61 id = i+j*m+1; 62 x2 = PetscRealPart(xy[2*id]); 63 y2 = y_1; 64 c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min)); 65 xmin = PetscMin(xmin,x2); 66 xmax = PetscMax(xmax,x2); 67 68 id = i+j*m+1+m; 69 x3 = x2; 70 y3 = PetscRealPart(xy[2*id+1]); 71 c3 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min)); 72 ymin = PetscMin(ymin,y3); 73 ymax = PetscMax(ymax,y3); 74 75 id = i+j*m+m; 76 x4 = x1; 77 y4 = y3; 78 c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min)); 79 80 ierr = PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);CHKERRQ(ierr); 81 ierr = PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);CHKERRQ(ierr); 82 if (zctx->showgrid) { 83 ierr = PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);CHKERRQ(ierr); 84 ierr = PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);CHKERRQ(ierr); 85 ierr = PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);CHKERRQ(ierr); 86 ierr = PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);CHKERRQ(ierr); 87 } 88 } 89 } 90 if (zctx->name0) { 91 PetscReal xl,yl,xr,yr,x,y; 92 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 93 x = xl + .3*(xr - xl); 94 xl = xl + .01*(xr - xl); 95 y = yr - .3*(yr - yl); 96 yl = yl + .01*(yr - yl); 97 ierr = PetscDrawString(draw,x,yl,PETSC_DRAW_BLACK,zctx->name0);CHKERRQ(ierr); 98 ierr = PetscDrawStringVertical(draw,xl,y,PETSC_DRAW_BLACK,zctx->name1);CHKERRQ(ierr); 99 } 100 /* 101 Ideally we would use the PetscDrawAxis object to manage displaying the coordinate limits 102 but that may require some refactoring. 103 */ 104 ierr = MPI_Allreduce(&xmin,&xminf,1,MPIU_REAL,MPIU_MAX,((PetscObject)draw)->comm);CHKERRQ(ierr); 105 ierr = MPI_Allreduce(&xmax,&xmaxf,1,MPIU_REAL,MPIU_MAX,((PetscObject)draw)->comm);CHKERRQ(ierr); 106 ierr = MPI_Allreduce(&ymin,&yminf,1,MPIU_REAL,MPIU_MAX,((PetscObject)draw)->comm);CHKERRQ(ierr); 107 ierr = MPI_Allreduce(&ymax,&ymaxf,1,MPIU_REAL,MPIU_MAX,((PetscObject)draw)->comm);CHKERRQ(ierr); 108 ierr = PetscSNPrintf(value,16,"%f",xminf);CHKERRQ(ierr); 109 ierr = PetscDrawString(draw,xminf,yminf - .05*(ymaxf - yminf),PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 110 ierr = PetscSNPrintf(value,16,"%f",xmaxf);CHKERRQ(ierr); 111 ierr = PetscStrlen(value,&len);CHKERRQ(ierr); 112 ierr = PetscDrawStringGetSize(draw,&w,PETSC_NULL);CHKERRQ(ierr); 113 ierr = PetscDrawString(draw,xmaxf - len*w,yminf - .05*(ymaxf - yminf),PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 114 ierr = PetscSNPrintf(value,16,"%f",yminf);CHKERRQ(ierr); 115 ierr = PetscDrawString(draw,xminf - .05*(xmaxf - xminf),yminf,PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 116 ierr = PetscSNPrintf(value,16,"%f",ymaxf);CHKERRQ(ierr); 117 ierr = PetscDrawString(draw,xminf - .05*(xmaxf - xminf),ymaxf,PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "VecView_MPI_Draw_DA2d" 123 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer) 124 { 125 DM da,dac,dag; 126 PetscErrorCode ierr; 127 PetscMPIInt rank; 128 PetscInt N,s,M,w; 129 const PetscInt *lx,*ly; 130 PetscReal coors[4],ymin,ymax,xmin,xmax; 131 PetscDraw draw,popup; 132 PetscBool isnull,useports = PETSC_FALSE; 133 MPI_Comm comm; 134 Vec xlocal,xcoor,xcoorl; 135 DMDABoundaryType bx,by; 136 DMDAStencilType st; 137 ZoomCtx zctx; 138 PetscDrawViewPorts *ports = PETSC_NULL; 139 PetscViewerFormat format; 140 PetscInt *displayfields; 141 PetscInt ndisplayfields,i,nbounds; 142 const PetscReal *bounds; 143 144 PetscFunctionBegin; 145 zctx.showgrid = PETSC_FALSE; 146 147 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 148 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 149 ierr = PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);CHKERRQ(ierr); 150 151 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 152 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 153 154 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 155 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 156 157 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr); 158 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,PETSC_NULL);CHKERRQ(ierr); 159 160 /* 161 Obtain a sequential vector that is going to contain the local values plus ONE layer of 162 ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will 163 update the local values pluse ONE layer of ghost values. 164 */ 165 ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr); 166 if (!xlocal) { 167 if (bx != DMDA_BOUNDARY_NONE || by != DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 168 /* 169 if original da is not of stencil width one, or periodic or not a box stencil then 170 create a special DMDA to handle one level of ghost points for graphics 171 */ 172 ierr = DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);CHKERRQ(ierr); 173 ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr); 174 } else { 175 /* otherwise we can use the da we already have */ 176 dac = da; 177 } 178 /* create local vector for holding ghosted values used in graphics */ 179 ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 180 if (dac != da) { 181 /* don't keep any public reference of this DMDA, is is only available through xlocal */ 182 ierr = PetscObjectDereference((PetscObject)dac);CHKERRQ(ierr); 183 } else { 184 /* remove association between xlocal and da, because below we compose in the opposite 185 direction and if we left this connect we'd get a loop, so the objects could 186 never be destroyed */ 187 ierr = PetscObjectRemoveReference((PetscObject)xlocal,"__PETSc_dm");CHKERRQ(ierr); 188 } 189 ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr); 190 ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr); 191 } else { 192 if (bx != DMDA_BOUNDARY_NONE || by != DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 193 ierr = VecGetDM(xlocal, &dac);CHKERRQ(ierr); 194 } else { 195 dac = da; 196 } 197 } 198 199 /* 200 Get local (ghosted) values of vector 201 */ 202 ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 203 ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 204 ierr = VecGetArray(xlocal,&zctx.v);CHKERRQ(ierr); 205 206 /* get coordinates of nodes */ 207 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 208 if (!xcoor) { 209 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 210 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 211 } 212 213 /* 214 Determine the min and max coordinates in plot 215 */ 216 ierr = VecStrideMin(xcoor,0,PETSC_NULL,&xmin);CHKERRQ(ierr); 217 ierr = VecStrideMax(xcoor,0,PETSC_NULL,&xmax);CHKERRQ(ierr); 218 ierr = VecStrideMin(xcoor,1,PETSC_NULL,&ymin);CHKERRQ(ierr); 219 ierr = VecStrideMax(xcoor,1,PETSC_NULL,&ymax);CHKERRQ(ierr); 220 coors[0] = xmin - .05*(xmax- xmin); coors[2] = xmax + .05*(xmax - xmin); 221 coors[1] = ymin - .05*(ymax- ymin); coors[3] = ymax + .05*(ymax - ymin); 222 ierr = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %G %G %G %G\n",coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 223 224 /* 225 get local ghosted version of coordinates 226 */ 227 ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr); 228 if (!xcoorl) { 229 /* create DMDA to get local version of graphics */ 230 ierr = DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);CHKERRQ(ierr); 231 ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr); 232 ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr); 233 ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr); 234 ierr = PetscObjectDereference((PetscObject)dag);CHKERRQ(ierr); 235 ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr); 236 } else { 237 ierr = VecGetDM(xcoorl,&dag);CHKERRQ(ierr); 238 } 239 ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 240 ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 241 ierr = VecGetArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 242 243 /* 244 Get information about size of area each processor must do graphics for 245 */ 246 ierr = DMDAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&bx,&by,0,0);CHKERRQ(ierr); 247 ierr = DMDAGetGhostCorners(dac,0,0,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr); 248 249 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_contour_grid",&zctx.showgrid,PETSC_NULL);CHKERRQ(ierr); 250 251 ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr); 252 253 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 254 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_ports",&useports,PETSC_NULL);CHKERRQ(ierr); 255 if (useports || format == PETSC_VIEWER_DRAW_PORTS) { 256 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 257 ierr = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);CHKERRQ(ierr); 258 zctx.name0 = 0; 259 zctx.name1 = 0; 260 } else { 261 ierr = DMDAGetCoordinateName(da,0,&zctx.name0);CHKERRQ(ierr); 262 ierr = DMDAGetCoordinateName(da,1,&zctx.name1);CHKERRQ(ierr); 263 } 264 265 /* 266 Loop over each field; drawing each in a different window 267 */ 268 for (i=0; i<ndisplayfields; i++) { 269 zctx.k = displayfields[i]; 270 if (useports) { 271 ierr = PetscDrawViewPortsSet(ports,i);CHKERRQ(ierr); 272 } else { 273 ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr); 274 } 275 276 /* 277 Determine the min and max color in plot 278 */ 279 ierr = VecStrideMin(xin,zctx.k,PETSC_NULL,&zctx.min);CHKERRQ(ierr); 280 ierr = VecStrideMax(xin,zctx.k,PETSC_NULL,&zctx.max);CHKERRQ(ierr); 281 if (zctx.k < nbounds) { 282 zctx.min = bounds[2*zctx.k]; 283 zctx.max = bounds[2*zctx.k+1]; 284 } 285 if (zctx.min == zctx.max) { 286 zctx.min -= 1.e-12; 287 zctx.max += 1.e-12; 288 } 289 290 if (!rank) { 291 const char *title; 292 293 ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr); 294 if (title) { 295 ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr); 296 } 297 } 298 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 299 ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr); 300 301 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 302 if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);} 303 304 zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min); 305 306 ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr); 307 } 308 ierr = PetscFree(displayfields);CHKERRQ(ierr); 309 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 310 311 ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 312 ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 317 #if defined(PETSC_HAVE_HDF5) 318 #undef __FUNCT__ 319 #define __FUNCT__ "VecView_MPI_HDF5_DA" 320 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 321 { 322 DM dm; 323 DM_DA *da; 324 hid_t filespace; /* file dataspace identifier */ 325 hid_t chunkspace; /* chunk dataset property identifier */ 326 hid_t plist_id; /* property list identifier */ 327 hid_t dset_id; /* dataset identifier */ 328 hid_t memspace; /* memory dataspace identifier */ 329 hid_t file_id; 330 hid_t group; 331 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 332 herr_t status; 333 hsize_t i, dim; 334 hsize_t maxDims[6], dims[6], chunkDims[6], count[6], offset[6]; 335 PetscInt timestep; 336 PetscScalar *x; 337 const char *vecname; 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 342 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 343 344 ierr = VecGetDM(xin,&dm);CHKERRQ(ierr); 345 if (!dm) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 346 da = (DM_DA*)dm->data; 347 348 /* Create the dataspace for the dataset. 349 * 350 * dims - holds the current dimensions of the dataset 351 * 352 * maxDims - holds the maximum dimensions of the dataset (unlimited 353 * for the number of time steps with the current dimensions for the 354 * other dimensions; so only additional time steps can be added). 355 * 356 * chunkDims - holds the size of a single time step (required to 357 * permit extending dataset). 358 */ 359 dim = 0; 360 if (timestep >= 0) { 361 dims[dim] = timestep+1; 362 maxDims[dim] = H5S_UNLIMITED; 363 chunkDims[dim] = 1; 364 ++dim; 365 } 366 if (da->dim == 3) { 367 ierr = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr); 368 maxDims[dim] = dims[dim]; 369 chunkDims[dim] = dims[dim]; 370 ++dim; 371 } 372 if (da->dim > 1) { 373 ierr = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr); 374 maxDims[dim] = dims[dim]; 375 chunkDims[dim] = dims[dim]; 376 ++dim; 377 } 378 ierr = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr); 379 maxDims[dim] = dims[dim]; 380 chunkDims[dim] = dims[dim]; 381 ++dim; 382 if (da->w > 1) { 383 ierr = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr); 384 maxDims[dim] = dims[dim]; 385 chunkDims[dim] = dims[dim]; 386 ++dim; 387 } 388 #if defined(PETSC_USE_COMPLEX) 389 dims[dim] = 2; 390 maxDims[dim] = dims[dim]; 391 chunkDims[dim] = dims[dim]; 392 ++dim; 393 #endif 394 for (i=0; i < dim; ++i) filespace = H5Screate_simple(dim, dims, maxDims); 395 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 396 397 #if defined(PETSC_USE_REAL_SINGLE) 398 scalartype = H5T_NATIVE_FLOAT; 399 #elif defined(PETSC_USE_REAL___FLOAT128) 400 #error "HDF5 output with 128 bit floats not supported." 401 #else 402 scalartype = H5T_NATIVE_DOUBLE; 403 #endif 404 405 /* Create the dataset with default properties and close filespace */ 406 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 407 if (!H5Lexists(group, vecname, H5P_DEFAULT)) { 408 /* Create chunk */ 409 chunkspace = H5Pcreate(H5P_DATASET_CREATE); 410 if (chunkspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 411 status = H5Pset_chunk(chunkspace, dim, chunkDims);CHKERRQ(status); 412 413 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 414 dset_id = H5Dcreate2(group, vecname, scalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT); 415 #else 416 dset_id = H5Dcreate(group, vecname, scalartype, filespace, H5P_DEFAULT); 417 #endif 418 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 419 } else { 420 dset_id = H5Dopen2(group, vecname, H5P_DEFAULT); 421 status = H5Dset_extent(dset_id, dims);CHKERRQ(status); 422 } 423 status = H5Sclose(filespace);CHKERRQ(status); 424 425 /* Each process defines a dataset and writes it to the hyperslab in the file */ 426 dim = 0; 427 if (timestep >= 0) { 428 offset[dim] = timestep; 429 ++dim; 430 } 431 if (da->dim == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);} 432 if (da->dim > 1) {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);} 433 ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr); 434 if (da->w > 1) offset[dim++] = 0; 435 #if defined(PETSC_USE_COMPLEX) 436 offset[dim++] = 0; 437 #endif 438 dim = 0; 439 if (timestep >= 0) { 440 count[dim] = 1; 441 ++dim; 442 } 443 if (da->dim == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);} 444 if (da->dim > 1) {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);} 445 ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr); 446 if (da->w > 1) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);} 447 #if defined(PETSC_USE_COMPLEX) 448 count[dim++] = 2; 449 #endif 450 memspace = H5Screate_simple(dim, count, NULL); 451 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 452 453 filespace = H5Dget_space(dset_id); 454 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 455 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 456 457 /* Create property list for collective dataset write */ 458 plist_id = H5Pcreate(H5P_DATASET_XFER); 459 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 460 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 461 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 462 #endif 463 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 464 465 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 466 status = H5Dwrite(dset_id, scalartype, memspace, filespace, plist_id, x);CHKERRQ(status); 467 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 468 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 469 470 /* Close/release resources */ 471 if (group != file_id) { 472 status = H5Gclose(group);CHKERRQ(status); 473 } 474 status = H5Pclose(plist_id);CHKERRQ(status); 475 status = H5Sclose(filespace);CHKERRQ(status); 476 status = H5Sclose(memspace);CHKERRQ(status); 477 status = H5Dclose(dset_id);CHKERRQ(status); 478 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 479 PetscFunctionReturn(0); 480 } 481 #endif 482 483 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 484 485 #if defined(PETSC_HAVE_MPIIO) 486 #undef __FUNCT__ 487 #define __FUNCT__ "DMDAArrayMPIIO" 488 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 489 { 490 PetscErrorCode ierr; 491 MPI_File mfdes; 492 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 493 MPI_Datatype view; 494 const PetscScalar *array; 495 MPI_Offset off; 496 MPI_Aint ub,ul; 497 PetscInt type,rows,vecrows,tr[2]; 498 DM_DA *dd = (DM_DA*)da->data; 499 500 PetscFunctionBegin; 501 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 502 if (!write) { 503 /* Read vector header. */ 504 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 505 type = tr[0]; 506 rows = tr[1]; 507 if (type != VEC_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not vector next in file"); 508 if (rows != vecrows) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 509 } else { 510 tr[0] = VEC_FILE_CLASSID; 511 tr[1] = vecrows; 512 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 513 } 514 515 ierr = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr); 516 gsizes[0] = dof; 517 ierr = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr); 518 ierr = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr); 519 ierr = PetscMPIIntCast(dd->P,gsizes+1);CHKERRQ(ierr); 520 lsizes[0] = dof; 521 ierr = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr); 522 ierr = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr); 523 ierr = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr); 524 lstarts[0] = 0; 525 ierr = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr); 526 ierr = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr); 527 ierr = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr); 528 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 529 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 530 531 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 532 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 533 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 534 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 535 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 536 if (write) { 537 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 538 } else { 539 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 540 } 541 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 542 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 543 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 544 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 545 PetscFunctionReturn(0); 546 } 547 #endif 548 549 EXTERN_C_BEGIN 550 #undef __FUNCT__ 551 #define __FUNCT__ "VecView_MPI_DA" 552 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 553 { 554 DM da; 555 PetscErrorCode ierr; 556 PetscInt dim; 557 Vec natural; 558 PetscBool isdraw,isvtk; 559 #if defined(PETSC_HAVE_HDF5) 560 PetscBool ishdf5; 561 #endif 562 const char *prefix,*name; 563 564 PetscFunctionBegin; 565 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 566 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 567 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 568 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 569 #if defined(PETSC_HAVE_HDF5) 570 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 571 #endif 572 if (isdraw) { 573 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 574 if (dim == 1) { 575 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 576 } else if (dim == 2) { 577 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 578 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 579 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 580 Vec Y; 581 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 582 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 583 if (((PetscObject)xin)->name) { 584 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 585 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 586 } 587 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 588 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 589 #if defined(PETSC_HAVE_HDF5) 590 } else if (ishdf5) { 591 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 592 #endif 593 } else { 594 #if defined(PETSC_HAVE_MPIIO) 595 PetscBool isbinary,isMPIIO; 596 597 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 598 if (isbinary) { 599 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 600 if (isMPIIO) { 601 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 602 PetscFunctionReturn(0); 603 } 604 } 605 #endif 606 607 /* call viewer on natural ordering */ 608 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 609 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 610 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 611 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 612 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 613 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 614 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 615 ierr = VecView(natural,viewer);CHKERRQ(ierr); 616 ierr = VecDestroy(&natural);CHKERRQ(ierr); 617 } 618 PetscFunctionReturn(0); 619 } 620 EXTERN_C_END 621 622 #if defined(PETSC_HAVE_HDF5) 623 #undef __FUNCT__ 624 #define __FUNCT__ "VecLoad_HDF5_DA" 625 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 626 { 627 DM da; 628 PetscErrorCode ierr; 629 hsize_t dim; 630 hsize_t count[5]; 631 hsize_t offset[5]; 632 PetscInt cnt = 0; 633 PetscScalar *x; 634 const char *vecname; 635 hid_t filespace; /* file dataspace identifier */ 636 hid_t plist_id; /* property list identifier */ 637 hid_t dset_id; /* dataset identifier */ 638 hid_t memspace; /* memory dataspace identifier */ 639 hid_t file_id; 640 herr_t status; 641 DM_DA *dd; 642 643 PetscFunctionBegin; 644 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 645 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 646 dd = (DM_DA*)da->data; 647 648 /* Create the dataspace for the dataset */ 649 ierr = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1),&dim);CHKERRQ(ierr); 650 #if defined(PETSC_USE_COMPLEX) 651 dim++; 652 #endif 653 654 /* Create the dataset with default properties and close filespace */ 655 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 656 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 657 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 658 #else 659 dset_id = H5Dopen(file_id, vecname); 660 #endif 661 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 662 filespace = H5Dget_space(dset_id); 663 664 /* Each process defines a dataset and reads it from the hyperslab in the file */ 665 cnt = 0; 666 if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->zs,offset + cnt++);CHKERRQ(ierr);} 667 if (dd->dim > 1) {ierr = PetscHDF5IntCast(dd->ys,offset + cnt++);CHKERRQ(ierr);} 668 ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + cnt++);CHKERRQ(ierr); 669 if (dd->w > 1) offset[cnt++] = 0; 670 #if defined(PETSC_USE_COMPLEX) 671 offset[cnt++] = 0; 672 #endif 673 cnt = 0; 674 if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + cnt++);CHKERRQ(ierr);} 675 if (dd->dim > 1) {ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + cnt++);CHKERRQ(ierr);} 676 ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + cnt++);CHKERRQ(ierr); 677 if (dd->w > 1) {ierr = PetscHDF5IntCast(dd->w,count + cnt++);CHKERRQ(ierr);} 678 #if defined(PETSC_USE_COMPLEX) 679 count[cnt++] = 2; 680 #endif 681 memspace = H5Screate_simple(dim, count, NULL); 682 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 683 684 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 685 686 /* Create property list for collective dataset write */ 687 plist_id = H5Pcreate(H5P_DATASET_XFER); 688 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 689 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 690 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 691 #endif 692 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 693 694 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 695 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 696 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 697 698 /* Close/release resources */ 699 status = H5Pclose(plist_id);CHKERRQ(status); 700 status = H5Sclose(filespace);CHKERRQ(status); 701 status = H5Sclose(memspace);CHKERRQ(status); 702 status = H5Dclose(dset_id);CHKERRQ(status); 703 PetscFunctionReturn(0); 704 } 705 #endif 706 707 #undef __FUNCT__ 708 #define __FUNCT__ "VecLoad_Binary_DA" 709 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 710 { 711 DM da; 712 PetscErrorCode ierr; 713 Vec natural; 714 const char *prefix; 715 PetscInt bs; 716 PetscBool flag; 717 DM_DA *dd; 718 #if defined(PETSC_HAVE_MPIIO) 719 PetscBool isMPIIO; 720 #endif 721 722 PetscFunctionBegin; 723 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 724 dd = (DM_DA*)da->data; 725 #if defined(PETSC_HAVE_MPIIO) 726 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 727 if (isMPIIO) { 728 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 729 PetscFunctionReturn(0); 730 } 731 #endif 732 733 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 734 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 735 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 736 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 737 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 738 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 739 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 740 ierr = VecDestroy(&natural);CHKERRQ(ierr); 741 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 742 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 743 if (flag && bs != dd->w) { 744 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 745 } 746 PetscFunctionReturn(0); 747 } 748 749 EXTERN_C_BEGIN 750 #undef __FUNCT__ 751 #define __FUNCT__ "VecLoad_Default_DA" 752 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 753 { 754 PetscErrorCode ierr; 755 DM da; 756 PetscBool isbinary; 757 #if defined(PETSC_HAVE_HDF5) 758 PetscBool ishdf5; 759 #endif 760 761 PetscFunctionBegin; 762 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 763 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 764 765 #if defined(PETSC_HAVE_HDF5) 766 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 767 #endif 768 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 769 770 if (isbinary) { 771 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 772 #if defined(PETSC_HAVE_HDF5) 773 } else if (ishdf5) { 774 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 775 #endif 776 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 777 PetscFunctionReturn(0); 778 } 779 EXTERN_C_END 780