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 = DMDAGetCoordinates(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 = DMDAGetCoordinates(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 } else { 221 ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr); 222 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 223 } 224 225 /* 226 Determine the min and max color in plot 227 */ 228 ierr = VecStrideMin(xin,zctx.k,PETSC_NULL,&zctx.min);CHKERRQ(ierr); 229 ierr = VecStrideMax(xin,zctx.k,PETSC_NULL,&zctx.max);CHKERRQ(ierr); 230 if (zctx.k < nbounds) { 231 zctx.min = PetscMin(zctx.min,bounds[2*zctx.k]); 232 zctx.max = PetscMax(zctx.max,bounds[2*zctx.k+1]); 233 } 234 if (zctx.min == zctx.max) { 235 zctx.min -= 1.e-12; 236 zctx.max += 1.e-12; 237 } 238 239 if (!rank) { 240 const char *title; 241 242 ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr); 243 if (title) { 244 ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr); 245 } 246 } 247 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 248 ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr); 249 250 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 251 if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);} 252 253 zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min); 254 255 ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr); 256 } 257 ierr = PetscFree(displayfields);CHKERRQ(ierr); 258 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 259 260 ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 261 ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 266 #if defined(PETSC_HAVE_HDF5) 267 #undef __FUNCT__ 268 #define __FUNCT__ "VecView_MPI_HDF5_DA" 269 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 270 { 271 DM dm; 272 DM_DA *da; 273 hid_t filespace; /* file dataspace identifier */ 274 hid_t chunkspace; /* chunk dataset property identifier */ 275 hid_t plist_id; /* property list identifier */ 276 hid_t dset_id; /* dataset identifier */ 277 hid_t memspace; /* memory dataspace identifier */ 278 hid_t file_id; 279 hid_t group; 280 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 281 herr_t status; 282 hsize_t i, dim; 283 hsize_t maxDims[6], dims[6], chunkDims[6], count[6], offset[6]; 284 PetscInt timestep; 285 PetscScalar *x; 286 const char *vecname; 287 PetscErrorCode ierr; 288 289 PetscFunctionBegin; 290 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 291 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 292 293 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&dm);CHKERRQ(ierr); 294 if (!dm) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 295 da = (DM_DA*)dm->data; 296 297 /* Create the dataspace for the dataset. 298 * 299 * dims - holds the current dimensions of the dataset 300 * 301 * maxDims - holds the maximum dimensions of the dataset (unlimited 302 * for the number of time steps with the current dimensions for the 303 * other dimensions; so only additional time steps can be added). 304 * 305 * chunkDims - holds the size of a single time step (required to 306 * permit extending dataset). 307 */ 308 dim = 0; 309 if (timestep >= 0) { 310 dims[dim] = timestep+1; 311 maxDims[dim] = H5S_UNLIMITED; 312 chunkDims[dim] = 1; 313 ++dim; 314 } 315 if (da->dim == 3) { 316 dims[dim] = PetscHDF5IntCast(da->P); 317 maxDims[dim] = dims[dim]; 318 chunkDims[dim] = dims[dim]; 319 ++dim; 320 } 321 if (da->dim > 1) { 322 dims[dim] = PetscHDF5IntCast(da->N); 323 maxDims[dim] = dims[dim]; 324 chunkDims[dim] = dims[dim]; 325 ++dim; 326 } 327 dims[dim] = PetscHDF5IntCast(da->M); 328 maxDims[dim] = dims[dim]; 329 chunkDims[dim] = dims[dim]; 330 ++dim; 331 if (da->w > 1) { 332 dims[dim] = PetscHDF5IntCast(da->w); 333 maxDims[dim] = dims[dim]; 334 chunkDims[dim] = dims[dim]; 335 ++dim; 336 } 337 #if defined(PETSC_USE_COMPLEX) 338 dims[dim] = 2; 339 maxDims[dim] = dims[dim]; 340 chunkDims[dim] = dims[dim]; 341 ++dim; 342 #endif 343 for (i=0; i < dim; ++i) 344 filespace = H5Screate_simple(dim, dims, maxDims); 345 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 346 347 #if defined(PETSC_USE_REAL_SINGLE) 348 scalartype = H5T_NATIVE_FLOAT; 349 #elif defined(PETSC_USE_REAL___FLOAT128) 350 #error "HDF5 output with 128 bit floats not supported." 351 #else 352 scalartype = H5T_NATIVE_DOUBLE; 353 #endif 354 355 /* Create the dataset with default properties and close filespace */ 356 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 357 if (!H5Lexists(group, vecname, H5P_DEFAULT)) { 358 /* Create chunk */ 359 chunkspace = H5Pcreate(H5P_DATASET_CREATE); 360 if (chunkspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 361 status = H5Pset_chunk(chunkspace, dim, chunkDims); CHKERRQ(status); 362 363 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 364 dset_id = H5Dcreate2(group, vecname, scalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT); 365 #else 366 dset_id = H5Dcreate(group, vecname, scalartype, filespace, H5P_DEFAULT); 367 #endif 368 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 369 } else { 370 dset_id = H5Dopen2(group, vecname, H5P_DEFAULT); 371 status = H5Dset_extent(dset_id, dims);CHKERRQ(status); 372 } 373 status = H5Sclose(filespace);CHKERRQ(status); 374 375 /* Each process defines a dataset and writes it to the hyperslab in the file */ 376 dim = 0; 377 if (timestep >= 0) { 378 offset[dim] = timestep; 379 ++dim; 380 } 381 if (da->dim == 3) offset[dim++] = PetscHDF5IntCast(da->zs); 382 if (da->dim > 1) offset[dim++] = PetscHDF5IntCast(da->ys); 383 offset[dim++] = PetscHDF5IntCast(da->xs/da->w); 384 if (da->w > 1) offset[dim++] = 0; 385 #if defined(PETSC_USE_COMPLEX) 386 offset[dim++] = 0; 387 #endif 388 dim = 0; 389 if (timestep >= 0) { 390 count[dim] = 1; 391 ++dim; 392 } 393 if (da->dim == 3) count[dim++] = PetscHDF5IntCast(da->ze - da->zs); 394 if (da->dim > 1) count[dim++] = PetscHDF5IntCast(da->ye - da->ys); 395 count[dim++] = PetscHDF5IntCast((da->xe - da->xs)/da->w); 396 if (da->w > 1) count[dim++] = PetscHDF5IntCast(da->w); 397 #if defined(PETSC_USE_COMPLEX) 398 count[dim++] = 2; 399 #endif 400 memspace = H5Screate_simple(dim, count, NULL); 401 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 402 403 filespace = H5Dget_space(dset_id); 404 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 405 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 406 407 /* Create property list for collective dataset write */ 408 plist_id = H5Pcreate(H5P_DATASET_XFER); 409 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 410 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 411 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 412 #endif 413 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 414 415 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 416 status = H5Dwrite(dset_id, scalartype, memspace, filespace, plist_id, x);CHKERRQ(status); 417 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 418 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 419 420 /* Close/release resources */ 421 if (group != file_id) { 422 status = H5Gclose(group);CHKERRQ(status); 423 } 424 status = H5Pclose(plist_id);CHKERRQ(status); 425 status = H5Sclose(filespace);CHKERRQ(status); 426 status = H5Sclose(memspace);CHKERRQ(status); 427 status = H5Dclose(dset_id);CHKERRQ(status); 428 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 429 PetscFunctionReturn(0); 430 } 431 #endif 432 433 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 434 435 #if defined(PETSC_HAVE_MPIIO) 436 #undef __FUNCT__ 437 #define __FUNCT__ "DMDAArrayMPIIO" 438 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 439 { 440 PetscErrorCode ierr; 441 MPI_File mfdes; 442 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 443 MPI_Datatype view; 444 const PetscScalar *array; 445 MPI_Offset off; 446 MPI_Aint ub,ul; 447 PetscInt type,rows,vecrows,tr[2]; 448 DM_DA *dd = (DM_DA*)da->data; 449 450 PetscFunctionBegin; 451 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 452 if (!write) { 453 /* Read vector header. */ 454 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 455 type = tr[0]; 456 rows = tr[1]; 457 if (type != VEC_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not vector next in file"); 458 if (rows != vecrows) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 459 } else { 460 tr[0] = VEC_FILE_CLASSID; 461 tr[1] = vecrows; 462 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 463 } 464 465 dof = PetscMPIIntCast(dd->w); 466 gsizes[0] = dof; gsizes[1] = PetscMPIIntCast(dd->M); gsizes[2] = PetscMPIIntCast(dd->N); gsizes[3] = PetscMPIIntCast(dd->P); 467 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); 468 lstarts[0] = 0; lstarts[1] = PetscMPIIntCast(dd->xs/dof); lstarts[2] = PetscMPIIntCast(dd->ys); lstarts[3] = PetscMPIIntCast(dd->zs); 469 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 470 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 471 472 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 473 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 474 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char *)"native",MPI_INFO_NULL);CHKERRQ(ierr); 475 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 476 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 477 if (write) { 478 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 479 } else { 480 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 481 } 482 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 483 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 484 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 485 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 #endif 489 490 EXTERN_C_BEGIN 491 #undef __FUNCT__ 492 #define __FUNCT__ "VecView_MPI_DA" 493 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 494 { 495 DM da; 496 PetscErrorCode ierr; 497 PetscInt dim; 498 Vec natural; 499 PetscBool isdraw,isvtk; 500 #if defined(PETSC_HAVE_HDF5) 501 PetscBool ishdf5; 502 #endif 503 const char *prefix,*name; 504 505 PetscFunctionBegin; 506 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 507 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 508 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 509 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 510 #if defined(PETSC_HAVE_HDF5) 511 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 512 #endif 513 if (isdraw) { 514 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 515 if (dim == 1) { 516 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 517 } else if (dim == 2) { 518 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 519 } else { 520 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 521 } 522 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 523 Vec Y; 524 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 525 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 526 if (((PetscObject)xin)->name) { 527 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 528 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 529 } 530 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 531 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 532 #if defined(PETSC_HAVE_HDF5) 533 } else if (ishdf5) { 534 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 535 #endif 536 } else { 537 #if defined(PETSC_HAVE_MPIIO) 538 PetscBool isbinary,isMPIIO; 539 540 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 541 if (isbinary) { 542 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 543 if (isMPIIO) { 544 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 545 PetscFunctionReturn(0); 546 } 547 } 548 #endif 549 550 /* call viewer on natural ordering */ 551 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 552 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 553 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 554 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 555 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 556 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 557 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 558 ierr = VecView(natural,viewer);CHKERRQ(ierr); 559 ierr = VecDestroy(&natural);CHKERRQ(ierr); 560 } 561 PetscFunctionReturn(0); 562 } 563 EXTERN_C_END 564 565 #if defined(PETSC_HAVE_HDF5) 566 #undef __FUNCT__ 567 #define __FUNCT__ "VecLoad_HDF5_DA" 568 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 569 { 570 DM da; 571 PetscErrorCode ierr; 572 hsize_t dim; 573 hsize_t count[5]; 574 hsize_t offset[5]; 575 PetscInt cnt = 0; 576 PetscScalar *x; 577 const char *vecname; 578 hid_t filespace; /* file dataspace identifier */ 579 hid_t plist_id; /* property list identifier */ 580 hid_t dset_id; /* dataset identifier */ 581 hid_t memspace; /* memory dataspace identifier */ 582 hid_t file_id; 583 herr_t status; 584 DM_DA *dd; 585 586 PetscFunctionBegin; 587 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 588 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 589 dd = (DM_DA*)da->data; 590 591 /* Create the dataspace for the dataset */ 592 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 593 #if defined(PETSC_USE_COMPLEX) 594 dim++; 595 #endif 596 597 /* Create the dataset with default properties and close filespace */ 598 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 599 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 600 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 601 #else 602 dset_id = H5Dopen(file_id, vecname); 603 #endif 604 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 605 filespace = H5Dget_space(dset_id); 606 607 /* Each process defines a dataset and reads it from the hyperslab in the file */ 608 cnt = 0; 609 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 610 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 611 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 612 if (dd->w > 1) offset[cnt++] = 0; 613 #if defined(PETSC_USE_COMPLEX) 614 offset[cnt++] = 0; 615 #endif 616 cnt = 0; 617 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 618 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 619 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 620 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 621 #if defined(PETSC_USE_COMPLEX) 622 count[cnt++] = 2; 623 #endif 624 memspace = H5Screate_simple(dim, count, NULL); 625 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 626 627 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 628 629 /* Create property list for collective dataset write */ 630 plist_id = H5Pcreate(H5P_DATASET_XFER); 631 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 632 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 633 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 634 #endif 635 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 636 637 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 638 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 639 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 640 641 /* Close/release resources */ 642 status = H5Pclose(plist_id);CHKERRQ(status); 643 status = H5Sclose(filespace);CHKERRQ(status); 644 status = H5Sclose(memspace);CHKERRQ(status); 645 status = H5Dclose(dset_id);CHKERRQ(status); 646 PetscFunctionReturn(0); 647 } 648 #endif 649 650 #undef __FUNCT__ 651 #define __FUNCT__ "VecLoad_Binary_DA" 652 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 653 { 654 DM da; 655 PetscErrorCode ierr; 656 Vec natural; 657 const char *prefix; 658 PetscInt bs; 659 PetscBool flag; 660 DM_DA *dd; 661 #if defined(PETSC_HAVE_MPIIO) 662 PetscBool isMPIIO; 663 #endif 664 665 PetscFunctionBegin; 666 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 667 dd = (DM_DA*)da->data; 668 #if defined(PETSC_HAVE_MPIIO) 669 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 670 if (isMPIIO) { 671 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 #endif 675 676 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 677 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 678 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 679 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 680 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 681 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 682 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 683 ierr = VecDestroy(&natural);CHKERRQ(ierr); 684 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 685 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 686 if (flag && bs != dd->w) { 687 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 688 } 689 PetscFunctionReturn(0); 690 } 691 692 EXTERN_C_BEGIN 693 #undef __FUNCT__ 694 #define __FUNCT__ "VecLoad_Default_DA" 695 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 696 { 697 PetscErrorCode ierr; 698 DM da; 699 PetscBool isbinary; 700 #if defined(PETSC_HAVE_HDF5) 701 PetscBool ishdf5; 702 #endif 703 704 PetscFunctionBegin; 705 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 706 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 707 708 #if defined(PETSC_HAVE_HDF5) 709 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 710 #endif 711 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 712 713 if (isbinary) { 714 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 715 #if defined(PETSC_HAVE_HDF5) 716 } else if (ishdf5) { 717 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 718 #endif 719 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 720 PetscFunctionReturn(0); 721 } 722 EXTERN_C_END 723