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