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 dimension, 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 /* following line uses sizeof(PetscReal) instead of sizeof(PetscScalar) because the last dimension of chunkDims[] captures the 2* when complex numbers are being used */ 352 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); 353 354 /* 355 if size/rank > max_chunk_size, we need radical measures: even going down to 356 avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter 357 what, composed in the most efficient way possible. 358 N.B. this minimises the number of chunks, which may or may not be the optimal 359 solution. In a BG, for example, the optimal solution is probably to make # chunks = # 360 IO nodes involved, but this author has no access to a BG to figure out how to 361 reliably find the right number. And even then it may or may not be enough. 362 */ 363 if (avg_local_vec_size > max_chunk_size) { 364 /* check if we can just split local z-axis: is that enough? */ 365 zslices = (PetscInt)ceil(vec_size*1.0/(da->p*max_chunk_size))*zslices; 366 if (zslices > da->P) { 367 /* lattice is too large in xy-directions, splitting z only is not enough */ 368 zslices = da->P; 369 yslices= (PetscInt)ceil(vec_size*1.0/(zslices*da->n*max_chunk_size))*yslices; 370 if (yslices > da->N) { 371 /* lattice is too large in x-direction, splitting along z, y is not enough */ 372 yslices = da->N; 373 xslices= (PetscInt)ceil(vec_size*1.0/(zslices*yslices*da->m*max_chunk_size))*xslices; 374 } 375 } 376 dim = 0; 377 if (timestep >= 0) { 378 ++dim; 379 } 380 /* prefer to split z-axis, even down to planar slices */ 381 if (dimension == 3) { 382 chunkDims[dim++] = (hsize_t) da->P/zslices; 383 chunkDims[dim++] = (hsize_t) da->N/yslices; 384 chunkDims[dim++] = (hsize_t) da->M/xslices; 385 } else { 386 /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */ 387 chunkDims[dim++] = (hsize_t) da->N/yslices; 388 chunkDims[dim++] = (hsize_t) da->M/xslices; 389 } 390 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); 391 } else { 392 if (target_size < chunk_size) { 393 /* only change the defaults if target_size < chunk_size */ 394 dim = 0; 395 if (timestep >= 0) { 396 ++dim; 397 } 398 /* prefer to split z-axis, even down to planar slices */ 399 if (dimension == 3) { 400 /* try splitting the z-axis to core-size bits, i.e. divide chunk size by # comm_size in z-direction */ 401 if (target_size >= chunk_size/da->p) { 402 /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */ 403 chunkDims[dim] = (hsize_t) ceil(da->P*1.0/da->p); 404 } else { 405 /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 406 radical and let everyone write all they've got */ 407 chunkDims[dim++] = (hsize_t) ceil(da->P*1.0/da->p); 408 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 409 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 410 } 411 } else { 412 /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */ 413 if (target_size >= chunk_size/da->n) { 414 /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */ 415 chunkDims[dim] = (hsize_t) ceil(da->N*1.0/da->n); 416 } else { 417 /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 418 radical and let everyone write all they've got */ 419 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 420 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 421 } 422 423 } 424 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); 425 } else { 426 /* precomputed chunks are fine, we don't need to do anything */ 427 } 428 } 429 PetscFunctionReturn(0); 430 } 431 #endif 432 433 #if defined(PETSC_HAVE_HDF5) 434 #undef __FUNCT__ 435 #define __FUNCT__ "VecView_MPI_HDF5_DA" 436 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 437 { 438 DM dm; 439 DM_DA *da; 440 hid_t filespace; /* file dataspace identifier */ 441 hid_t chunkspace; /* chunk dataset property identifier */ 442 hid_t plist_id; /* property list identifier */ 443 hid_t dset_id; /* dataset identifier */ 444 hid_t memspace; /* memory dataspace identifier */ 445 hid_t file_id; 446 hid_t group; 447 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 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, dimension; 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 ierr = DMGetDimension(dm, &dimension);CHKERRQ(ierr); 463 464 /* Create the dataspace for the dataset. 465 * 466 * dims - holds the current dimensions of the dataset 467 * 468 * maxDims - holds the maximum dimensions of the dataset (unlimited 469 * for the number of time steps with the current dimensions for the 470 * other dimensions; so only additional time steps can be added). 471 * 472 * chunkDims - holds the size of a single time step (required to 473 * permit extending dataset). 474 */ 475 dim = 0; 476 if (timestep >= 0) { 477 dims[dim] = timestep+1; 478 maxDims[dim] = H5S_UNLIMITED; 479 chunkDims[dim] = 1; 480 ++dim; 481 } 482 if (dimension == 3) { 483 ierr = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr); 484 maxDims[dim] = dims[dim]; 485 chunkDims[dim] = dims[dim]; 486 ++dim; 487 } 488 if (dimension > 1) { 489 ierr = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr); 490 maxDims[dim] = dims[dim]; 491 chunkDims[dim] = dims[dim]; 492 ++dim; 493 } 494 ierr = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr); 495 maxDims[dim] = dims[dim]; 496 chunkDims[dim] = dims[dim]; 497 ++dim; 498 if (da->w > 1) { 499 ierr = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr); 500 maxDims[dim] = dims[dim]; 501 chunkDims[dim] = dims[dim]; 502 ++dim; 503 } 504 #if defined(PETSC_USE_COMPLEX) 505 dims[dim] = 2; 506 maxDims[dim] = dims[dim]; 507 chunkDims[dim] = dims[dim]; 508 ++dim; 509 #endif 510 511 ierr = VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims); CHKERRQ(ierr); 512 513 PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims)); 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 PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE)); 528 PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims)); 529 530 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 531 PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, vecname, scalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT)); 532 #else 533 PetscStackCallHDF5Return(dset_id,H5Dcreate,(group, vecname, scalartype, filespace, H5P_DEFAULT)); 534 #endif 535 } else { 536 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 537 PetscStackCallHDF5(H5Dset_extent,(dset_id, dims)); 538 } 539 PetscStackCallHDF5(H5Sclose,(filespace)); 540 541 /* Each process defines a dataset and writes it to the hyperslab in the file */ 542 dim = 0; 543 if (timestep >= 0) { 544 offset[dim] = timestep; 545 ++dim; 546 } 547 if (dimension == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);} 548 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);} 549 ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr); 550 if (da->w > 1) offset[dim++] = 0; 551 #if defined(PETSC_USE_COMPLEX) 552 offset[dim++] = 0; 553 #endif 554 dim = 0; 555 if (timestep >= 0) { 556 count[dim] = 1; 557 ++dim; 558 } 559 if (dimension == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);} 560 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);} 561 ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr); 562 if (da->w > 1) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);} 563 #if defined(PETSC_USE_COMPLEX) 564 count[dim++] = 2; 565 #endif 566 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 567 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 568 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 569 570 /* Create property list for collective dataset write */ 571 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 572 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 573 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 574 #endif 575 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 576 577 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 578 PetscStackCallHDF5(H5Dwrite,(dset_id, scalartype, memspace, filespace, plist_id, x)); 579 PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL)); 580 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 581 582 /* Close/release resources */ 583 if (group != file_id) { 584 PetscStackCallHDF5(H5Gclose,(group)); 585 } 586 PetscStackCallHDF5(H5Pclose,(plist_id)); 587 PetscStackCallHDF5(H5Sclose,(filespace)); 588 PetscStackCallHDF5(H5Sclose,(memspace)); 589 PetscStackCallHDF5(H5Dclose,(dset_id)); 590 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 591 PetscFunctionReturn(0); 592 } 593 #endif 594 595 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 596 597 #if defined(PETSC_HAVE_MPIIO) 598 #undef __FUNCT__ 599 #define __FUNCT__ "DMDAArrayMPIIO" 600 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 601 { 602 PetscErrorCode ierr; 603 MPI_File mfdes; 604 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 605 MPI_Datatype view; 606 const PetscScalar *array; 607 MPI_Offset off; 608 MPI_Aint ub,ul; 609 PetscInt type,rows,vecrows,tr[2]; 610 DM_DA *dd = (DM_DA*)da->data; 611 612 PetscFunctionBegin; 613 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 614 if (!write) { 615 /* Read vector header. */ 616 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 617 type = tr[0]; 618 rows = tr[1]; 619 if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file"); 620 if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 621 } else { 622 tr[0] = VEC_FILE_CLASSID; 623 tr[1] = vecrows; 624 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 625 } 626 627 ierr = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr); 628 gsizes[0] = dof; 629 ierr = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr); 630 ierr = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr); 631 ierr = PetscMPIIntCast(dd->P,gsizes+3);CHKERRQ(ierr); 632 lsizes[0] = dof; 633 ierr = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr); 634 ierr = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr); 635 ierr = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr); 636 lstarts[0] = 0; 637 ierr = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr); 638 ierr = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr); 639 ierr = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr); 640 ierr = MPI_Type_create_subarray(da->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 641 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 642 643 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 644 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 645 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 646 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 647 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 648 if (write) { 649 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 650 } else { 651 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 652 } 653 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 654 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 655 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 656 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 657 PetscFunctionReturn(0); 658 } 659 #endif 660 661 #undef __FUNCT__ 662 #define __FUNCT__ "VecView_MPI_DA" 663 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 664 { 665 DM da; 666 PetscErrorCode ierr; 667 PetscInt dim; 668 Vec natural; 669 PetscBool isdraw,isvtk; 670 #if defined(PETSC_HAVE_HDF5) 671 PetscBool ishdf5; 672 #endif 673 const char *prefix,*name; 674 PetscViewerFormat format; 675 676 PetscFunctionBegin; 677 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 678 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 679 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 680 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 681 #if defined(PETSC_HAVE_HDF5) 682 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 683 #endif 684 if (isdraw) { 685 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 686 if (dim == 1) { 687 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 688 } else if (dim == 2) { 689 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 690 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 691 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 692 Vec Y; 693 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 694 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 695 if (((PetscObject)xin)->name) { 696 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 697 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 698 } 699 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 700 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 701 #if defined(PETSC_HAVE_HDF5) 702 } else if (ishdf5) { 703 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 704 #endif 705 } else { 706 #if defined(PETSC_HAVE_MPIIO) 707 PetscBool isbinary,isMPIIO; 708 709 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 710 if (isbinary) { 711 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 712 if (isMPIIO) { 713 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 714 PetscFunctionReturn(0); 715 } 716 } 717 #endif 718 719 /* call viewer on natural ordering */ 720 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 721 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 722 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 723 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 724 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 725 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 726 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 727 728 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 729 if (format == PETSC_VIEWER_BINARY_MATLAB) { 730 /* temporarily remove viewer format so it won't trigger in the VecView() */ 731 ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 732 } 733 734 ierr = VecView(natural,viewer);CHKERRQ(ierr); 735 736 if (format == PETSC_VIEWER_BINARY_MATLAB) { 737 MPI_Comm comm; 738 FILE *info; 739 const char *fieldname; 740 char fieldbuf[256]; 741 PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n; 742 743 /* set the viewer format back into the viewer */ 744 ierr = PetscViewerSetFormat(viewer,format);CHKERRQ(ierr); 745 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 746 ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr); 747 ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr); 748 ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr); 749 ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr); 750 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); } 751 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); } 752 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); } 753 754 for (n=0; n<dof; n++) { 755 ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr); 756 if (!fieldname || !fieldname[0]) { 757 ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr); 758 fieldname = fieldbuf; 759 } 760 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 761 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 762 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);} 763 } 764 ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr); 765 ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr); 766 } 767 768 ierr = VecDestroy(&natural);CHKERRQ(ierr); 769 } 770 PetscFunctionReturn(0); 771 } 772 773 #if defined(PETSC_HAVE_HDF5) 774 #undef __FUNCT__ 775 #define __FUNCT__ "VecLoad_HDF5_DA" 776 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 777 { 778 DM da; 779 PetscErrorCode ierr; 780 hsize_t dim; 781 hsize_t count[5]; 782 hsize_t offset[5]; 783 PetscInt cnt = 0, dimension; 784 PetscScalar *x; 785 const char *vecname; 786 hid_t filespace; /* file dataspace identifier */ 787 hid_t plist_id; /* property list identifier */ 788 hid_t dset_id; /* dataset identifier */ 789 hid_t memspace; /* memory dataspace identifier */ 790 hid_t file_id,group; 791 DM_DA *dd; 792 793 PetscFunctionBegin; 794 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 795 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 796 dd = (DM_DA*)da->data; 797 ierr = DMGetDimension(da, &dimension);CHKERRQ(ierr); 798 799 /* Create the dataspace for the dataset */ 800 ierr = PetscHDF5IntCast(dimension + ((dd->w == 1) ? 0 : 1),&dim);CHKERRQ(ierr); 801 #if defined(PETSC_USE_COMPLEX) 802 dim++; 803 #endif 804 805 /* Create the dataset with default properties and close filespace */ 806 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 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 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 813 814 /* Each process defines a dataset and reads it from the hyperslab in the file */ 815 cnt = 0; 816 if (dimension == 3) {ierr = PetscHDF5IntCast(dd->zs,offset + cnt++);CHKERRQ(ierr);} 817 if (dimension > 1) {ierr = PetscHDF5IntCast(dd->ys,offset + cnt++);CHKERRQ(ierr);} 818 ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + cnt++);CHKERRQ(ierr); 819 if (dd->w > 1) offset[cnt++] = 0; 820 #if defined(PETSC_USE_COMPLEX) 821 offset[cnt++] = 0; 822 #endif 823 cnt = 0; 824 if (dimension == 3) {ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + cnt++);CHKERRQ(ierr);} 825 if (dimension > 1) {ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + cnt++);CHKERRQ(ierr);} 826 ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + cnt++);CHKERRQ(ierr); 827 if (dd->w > 1) {ierr = PetscHDF5IntCast(dd->w,count + cnt++);CHKERRQ(ierr);} 828 #if defined(PETSC_USE_COMPLEX) 829 count[cnt++] = 2; 830 #endif 831 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 832 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 833 834 /* Create property list for collective dataset write */ 835 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 836 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 837 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 838 #endif 839 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 840 841 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 842 PetscStackCallHDF5(H5Dread,(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x)); 843 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 844 845 /* Close/release resources */ 846 if (group != file_id) { 847 PetscStackCallHDF5(H5Gclose,(group)); 848 } 849 PetscStackCallHDF5(H5Pclose,(plist_id)); 850 PetscStackCallHDF5(H5Sclose,(filespace)); 851 PetscStackCallHDF5(H5Sclose,(memspace)); 852 PetscStackCallHDF5(H5Dclose,(dset_id)); 853 PetscFunctionReturn(0); 854 } 855 #endif 856 857 #undef __FUNCT__ 858 #define __FUNCT__ "VecLoad_Binary_DA" 859 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 860 { 861 DM da; 862 PetscErrorCode ierr; 863 Vec natural; 864 const char *prefix; 865 PetscInt bs; 866 PetscBool flag; 867 DM_DA *dd; 868 #if defined(PETSC_HAVE_MPIIO) 869 PetscBool isMPIIO; 870 #endif 871 872 PetscFunctionBegin; 873 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 874 dd = (DM_DA*)da->data; 875 #if defined(PETSC_HAVE_MPIIO) 876 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 877 if (isMPIIO) { 878 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 879 PetscFunctionReturn(0); 880 } 881 #endif 882 883 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 884 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 885 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 886 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 887 ierr = VecLoad(natural,viewer);CHKERRQ(ierr); 888 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 889 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 890 ierr = VecDestroy(&natural);CHKERRQ(ierr); 891 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 892 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 893 if (flag && bs != dd->w) { 894 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 895 } 896 PetscFunctionReturn(0); 897 } 898 899 #undef __FUNCT__ 900 #define __FUNCT__ "VecLoad_Default_DA" 901 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 902 { 903 PetscErrorCode ierr; 904 DM da; 905 PetscBool isbinary; 906 #if defined(PETSC_HAVE_HDF5) 907 PetscBool ishdf5; 908 #endif 909 910 PetscFunctionBegin; 911 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 912 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 913 914 #if defined(PETSC_HAVE_HDF5) 915 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 916 #endif 917 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 918 919 if (isbinary) { 920 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 921 #if defined(PETSC_HAVE_HDF5) 922 } else if (ishdf5) { 923 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 924 #endif 925 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 926 PetscFunctionReturn(0); 927 } 928