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