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