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 { 526 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 527 } 528 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 529 Vec Y; 530 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 531 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 532 if (((PetscObject)xin)->name) { 533 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 534 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 535 } 536 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 537 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 538 #if defined(PETSC_HAVE_HDF5) 539 } else if (ishdf5) { 540 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 541 #endif 542 } else { 543 #if defined(PETSC_HAVE_MPIIO) 544 PetscBool isbinary,isMPIIO; 545 546 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 547 if (isbinary) { 548 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 549 if (isMPIIO) { 550 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 551 PetscFunctionReturn(0); 552 } 553 } 554 #endif 555 556 /* call viewer on natural ordering */ 557 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 558 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 559 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 560 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 561 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 562 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 563 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 564 ierr = VecView(natural,viewer);CHKERRQ(ierr); 565 ierr = VecDestroy(&natural);CHKERRQ(ierr); 566 } 567 PetscFunctionReturn(0); 568 } 569 EXTERN_C_END 570 571 #if defined(PETSC_HAVE_HDF5) 572 #undef __FUNCT__ 573 #define __FUNCT__ "VecLoad_HDF5_DA" 574 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 575 { 576 DM da; 577 PetscErrorCode ierr; 578 hsize_t dim; 579 hsize_t count[5]; 580 hsize_t offset[5]; 581 PetscInt cnt = 0; 582 PetscScalar *x; 583 const char *vecname; 584 hid_t filespace; /* file dataspace identifier */ 585 hid_t plist_id; /* property list identifier */ 586 hid_t dset_id; /* dataset identifier */ 587 hid_t memspace; /* memory dataspace identifier */ 588 hid_t file_id; 589 herr_t status; 590 DM_DA *dd; 591 592 PetscFunctionBegin; 593 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 594 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 595 dd = (DM_DA*)da->data; 596 597 /* Create the dataspace for the dataset */ 598 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 599 #if defined(PETSC_USE_COMPLEX) 600 dim++; 601 #endif 602 603 /* Create the dataset with default properties and close filespace */ 604 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 605 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 606 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 607 #else 608 dset_id = H5Dopen(file_id, vecname); 609 #endif 610 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 611 filespace = H5Dget_space(dset_id); 612 613 /* Each process defines a dataset and reads it from the hyperslab in the file */ 614 cnt = 0; 615 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 616 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 617 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 618 if (dd->w > 1) offset[cnt++] = 0; 619 #if defined(PETSC_USE_COMPLEX) 620 offset[cnt++] = 0; 621 #endif 622 cnt = 0; 623 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 624 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 625 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 626 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 627 #if defined(PETSC_USE_COMPLEX) 628 count[cnt++] = 2; 629 #endif 630 memspace = H5Screate_simple(dim, count, NULL); 631 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 632 633 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 634 635 /* Create property list for collective dataset write */ 636 plist_id = H5Pcreate(H5P_DATASET_XFER); 637 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 638 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 639 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 640 #endif 641 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 642 643 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 644 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 645 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 646 647 /* Close/release resources */ 648 status = H5Pclose(plist_id);CHKERRQ(status); 649 status = H5Sclose(filespace);CHKERRQ(status); 650 status = H5Sclose(memspace);CHKERRQ(status); 651 status = H5Dclose(dset_id);CHKERRQ(status); 652 PetscFunctionReturn(0); 653 } 654 #endif 655 656 #undef __FUNCT__ 657 #define __FUNCT__ "VecLoad_Binary_DA" 658 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 659 { 660 DM da; 661 PetscErrorCode ierr; 662 Vec natural; 663 const char *prefix; 664 PetscInt bs; 665 PetscBool flag; 666 DM_DA *dd; 667 #if defined(PETSC_HAVE_MPIIO) 668 PetscBool isMPIIO; 669 #endif 670 671 PetscFunctionBegin; 672 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 673 dd = (DM_DA*)da->data; 674 #if defined(PETSC_HAVE_MPIIO) 675 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 676 if (isMPIIO) { 677 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 678 PetscFunctionReturn(0); 679 } 680 #endif 681 682 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 683 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 684 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 685 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 686 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 687 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 688 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 689 ierr = VecDestroy(&natural);CHKERRQ(ierr); 690 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 691 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 692 if (flag && bs != dd->w) { 693 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 694 } 695 PetscFunctionReturn(0); 696 } 697 698 EXTERN_C_BEGIN 699 #undef __FUNCT__ 700 #define __FUNCT__ "VecLoad_Default_DA" 701 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 702 { 703 PetscErrorCode ierr; 704 DM da; 705 PetscBool isbinary; 706 #if defined(PETSC_HAVE_HDF5) 707 PetscBool ishdf5; 708 #endif 709 710 PetscFunctionBegin; 711 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 712 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 713 714 #if defined(PETSC_HAVE_HDF5) 715 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 716 #endif 717 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 718 719 if (isbinary) { 720 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 721 #if defined(PETSC_HAVE_HDF5) 722 } else if (ishdf5) { 723 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 724 #endif 725 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 726 PetscFunctionReturn(0); 727 } 728 EXTERN_C_END 729