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