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