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