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