1 2 /* 3 Plots vectors obtained with DMDACreate2d() 4 */ 5 6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 #include <petsc/private/vecimpl.h> 8 #include <petscdraw.h> 9 #include <petscviewerhdf5.h> 10 11 /* 12 The data that is passed into the graphics callback 13 */ 14 typedef struct { 15 PetscMPIInt rank; 16 PetscInt m,n,dof,k; 17 PetscReal xmin,xmax,ymin,ymax,min,max; 18 const PetscScalar *xy,*v; 19 PetscBool showgrid; 20 const char *name0,*name1; 21 } ZoomCtx; 22 23 /* 24 This does the drawing for one particular field 25 in one particular set of coordinates. It is a callback 26 called from PetscDrawZoom() 27 */ 28 #undef __FUNCT__ 29 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom" 30 PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx) 31 { 32 ZoomCtx *zctx = (ZoomCtx*)ctx; 33 PetscErrorCode ierr; 34 PetscInt m,n,i,j,k,dof,id,c1,c2,c3,c4; 35 PetscReal min,max,x1,x2,x3,x4,y_1,y2,y3,y4; 36 const PetscScalar *xy,*v; 37 38 PetscFunctionBegin; 39 m = zctx->m; 40 n = zctx->n; 41 dof = zctx->dof; 42 k = zctx->k; 43 xy = zctx->xy; 44 v = zctx->v; 45 min = zctx->min; 46 max = zctx->max; 47 48 /* PetscDraw the contour plot patch */ 49 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 50 for (j=0; j<n-1; j++) { 51 for (i=0; i<m-1; i++) { 52 id = i+j*m; 53 x1 = PetscRealPart(xy[2*id]); 54 y_1 = PetscRealPart(xy[2*id+1]); 55 c1 = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max); 56 57 id = i+j*m+1; 58 x2 = PetscRealPart(xy[2*id]); 59 y2 = PetscRealPart(xy[2*id+1]); 60 c2 = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max); 61 62 id = i+j*m+1+m; 63 x3 = PetscRealPart(xy[2*id]); 64 y3 = PetscRealPart(xy[2*id+1]); 65 c3 = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max); 66 67 id = i+j*m+m; 68 x4 = PetscRealPart(xy[2*id]); 69 y4 = PetscRealPart(xy[2*id+1]); 70 c4 = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max); 71 72 ierr = PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);CHKERRQ(ierr); 73 ierr = PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);CHKERRQ(ierr); 74 if (zctx->showgrid) { 75 ierr = PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);CHKERRQ(ierr); 76 ierr = PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);CHKERRQ(ierr); 77 ierr = PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);CHKERRQ(ierr); 78 ierr = PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);CHKERRQ(ierr); 79 } 80 } 81 } 82 if (!zctx->rank) { 83 if (zctx->name0 || zctx->name1) { 84 PetscReal xl,yl,xr,yr,x,y; 85 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 86 x = xl + .30*(xr - xl); 87 xl = xl + .01*(xr - xl); 88 y = yr - .30*(yr - yl); 89 yl = yl + .01*(yr - yl); 90 if (zctx->name0) {ierr = PetscDrawString(draw,x,yl,PETSC_DRAW_BLACK,zctx->name0);CHKERRQ(ierr);} 91 if (zctx->name1) {ierr = PetscDrawStringVertical(draw,xl,y,PETSC_DRAW_BLACK,zctx->name1);CHKERRQ(ierr);} 92 } 93 /* 94 Ideally we would use the PetscDrawAxis object to manage displaying the coordinate limits 95 but that may require some refactoring. 96 */ 97 { 98 PetscReal xmin = zctx->xmin, ymin = zctx->ymin; 99 PetscReal xmax = zctx->xmax, ymax = zctx->ymax; 100 char value[16]; size_t len; PetscReal w; 101 ierr = PetscSNPrintf(value,16,"%f",xmin);CHKERRQ(ierr); 102 ierr = PetscDrawString(draw,xmin,ymin - .05*(ymax - ymin),PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 103 ierr = PetscSNPrintf(value,16,"%f",xmax);CHKERRQ(ierr); 104 ierr = PetscStrlen(value,&len);CHKERRQ(ierr); 105 ierr = PetscDrawStringGetSize(draw,&w,NULL);CHKERRQ(ierr); 106 ierr = PetscDrawString(draw,xmax - len*w,ymin - .05*(ymax - ymin),PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 107 ierr = PetscSNPrintf(value,16,"%f",ymin);CHKERRQ(ierr); 108 ierr = PetscDrawString(draw,xmin - .05*(xmax - xmin),ymin,PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 109 ierr = PetscSNPrintf(value,16,"%f",ymax);CHKERRQ(ierr); 110 ierr = PetscDrawString(draw,xmin - .05*(xmax - xmin),ymax,PETSC_DRAW_BLACK,value);CHKERRQ(ierr); 111 } 112 } 113 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "VecView_MPI_Draw_DA2d" 119 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer) 120 { 121 DM da,dac,dag; 122 PetscErrorCode ierr; 123 PetscInt N,s,M,w,ncoors = 4; 124 const PetscInt *lx,*ly; 125 PetscReal coors[4]; 126 PetscDraw draw,popup; 127 PetscBool isnull,useports = PETSC_FALSE; 128 MPI_Comm comm; 129 Vec xlocal,xcoor,xcoorl; 130 DMBoundaryType bx,by; 131 DMDAStencilType st; 132 ZoomCtx zctx; 133 PetscDrawViewPorts *ports = NULL; 134 PetscViewerFormat format; 135 PetscInt *displayfields; 136 PetscInt ndisplayfields,i,nbounds; 137 const PetscReal *bounds; 138 139 PetscFunctionBegin; 140 zctx.showgrid = PETSC_FALSE; 141 142 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 143 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 144 if (isnull) PetscFunctionReturn(0); 145 146 ierr = PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);CHKERRQ(ierr); 147 148 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 149 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 150 151 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 152 ierr = MPI_Comm_rank(comm,&zctx.rank);CHKERRQ(ierr); 153 154 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr); 155 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,NULL);CHKERRQ(ierr); 156 157 /* 158 Obtain a sequential vector that is going to contain the local values plus ONE layer of 159 ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will 160 update the local values pluse ONE layer of ghost values. 161 */ 162 ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr); 163 if (!xlocal) { 164 if (bx != DM_BOUNDARY_NONE || by != DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 165 /* 166 if original da is not of stencil width one, or periodic or not a box stencil then 167 create a special DMDA to handle one level of ghost points for graphics 168 */ 169 ierr = DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);CHKERRQ(ierr); 170 ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr); 171 } else { 172 /* otherwise we can use the da we already have */ 173 dac = da; 174 } 175 /* create local vector for holding ghosted values used in graphics */ 176 ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 177 if (dac != da) { 178 /* don't keep any public reference of this DMDA, is is only available through xlocal */ 179 ierr = PetscObjectDereference((PetscObject)dac);CHKERRQ(ierr); 180 } else { 181 /* remove association between xlocal and da, because below we compose in the opposite 182 direction and if we left this connect we'd get a loop, so the objects could 183 never be destroyed */ 184 ierr = PetscObjectRemoveReference((PetscObject)xlocal,"__PETSc_dm");CHKERRQ(ierr); 185 } 186 ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr); 187 ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr); 188 } else { 189 if (bx != DM_BOUNDARY_NONE || by != DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 190 ierr = VecGetDM(xlocal, &dac);CHKERRQ(ierr); 191 } else { 192 dac = da; 193 } 194 } 195 196 /* 197 Get local (ghosted) values of vector 198 */ 199 ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 200 ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 201 ierr = VecGetArrayRead(xlocal,&zctx.v);CHKERRQ(ierr); 202 203 /* get coordinates of nodes */ 204 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 205 if (!xcoor) { 206 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 207 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 208 } 209 210 /* 211 Determine the min and max coordinates in plot 212 */ 213 ierr = VecStrideMin(xcoor,0,NULL,&zctx.xmin);CHKERRQ(ierr); 214 ierr = VecStrideMax(xcoor,0,NULL,&zctx.xmax);CHKERRQ(ierr); 215 ierr = VecStrideMin(xcoor,1,NULL,&zctx.ymin);CHKERRQ(ierr); 216 ierr = VecStrideMax(xcoor,1,NULL,&zctx.ymax);CHKERRQ(ierr); 217 coors[0] = zctx.xmin - .05*(zctx.xmax- zctx.xmin); coors[2] = zctx.xmax + .05*(zctx.xmax - zctx.xmin); 218 coors[1] = zctx.ymin - .05*(zctx.ymax- zctx.ymin); coors[3] = zctx.ymax + .05*(zctx.ymax - zctx.ymin); 219 ierr = PetscOptionsGetRealArray(NULL,NULL,"-draw_coordinates",coors,&ncoors,NULL);CHKERRQ(ierr); 220 ierr = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %g %g %g %g\n",(double)coors[0],(double)coors[1],(double)coors[2],(double)coors[3]);CHKERRQ(ierr); 221 222 /* 223 get local ghosted version of coordinates 224 */ 225 ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr); 226 if (!xcoorl) { 227 /* create DMDA to get local version of graphics */ 228 ierr = DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);CHKERRQ(ierr); 229 ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr); 230 ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr); 231 ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr); 232 ierr = PetscObjectDereference((PetscObject)dag);CHKERRQ(ierr); 233 ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr); 234 } else { 235 ierr = VecGetDM(xcoorl,&dag);CHKERRQ(ierr); 236 } 237 ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 238 ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 239 ierr = VecGetArrayRead(xcoorl,&zctx.xy);CHKERRQ(ierr); 240 241 /* 242 Get information about size of area each processor must do graphics for 243 */ 244 ierr = DMDAGetInfo(dac,NULL,&M,&N,NULL,NULL,NULL,NULL,&zctx.dof,NULL,&bx,&by,NULL,NULL);CHKERRQ(ierr); 245 ierr = DMDAGetGhostCorners(dac,NULL,NULL,NULL,&zctx.m,&zctx.n,NULL);CHKERRQ(ierr); 246 ierr = PetscOptionsGetBool(NULL,NULL,"-draw_contour_grid",&zctx.showgrid,NULL);CHKERRQ(ierr); 247 248 ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr); 249 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 250 ierr = PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);CHKERRQ(ierr); 251 if (useports || format == PETSC_VIEWER_DRAW_PORTS) { 252 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 253 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 254 ierr = PetscDrawClear(draw);CHKERRQ(ierr); 255 ierr = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);CHKERRQ(ierr); 256 } 257 ierr = DMDAGetCoordinateName(da,0,&zctx.name0);CHKERRQ(ierr); 258 ierr = DMDAGetCoordinateName(da,1,&zctx.name1);CHKERRQ(ierr); 259 260 /* loop over each field; drawing each in a different window */ 261 for (i=0; i<ndisplayfields; i++) { 262 zctx.k = displayfields[i]; 263 if (useports || format == PETSC_VIEWER_DRAW_PORTS) { 264 ierr = PetscDrawViewPortsSet(ports,i);CHKERRQ(ierr); 265 } else { 266 const char *title; 267 ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr); 268 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 269 ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr); 270 if (title) {ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);} 271 } 272 273 /* determine the min and max value in plot */ 274 ierr = VecStrideMin(xin,zctx.k,NULL,&zctx.min);CHKERRQ(ierr); 275 ierr = VecStrideMax(xin,zctx.k,NULL,&zctx.max);CHKERRQ(ierr); 276 if (zctx.k < nbounds) { 277 zctx.min = bounds[2*zctx.k]; 278 zctx.max = bounds[2*zctx.k+1]; 279 } 280 if (zctx.min == zctx.max) { 281 zctx.min -= 1.e-12; 282 zctx.max += 1.e-12; 283 } 284 ierr = PetscInfo2(da,"DMDA 2d contour plot min %g max %g\n",(double)zctx.min,(double)zctx.max);CHKERRQ(ierr); 285 286 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 287 if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);} 288 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 289 ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr); 290 } 291 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 292 ierr = PetscFree(displayfields);CHKERRQ(ierr); 293 294 ierr = VecRestoreArrayRead(xcoorl,&zctx.xy);CHKERRQ(ierr); 295 ierr = VecRestoreArrayRead(xlocal,&zctx.v);CHKERRQ(ierr); 296 PetscFunctionReturn(0); 297 } 298 299 #if defined(PETSC_HAVE_HDF5) 300 #undef __FUNCT__ 301 #define __FUNCT__ "VecGetHDF5ChunkSize" 302 static PetscErrorCode VecGetHDF5ChunkSize(DM_DA *da, Vec xin, PetscInt dimension, PetscInt timestep, hsize_t *chunkDims) 303 { 304 PetscMPIInt comm_size; 305 PetscErrorCode ierr; 306 hsize_t chunk_size, target_size, dim; 307 hsize_t vec_size = sizeof(PetscScalar)*da->M*da->N*da->P*da->w; 308 hsize_t avg_local_vec_size,KiB = 1024,MiB = KiB*KiB,GiB = MiB*KiB,min_size = MiB; 309 hsize_t max_chunks = 64*KiB; /* HDF5 internal limitation */ 310 hsize_t max_chunk_size = 4*GiB; /* HDF5 internal limitation */ 311 PetscInt zslices=da->p, yslices=da->n, xslices=da->m; 312 313 PetscFunctionBegin; 314 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)xin), &comm_size);CHKERRQ(ierr); 315 avg_local_vec_size = (hsize_t) ceil(vec_size*1.0/comm_size); /* we will attempt to use this as the chunk size */ 316 317 target_size = (hsize_t) ceil(PetscMin(vec_size,PetscMin(max_chunk_size,PetscMax(avg_local_vec_size,PetscMax(ceil(vec_size*1.0/max_chunks),min_size))))); 318 /* following line uses sizeof(PetscReal) instead of sizeof(PetscScalar) because the last dimension of chunkDims[] captures the 2* when complex numbers are being used */ 319 chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(PetscReal); 320 321 /* 322 if size/rank > max_chunk_size, we need radical measures: even going down to 323 avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter 324 what, composed in the most efficient way possible. 325 N.B. this minimises the number of chunks, which may or may not be the optimal 326 solution. In a BG, for example, the optimal solution is probably to make # chunks = # 327 IO nodes involved, but this author has no access to a BG to figure out how to 328 reliably find the right number. And even then it may or may not be enough. 329 */ 330 if (avg_local_vec_size > max_chunk_size) { 331 /* check if we can just split local z-axis: is that enough? */ 332 zslices = (PetscInt)ceil(vec_size*1.0/(da->p*max_chunk_size))*zslices; 333 if (zslices > da->P) { 334 /* lattice is too large in xy-directions, splitting z only is not enough */ 335 zslices = da->P; 336 yslices= (PetscInt)ceil(vec_size*1.0/(zslices*da->n*max_chunk_size))*yslices; 337 if (yslices > da->N) { 338 /* lattice is too large in x-direction, splitting along z, y is not enough */ 339 yslices = da->N; 340 xslices= (PetscInt)ceil(vec_size*1.0/(zslices*yslices*da->m*max_chunk_size))*xslices; 341 } 342 } 343 dim = 0; 344 if (timestep >= 0) { 345 ++dim; 346 } 347 /* prefer to split z-axis, even down to planar slices */ 348 if (dimension == 3) { 349 chunkDims[dim++] = (hsize_t) da->P/zslices; 350 chunkDims[dim++] = (hsize_t) da->N/yslices; 351 chunkDims[dim++] = (hsize_t) da->M/xslices; 352 } else { 353 /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */ 354 chunkDims[dim++] = (hsize_t) da->N/yslices; 355 chunkDims[dim++] = (hsize_t) da->M/xslices; 356 } 357 chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(double); 358 } else { 359 if (target_size < chunk_size) { 360 /* only change the defaults if target_size < chunk_size */ 361 dim = 0; 362 if (timestep >= 0) { 363 ++dim; 364 } 365 /* prefer to split z-axis, even down to planar slices */ 366 if (dimension == 3) { 367 /* try splitting the z-axis to core-size bits, i.e. divide chunk size by # comm_size in z-direction */ 368 if (target_size >= chunk_size/da->p) { 369 /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */ 370 chunkDims[dim] = (hsize_t) ceil(da->P*1.0/da->p); 371 } else { 372 /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 373 radical and let everyone write all they've got */ 374 chunkDims[dim++] = (hsize_t) ceil(da->P*1.0/da->p); 375 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 376 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 377 } 378 } else { 379 /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */ 380 if (target_size >= chunk_size/da->n) { 381 /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */ 382 chunkDims[dim] = (hsize_t) ceil(da->N*1.0/da->n); 383 } else { 384 /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 385 radical and let everyone write all they've got */ 386 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 387 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 388 } 389 390 } 391 chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(double); 392 } else { 393 /* precomputed chunks are fine, we don't need to do anything */ 394 } 395 } 396 PetscFunctionReturn(0); 397 } 398 #endif 399 400 #if defined(PETSC_HAVE_HDF5) 401 #undef __FUNCT__ 402 #define __FUNCT__ "VecView_MPI_HDF5_DA" 403 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 404 { 405 DM dm; 406 DM_DA *da; 407 hid_t filespace; /* file dataspace identifier */ 408 hid_t chunkspace; /* chunk dataset property identifier */ 409 hid_t plist_id; /* property list identifier */ 410 hid_t dset_id; /* dataset identifier */ 411 hid_t memspace; /* memory dataspace identifier */ 412 hid_t file_id; 413 hid_t group; 414 hid_t memscalartype; /* scalar type for mem (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 415 hid_t filescalartype; /* scalar type for file (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 416 hsize_t dim; 417 hsize_t maxDims[6]={0}, dims[6]={0}, chunkDims[6]={0}, count[6]={0}, offset[6]={0}; /* we depend on these being sane later on */ 418 PetscInt timestep, dimension; 419 const PetscScalar *x; 420 const char *vecname; 421 PetscErrorCode ierr; 422 PetscBool dim2; 423 PetscBool spoutput; 424 425 PetscFunctionBegin; 426 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 427 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 428 ierr = PetscViewerHDF5GetBaseDimension2(viewer,&dim2);CHKERRQ(ierr); 429 ierr = PetscViewerHDF5GetSPOutput(viewer,&spoutput);CHKERRQ(ierr); 430 431 ierr = VecGetDM(xin,&dm);CHKERRQ(ierr); 432 if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 433 da = (DM_DA*)dm->data; 434 ierr = DMGetDimension(dm, &dimension);CHKERRQ(ierr); 435 436 /* Create the dataspace for the dataset. 437 * 438 * dims - holds the current dimensions of the dataset 439 * 440 * maxDims - holds the maximum dimensions of the dataset (unlimited 441 * for the number of time steps with the current dimensions for the 442 * other dimensions; so only additional time steps can be added). 443 * 444 * chunkDims - holds the size of a single time step (required to 445 * permit extending dataset). 446 */ 447 dim = 0; 448 if (timestep >= 0) { 449 dims[dim] = timestep+1; 450 maxDims[dim] = H5S_UNLIMITED; 451 chunkDims[dim] = 1; 452 ++dim; 453 } 454 if (dimension == 3) { 455 ierr = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr); 456 maxDims[dim] = dims[dim]; 457 chunkDims[dim] = dims[dim]; 458 ++dim; 459 } 460 if (dimension > 1) { 461 ierr = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr); 462 maxDims[dim] = dims[dim]; 463 chunkDims[dim] = dims[dim]; 464 ++dim; 465 } 466 ierr = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr); 467 maxDims[dim] = dims[dim]; 468 chunkDims[dim] = dims[dim]; 469 ++dim; 470 if (da->w > 1 || dim2) { 471 ierr = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr); 472 maxDims[dim] = dims[dim]; 473 chunkDims[dim] = dims[dim]; 474 ++dim; 475 } 476 #if defined(PETSC_USE_COMPLEX) 477 dims[dim] = 2; 478 maxDims[dim] = dims[dim]; 479 chunkDims[dim] = dims[dim]; 480 ++dim; 481 #endif 482 483 ierr = VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims); CHKERRQ(ierr); 484 485 PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims)); 486 487 #if defined(PETSC_USE_REAL_SINGLE) 488 memscalartype = H5T_NATIVE_FLOAT; 489 filescalartype = H5T_NATIVE_FLOAT; 490 #elif defined(PETSC_USE_REAL___FLOAT128) 491 #error "HDF5 output with 128 bit floats not supported." 492 #else 493 memscalartype = H5T_NATIVE_DOUBLE; 494 if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT; 495 else filescalartype = H5T_NATIVE_DOUBLE; 496 #endif 497 498 /* Create the dataset with default properties and close filespace */ 499 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 500 if (!H5Lexists(group, vecname, H5P_DEFAULT)) { 501 /* Create chunk */ 502 PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE)); 503 PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims)); 504 505 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 506 PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT)); 507 #else 508 PetscStackCallHDF5Return(dset_id,H5Dcreate,(group, vecname, filescalartype, filespace, H5P_DEFAULT)); 509 #endif 510 } else { 511 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 512 PetscStackCallHDF5(H5Dset_extent,(dset_id, dims)); 513 } 514 PetscStackCallHDF5(H5Sclose,(filespace)); 515 516 /* Each process defines a dataset and writes it to the hyperslab in the file */ 517 dim = 0; 518 if (timestep >= 0) { 519 offset[dim] = timestep; 520 ++dim; 521 } 522 if (dimension == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);} 523 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);} 524 ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr); 525 if (da->w > 1 || dim2) offset[dim++] = 0; 526 #if defined(PETSC_USE_COMPLEX) 527 offset[dim++] = 0; 528 #endif 529 dim = 0; 530 if (timestep >= 0) { 531 count[dim] = 1; 532 ++dim; 533 } 534 if (dimension == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);} 535 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);} 536 ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr); 537 if (da->w > 1 || dim2) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);} 538 #if defined(PETSC_USE_COMPLEX) 539 count[dim++] = 2; 540 #endif 541 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 542 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 543 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 544 545 /* Create property list for collective dataset write */ 546 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 547 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 548 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 549 #endif 550 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 551 552 ierr = VecGetArrayRead(xin, &x);CHKERRQ(ierr); 553 PetscStackCallHDF5(H5Dwrite,(dset_id, memscalartype, memspace, filespace, plist_id, x)); 554 PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL)); 555 ierr = VecRestoreArrayRead(xin, &x);CHKERRQ(ierr); 556 557 /* Close/release resources */ 558 if (group != file_id) { 559 PetscStackCallHDF5(H5Gclose,(group)); 560 } 561 PetscStackCallHDF5(H5Pclose,(plist_id)); 562 PetscStackCallHDF5(H5Sclose,(filespace)); 563 PetscStackCallHDF5(H5Sclose,(memspace)); 564 PetscStackCallHDF5(H5Dclose,(dset_id)); 565 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 566 PetscFunctionReturn(0); 567 } 568 #endif 569 570 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 571 572 #if defined(PETSC_HAVE_MPIIO) 573 #undef __FUNCT__ 574 #define __FUNCT__ "DMDAArrayMPIIO" 575 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 576 { 577 PetscErrorCode ierr; 578 MPI_File mfdes; 579 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 580 MPI_Datatype view; 581 const PetscScalar *array; 582 MPI_Offset off; 583 MPI_Aint ub,ul; 584 PetscInt type,rows,vecrows,tr[2]; 585 DM_DA *dd = (DM_DA*)da->data; 586 PetscBool skipheader; 587 588 PetscFunctionBegin; 589 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 590 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr); 591 if (!write) { 592 /* Read vector header. */ 593 if (!skipheader) { 594 ierr = PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);CHKERRQ(ierr); 595 type = tr[0]; 596 rows = tr[1]; 597 if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file"); 598 if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 599 } 600 } else { 601 tr[0] = VEC_FILE_CLASSID; 602 tr[1] = vecrows; 603 if (!skipheader) { 604 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 605 } 606 } 607 608 ierr = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr); 609 gsizes[0] = dof; 610 ierr = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr); 611 ierr = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr); 612 ierr = PetscMPIIntCast(dd->P,gsizes+3);CHKERRQ(ierr); 613 lsizes[0] = dof; 614 ierr = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr); 615 ierr = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr); 616 ierr = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr); 617 lstarts[0] = 0; 618 ierr = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr); 619 ierr = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr); 620 ierr = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr); 621 ierr = MPI_Type_create_subarray(da->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 622 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 623 624 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 625 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 626 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 627 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 628 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 629 if (write) { 630 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 631 } else { 632 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 633 } 634 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 635 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 636 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 637 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 638 PetscFunctionReturn(0); 639 } 640 #endif 641 642 #undef __FUNCT__ 643 #define __FUNCT__ "VecView_MPI_DA" 644 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 645 { 646 DM da; 647 PetscErrorCode ierr; 648 PetscInt dim; 649 Vec natural; 650 PetscBool isdraw,isvtk; 651 #if defined(PETSC_HAVE_HDF5) 652 PetscBool ishdf5; 653 #endif 654 const char *prefix,*name; 655 PetscViewerFormat format; 656 657 PetscFunctionBegin; 658 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 659 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 660 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 661 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 662 #if defined(PETSC_HAVE_HDF5) 663 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 664 #endif 665 if (isdraw) { 666 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 667 if (dim == 1) { 668 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 669 } else if (dim == 2) { 670 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 671 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 672 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 673 Vec Y; 674 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 675 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 676 if (((PetscObject)xin)->name) { 677 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 678 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 679 } 680 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 681 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 682 #if defined(PETSC_HAVE_HDF5) 683 } else if (ishdf5) { 684 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 685 #endif 686 } else { 687 #if defined(PETSC_HAVE_MPIIO) 688 PetscBool isbinary,isMPIIO; 689 690 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 691 if (isbinary) { 692 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 693 if (isMPIIO) { 694 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 } 698 #endif 699 700 /* call viewer on natural ordering */ 701 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 702 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 703 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 704 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 705 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 706 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 707 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 708 709 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 710 if (format == PETSC_VIEWER_BINARY_MATLAB) { 711 /* temporarily remove viewer format so it won't trigger in the VecView() */ 712 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 713 } 714 715 ierr = VecView(natural,viewer);CHKERRQ(ierr); 716 717 if (format == PETSC_VIEWER_BINARY_MATLAB) { 718 MPI_Comm comm; 719 FILE *info; 720 const char *fieldname; 721 char fieldbuf[256]; 722 PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n; 723 724 /* set the viewer format back into the viewer */ 725 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 726 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 727 ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr); 728 ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr); 729 ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr); 730 ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr); 731 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); } 732 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); } 733 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); } 734 735 for (n=0; n<dof; n++) { 736 ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr); 737 if (!fieldname || !fieldname[0]) { 738 ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr); 739 fieldname = fieldbuf; 740 } 741 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 742 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 743 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);} 744 } 745 ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr); 746 ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr); 747 } 748 749 ierr = VecDestroy(&natural);CHKERRQ(ierr); 750 } 751 PetscFunctionReturn(0); 752 } 753 754 #if defined(PETSC_HAVE_HDF5) 755 #undef __FUNCT__ 756 #define __FUNCT__ "VecLoad_HDF5_DA" 757 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 758 { 759 DM da; 760 PetscErrorCode ierr; 761 hsize_t dim,rdim; 762 hsize_t dims[6]={0},count[6]={0},offset[6]={0}; 763 PetscInt dimension,timestep,dofInd; 764 PetscScalar *x; 765 const char *vecname; 766 hid_t filespace; /* file dataspace identifier */ 767 hid_t plist_id; /* property list identifier */ 768 hid_t dset_id; /* dataset identifier */ 769 hid_t memspace; /* memory dataspace identifier */ 770 hid_t file_id,group; 771 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 772 DM_DA *dd; 773 PetscBool dim2 = PETSC_FALSE; 774 775 PetscFunctionBegin; 776 #if defined(PETSC_USE_REAL_SINGLE) 777 scalartype = H5T_NATIVE_FLOAT; 778 #elif defined(PETSC_USE_REAL___FLOAT128) 779 #error "HDF5 output with 128 bit floats not supported." 780 #else 781 scalartype = H5T_NATIVE_DOUBLE; 782 #endif 783 784 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 785 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 786 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 787 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 788 dd = (DM_DA*)da->data; 789 ierr = DMGetDimension(da, &dimension);CHKERRQ(ierr); 790 791 /* Open dataset */ 792 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 793 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 794 #else 795 PetscStackCallHDF5Return(dset_id,H5Dopen,(group, vecname)); 796 #endif 797 798 /* Retrieve the dataspace for the dataset */ 799 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 800 PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL)); 801 802 /* Expected dimension for holding the dof's */ 803 #if defined(PETSC_USE_COMPLEX) 804 dofInd = rdim-2; 805 #else 806 dofInd = rdim-1; 807 #endif 808 809 /* The expected number of dimensions, assuming basedimension2 = false */ 810 dim = dimension; 811 if (dd->w > 1) ++dim; 812 if (timestep >= 0) ++dim; 813 #if defined(PETSC_USE_COMPLEX) 814 ++dim; 815 #endif 816 817 /* In this case the input dataset have one extra, unexpected dimension. */ 818 if (rdim == dim+1) { 819 /* In this case the block size unity */ 820 if (dd->w == 1 && dims[dofInd] == 1) dim2 = PETSC_TRUE; 821 822 /* Special error message for the case where dof does not match the input file */ 823 else if (dd->w != (PetscInt) dims[dofInd]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Number of dofs in file is %D, not %D as expected",(PetscInt)dims[dofInd],dd->w); 824 825 /* Other cases where rdim != dim cannot be handled currently */ 826 } else if (rdim != dim) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file is %d, not %d as expected with dof = %D",rdim,dim,dd->w); 827 828 /* Set up the hyperslab size */ 829 dim = 0; 830 if (timestep >= 0) { 831 offset[dim] = timestep; 832 count[dim] = 1; 833 ++dim; 834 } 835 if (dimension == 3) { 836 ierr = PetscHDF5IntCast(dd->zs,offset + dim);CHKERRQ(ierr); 837 ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + dim);CHKERRQ(ierr); 838 ++dim; 839 } 840 if (dimension > 1) { 841 ierr = PetscHDF5IntCast(dd->ys,offset + dim);CHKERRQ(ierr); 842 ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + dim);CHKERRQ(ierr); 843 ++dim; 844 } 845 ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + dim);CHKERRQ(ierr); 846 ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + dim);CHKERRQ(ierr); 847 ++dim; 848 if (dd->w > 1 || dim2) { 849 offset[dim] = 0; 850 ierr = PetscHDF5IntCast(dd->w,count + dim);CHKERRQ(ierr); 851 ++dim; 852 } 853 #if defined(PETSC_USE_COMPLEX) 854 offset[dim] = 0; 855 count[dim] = 2; 856 ++dim; 857 #endif 858 859 /* Create the memory and filespace */ 860 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 861 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 862 863 /* Create property list for collective dataset write */ 864 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 865 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 866 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 867 #endif 868 /* To read dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 869 870 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 871 PetscStackCallHDF5(H5Dread,(dset_id, scalartype, memspace, filespace, plist_id, x)); 872 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 873 874 /* Close/release resources */ 875 if (group != file_id) { 876 PetscStackCallHDF5(H5Gclose,(group)); 877 } 878 PetscStackCallHDF5(H5Pclose,(plist_id)); 879 PetscStackCallHDF5(H5Sclose,(filespace)); 880 PetscStackCallHDF5(H5Sclose,(memspace)); 881 PetscStackCallHDF5(H5Dclose,(dset_id)); 882 PetscFunctionReturn(0); 883 } 884 #endif 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "VecLoad_Binary_DA" 888 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 889 { 890 DM da; 891 PetscErrorCode ierr; 892 Vec natural; 893 const char *prefix; 894 PetscInt bs; 895 PetscBool flag; 896 DM_DA *dd; 897 #if defined(PETSC_HAVE_MPIIO) 898 PetscBool isMPIIO; 899 #endif 900 901 PetscFunctionBegin; 902 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 903 dd = (DM_DA*)da->data; 904 #if defined(PETSC_HAVE_MPIIO) 905 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 906 if (isMPIIO) { 907 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 908 PetscFunctionReturn(0); 909 } 910 #endif 911 912 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 913 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 914 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 915 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 916 ierr = VecLoad(natural,viewer);CHKERRQ(ierr); 917 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 918 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 919 ierr = VecDestroy(&natural);CHKERRQ(ierr); 920 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 921 ierr = PetscOptionsGetInt(NULL,((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 922 if (flag && bs != dd->w) { 923 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 924 } 925 PetscFunctionReturn(0); 926 } 927 928 #undef __FUNCT__ 929 #define __FUNCT__ "VecLoad_Default_DA" 930 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 931 { 932 PetscErrorCode ierr; 933 DM da; 934 PetscBool isbinary; 935 #if defined(PETSC_HAVE_HDF5) 936 PetscBool ishdf5; 937 #endif 938 939 PetscFunctionBegin; 940 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 941 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 942 943 #if defined(PETSC_HAVE_HDF5) 944 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 945 #endif 946 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 947 948 if (isbinary) { 949 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 950 #if defined(PETSC_HAVE_HDF5) 951 } else if (ishdf5) { 952 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 953 #endif 954 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 955 PetscFunctionReturn(0); 956 } 957