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 103 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 104 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 105 ierr = PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);CHKERRQ(ierr); 106 107 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 108 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 109 110 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 111 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 112 113 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr); 114 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,PETSC_NULL);CHKERRQ(ierr); 115 116 /* 117 Obtain a sequential vector that is going to contain the local values plus ONE layer of 118 ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will 119 update the local values pluse ONE layer of ghost values. 120 */ 121 ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr); 122 if (!xlocal) { 123 if (bx != DMDA_BOUNDARY_NONE || by != DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 124 /* 125 if original da is not of stencil width one, or periodic or not a box stencil then 126 create a special DMDA to handle one level of ghost points for graphics 127 */ 128 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); 129 ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr); 130 } else { 131 /* otherwise we can use the da we already have */ 132 dac = da; 133 } 134 /* create local vector for holding ghosted values used in graphics */ 135 ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 136 if (dac != da) { 137 /* don't keep any public reference of this DMDA, is is only available through xlocal */ 138 ierr = PetscObjectDereference((PetscObject)dac);CHKERRQ(ierr); 139 } else { 140 /* remove association between xlocal and da, because below we compose in the opposite 141 direction and if we left this connect we'd get a loop, so the objects could 142 never be destroyed */ 143 ierr = PetscObjectRemoveReference((PetscObject)xlocal,"__PETSc_dm");CHKERRQ(ierr); 144 } 145 ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr); 146 ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr); 147 } else { 148 if (bx != DMDA_BOUNDARY_NONE || by != DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 149 ierr = VecGetDM(xlocal, &dac);CHKERRQ(ierr); 150 } else { 151 dac = da; 152 } 153 } 154 155 /* 156 Get local (ghosted) values of vector 157 */ 158 ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 159 ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 160 ierr = VecGetArray(xlocal,&zctx.v);CHKERRQ(ierr); 161 162 /* get coordinates of nodes */ 163 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 164 if (!xcoor) { 165 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 166 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 167 } 168 169 /* 170 Determine the min and max coordinates in plot 171 */ 172 ierr = VecStrideMin(xcoor,0,PETSC_NULL,&xmin);CHKERRQ(ierr); 173 ierr = VecStrideMax(xcoor,0,PETSC_NULL,&xmax);CHKERRQ(ierr); 174 ierr = VecStrideMin(xcoor,1,PETSC_NULL,&ymin);CHKERRQ(ierr); 175 ierr = VecStrideMax(xcoor,1,PETSC_NULL,&ymax);CHKERRQ(ierr); 176 coors[0] = xmin - .05*(xmax- xmin); coors[2] = xmax + .05*(xmax - xmin); 177 coors[1] = ymin - .05*(ymax- ymin); coors[3] = ymax + .05*(ymax - ymin); 178 ierr = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %G %G %G %G\n",coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 179 180 /* 181 get local ghosted version of coordinates 182 */ 183 ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr); 184 if (!xcoorl) { 185 /* create DMDA to get local version of graphics */ 186 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); 187 ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr); 188 ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr); 189 ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr); 190 ierr = PetscObjectDereference((PetscObject)dag);CHKERRQ(ierr); 191 ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr); 192 } else { 193 ierr = VecGetDM(xcoorl,&dag);CHKERRQ(ierr); 194 } 195 ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 196 ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 197 ierr = VecGetArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 198 199 /* 200 Get information about size of area each processor must do graphics for 201 */ 202 ierr = DMDAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&bx,&by,0,0);CHKERRQ(ierr); 203 ierr = DMDAGetGhostCorners(dac,0,0,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr); 204 205 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_contour_grid",&zctx.showgrid,PETSC_NULL);CHKERRQ(ierr); 206 207 ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr); 208 209 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 210 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_ports",&useports,PETSC_NULL);CHKERRQ(ierr); 211 if (useports || format == PETSC_VIEWER_DRAW_PORTS) { 212 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 213 ierr = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);CHKERRQ(ierr); 214 zctx.name0 = 0; 215 zctx.name1 = 0; 216 } else { 217 ierr = DMDAGetCoordinateName(da,0,&zctx.name0);CHKERRQ(ierr); 218 ierr = DMDAGetCoordinateName(da,1,&zctx.name1);CHKERRQ(ierr); 219 } 220 221 /* 222 Loop over each field; drawing each in a different window 223 */ 224 for (i=0; i<ndisplayfields; i++) { 225 zctx.k = displayfields[i]; 226 if (useports) { 227 ierr = PetscDrawViewPortsSet(ports,i);CHKERRQ(ierr); 228 } else { 229 ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr); 230 } 231 232 /* 233 Determine the min and max color in plot 234 */ 235 ierr = VecStrideMin(xin,zctx.k,PETSC_NULL,&zctx.min);CHKERRQ(ierr); 236 ierr = VecStrideMax(xin,zctx.k,PETSC_NULL,&zctx.max);CHKERRQ(ierr); 237 if (zctx.k < nbounds) { 238 zctx.min = bounds[2*zctx.k]; 239 zctx.max = bounds[2*zctx.k+1]; 240 } 241 if (zctx.min == zctx.max) { 242 zctx.min -= 1.e-12; 243 zctx.max += 1.e-12; 244 } 245 246 if (!rank) { 247 const char *title; 248 249 ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr); 250 if (title) { 251 ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr); 252 } 253 } 254 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 255 ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr); 256 257 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 258 if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);} 259 260 zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min); 261 262 ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr); 263 } 264 ierr = PetscFree(displayfields);CHKERRQ(ierr); 265 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 266 267 ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 268 ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 273 #if defined(PETSC_HAVE_HDF5) 274 #undef __FUNCT__ 275 #define __FUNCT__ "VecView_MPI_HDF5_DA" 276 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 277 { 278 DM dm; 279 DM_DA *da; 280 hid_t filespace; /* file dataspace identifier */ 281 hid_t chunkspace; /* chunk dataset property identifier */ 282 hid_t plist_id; /* property list identifier */ 283 hid_t dset_id; /* dataset identifier */ 284 hid_t memspace; /* memory dataspace identifier */ 285 hid_t file_id; 286 hid_t group; 287 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 288 herr_t status; 289 hsize_t i, dim; 290 hsize_t maxDims[6], dims[6], chunkDims[6], count[6], offset[6]; 291 PetscInt timestep; 292 PetscScalar *x; 293 const char *vecname; 294 PetscErrorCode ierr; 295 296 PetscFunctionBegin; 297 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 298 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 299 300 ierr = VecGetDM(xin,&dm);CHKERRQ(ierr); 301 if (!dm) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 302 da = (DM_DA*)dm->data; 303 304 /* Create the dataspace for the dataset. 305 * 306 * dims - holds the current dimensions of the dataset 307 * 308 * maxDims - holds the maximum dimensions of the dataset (unlimited 309 * for the number of time steps with the current dimensions for the 310 * other dimensions; so only additional time steps can be added). 311 * 312 * chunkDims - holds the size of a single time step (required to 313 * permit extending dataset). 314 */ 315 dim = 0; 316 if (timestep >= 0) { 317 dims[dim] = timestep+1; 318 maxDims[dim] = H5S_UNLIMITED; 319 chunkDims[dim] = 1; 320 ++dim; 321 } 322 if (da->dim == 3) { 323 ierr = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr); 324 maxDims[dim] = dims[dim]; 325 chunkDims[dim] = dims[dim]; 326 ++dim; 327 } 328 if (da->dim > 1) { 329 ierr = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr); 330 maxDims[dim] = dims[dim]; 331 chunkDims[dim] = dims[dim]; 332 ++dim; 333 } 334 ierr = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr); 335 maxDims[dim] = dims[dim]; 336 chunkDims[dim] = dims[dim]; 337 ++dim; 338 if (da->w > 1) { 339 ierr = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr); 340 maxDims[dim] = dims[dim]; 341 chunkDims[dim] = dims[dim]; 342 ++dim; 343 } 344 #if defined(PETSC_USE_COMPLEX) 345 dims[dim] = 2; 346 maxDims[dim] = dims[dim]; 347 chunkDims[dim] = dims[dim]; 348 ++dim; 349 #endif 350 for (i=0; i < dim; ++i) 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