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