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 dims[dim] = PetscHDF5IntCast(da->P); 323 maxDims[dim] = dims[dim]; 324 chunkDims[dim] = dims[dim]; 325 ++dim; 326 } 327 if (da->dim > 1) { 328 dims[dim] = PetscHDF5IntCast(da->N); 329 maxDims[dim] = dims[dim]; 330 chunkDims[dim] = dims[dim]; 331 ++dim; 332 } 333 dims[dim] = PetscHDF5IntCast(da->M); 334 maxDims[dim] = dims[dim]; 335 chunkDims[dim] = dims[dim]; 336 ++dim; 337 if (da->w > 1) { 338 dims[dim] = PetscHDF5IntCast(da->w); 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) offset[dim++] = PetscHDF5IntCast(da->zs); 388 if (da->dim > 1) offset[dim++] = PetscHDF5IntCast(da->ys); 389 offset[dim++] = PetscHDF5IntCast(da->xs/da->w); 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) count[dim++] = PetscHDF5IntCast(da->ze - da->zs); 400 if (da->dim > 1) count[dim++] = PetscHDF5IntCast(da->ye - da->ys); 401 count[dim++] = PetscHDF5IntCast((da->xe - da->xs)/da->w); 402 if (da->w > 1) count[dim++] = PetscHDF5IntCast(da->w); 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 dof = PetscMPIIntCast(dd->w); 472 gsizes[0] = dof; gsizes[1] = PetscMPIIntCast(dd->M); gsizes[2] = PetscMPIIntCast(dd->N); gsizes[3] = PetscMPIIntCast(dd->P); 473 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); 474 lstarts[0] = 0; lstarts[1] = PetscMPIIntCast(dd->xs/dof); lstarts[2] = PetscMPIIntCast(dd->ys); lstarts[3] = PetscMPIIntCast(dd->zs); 475 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 476 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 477 478 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 479 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 480 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char *)"native",MPI_INFO_NULL);CHKERRQ(ierr); 481 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 482 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 483 if (write) { 484 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 485 } else { 486 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 487 } 488 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 489 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 490 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 491 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 492 PetscFunctionReturn(0); 493 } 494 #endif 495 496 EXTERN_C_BEGIN 497 #undef __FUNCT__ 498 #define __FUNCT__ "VecView_MPI_DA" 499 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 500 { 501 DM da; 502 PetscErrorCode ierr; 503 PetscInt dim; 504 Vec natural; 505 PetscBool isdraw,isvtk; 506 #if defined(PETSC_HAVE_HDF5) 507 PetscBool ishdf5; 508 #endif 509 const char *prefix,*name; 510 511 PetscFunctionBegin; 512 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 513 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 514 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 515 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 516 #if defined(PETSC_HAVE_HDF5) 517 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 518 #endif 519 if (isdraw) { 520 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 521 if (dim == 1) { 522 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 523 } else if (dim == 2) { 524 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 525 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 526 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 527 Vec Y; 528 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 529 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 530 if (((PetscObject)xin)->name) { 531 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 532 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 533 } 534 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 535 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 536 #if defined(PETSC_HAVE_HDF5) 537 } else if (ishdf5) { 538 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 539 #endif 540 } else { 541 #if defined(PETSC_HAVE_MPIIO) 542 PetscBool isbinary,isMPIIO; 543 544 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 545 if (isbinary) { 546 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 547 if (isMPIIO) { 548 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 549 PetscFunctionReturn(0); 550 } 551 } 552 #endif 553 554 /* call viewer on natural ordering */ 555 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 556 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 557 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 558 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 559 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 560 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 561 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 562 ierr = VecView(natural,viewer);CHKERRQ(ierr); 563 ierr = VecDestroy(&natural);CHKERRQ(ierr); 564 } 565 PetscFunctionReturn(0); 566 } 567 EXTERN_C_END 568 569 #if defined(PETSC_HAVE_HDF5) 570 #undef __FUNCT__ 571 #define __FUNCT__ "VecLoad_HDF5_DA" 572 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 573 { 574 DM da; 575 PetscErrorCode ierr; 576 hsize_t dim; 577 hsize_t count[5]; 578 hsize_t offset[5]; 579 PetscInt cnt = 0; 580 PetscScalar *x; 581 const char *vecname; 582 hid_t filespace; /* file dataspace identifier */ 583 hid_t plist_id; /* property list identifier */ 584 hid_t dset_id; /* dataset identifier */ 585 hid_t memspace; /* memory dataspace identifier */ 586 hid_t file_id; 587 herr_t status; 588 DM_DA *dd; 589 590 PetscFunctionBegin; 591 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 592 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 593 dd = (DM_DA*)da->data; 594 595 /* Create the dataspace for the dataset */ 596 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 597 #if defined(PETSC_USE_COMPLEX) 598 dim++; 599 #endif 600 601 /* Create the dataset with default properties and close filespace */ 602 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 603 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 604 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 605 #else 606 dset_id = H5Dopen(file_id, vecname); 607 #endif 608 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 609 filespace = H5Dget_space(dset_id); 610 611 /* Each process defines a dataset and reads it from the hyperslab in the file */ 612 cnt = 0; 613 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 614 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 615 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 616 if (dd->w > 1) offset[cnt++] = 0; 617 #if defined(PETSC_USE_COMPLEX) 618 offset[cnt++] = 0; 619 #endif 620 cnt = 0; 621 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 622 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 623 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 624 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 625 #if defined(PETSC_USE_COMPLEX) 626 count[cnt++] = 2; 627 #endif 628 memspace = H5Screate_simple(dim, count, NULL); 629 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 630 631 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 632 633 /* Create property list for collective dataset write */ 634 plist_id = H5Pcreate(H5P_DATASET_XFER); 635 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 636 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 637 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 638 #endif 639 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 640 641 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 642 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 643 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 644 645 /* Close/release resources */ 646 status = H5Pclose(plist_id);CHKERRQ(status); 647 status = H5Sclose(filespace);CHKERRQ(status); 648 status = H5Sclose(memspace);CHKERRQ(status); 649 status = H5Dclose(dset_id);CHKERRQ(status); 650 PetscFunctionReturn(0); 651 } 652 #endif 653 654 #undef __FUNCT__ 655 #define __FUNCT__ "VecLoad_Binary_DA" 656 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 657 { 658 DM da; 659 PetscErrorCode ierr; 660 Vec natural; 661 const char *prefix; 662 PetscInt bs; 663 PetscBool flag; 664 DM_DA *dd; 665 #if defined(PETSC_HAVE_MPIIO) 666 PetscBool isMPIIO; 667 #endif 668 669 PetscFunctionBegin; 670 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 671 dd = (DM_DA*)da->data; 672 #if defined(PETSC_HAVE_MPIIO) 673 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 674 if (isMPIIO) { 675 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 #endif 679 680 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 681 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 682 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 683 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 684 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 685 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 686 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 687 ierr = VecDestroy(&natural);CHKERRQ(ierr); 688 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 689 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 690 if (flag && bs != dd->w) { 691 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 692 } 693 PetscFunctionReturn(0); 694 } 695 696 EXTERN_C_BEGIN 697 #undef __FUNCT__ 698 #define __FUNCT__ "VecLoad_Default_DA" 699 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 700 { 701 PetscErrorCode ierr; 702 DM da; 703 PetscBool isbinary; 704 #if defined(PETSC_HAVE_HDF5) 705 PetscBool ishdf5; 706 #endif 707 708 PetscFunctionBegin; 709 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 710 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 711 712 #if defined(PETSC_HAVE_HDF5) 713 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 714 #endif 715 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 716 717 if (isbinary) { 718 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 719 #if defined(PETSC_HAVE_HDF5) 720 } else if (ishdf5) { 721 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 722 #endif 723 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 724 PetscFunctionReturn(0); 725 } 726 EXTERN_C_END 727