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 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 (dimension == 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 (dimension == 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 hsize_t dim; 448 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 */ 449 PetscInt timestep, dimension; 450 PetscScalar *x; 451 const char *vecname; 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 456 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 457 458 ierr = VecGetDM(xin,&dm);CHKERRQ(ierr); 459 if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 460 da = (DM_DA*)dm->data; 461 ierr = DMGetDimension(dm, &dimension);CHKERRQ(ierr); 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 (dimension == 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 (dimension > 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, dimension, timestep, chunkDims); CHKERRQ(ierr); 511 512 PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims)); 513 514 #if defined(PETSC_USE_REAL_SINGLE) 515 scalartype = H5T_NATIVE_FLOAT; 516 #elif defined(PETSC_USE_REAL___FLOAT128) 517 #error "HDF5 output with 128 bit floats not supported." 518 #else 519 scalartype = H5T_NATIVE_DOUBLE; 520 #endif 521 522 /* Create the dataset with default properties and close filespace */ 523 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 524 if (!H5Lexists(group, vecname, H5P_DEFAULT)) { 525 /* Create chunk */ 526 PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE)); 527 PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims)); 528 529 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 530 PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, vecname, scalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT)); 531 #else 532 PetscStackCallHDF5Return(dset_id,H5Dcreate,(group, vecname, scalartype, filespace, H5P_DEFAULT)); 533 #endif 534 } else { 535 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 536 PetscStackCallHDF5(H5Dset_extent,(dset_id, dims)); 537 } 538 PetscStackCallHDF5(H5Sclose,(filespace)); 539 540 /* Each process defines a dataset and writes it to the hyperslab in the file */ 541 dim = 0; 542 if (timestep >= 0) { 543 offset[dim] = timestep; 544 ++dim; 545 } 546 if (dimension == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);} 547 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);} 548 ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr); 549 if (da->w > 1) offset[dim++] = 0; 550 #if defined(PETSC_USE_COMPLEX) 551 offset[dim++] = 0; 552 #endif 553 dim = 0; 554 if (timestep >= 0) { 555 count[dim] = 1; 556 ++dim; 557 } 558 if (dimension == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);} 559 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);} 560 ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr); 561 if (da->w > 1) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);} 562 #if defined(PETSC_USE_COMPLEX) 563 count[dim++] = 2; 564 #endif 565 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 566 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 567 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 568 569 /* Create property list for collective dataset write */ 570 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 571 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 572 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 573 #endif 574 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 575 576 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 577 PetscStackCallHDF5(H5Dwrite,(dset_id, scalartype, memspace, filespace, plist_id, x)); 578 PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL)); 579 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 580 581 /* Close/release resources */ 582 if (group != file_id) { 583 PetscStackCallHDF5(H5Gclose,(group)); 584 } 585 PetscStackCallHDF5(H5Pclose,(plist_id)); 586 PetscStackCallHDF5(H5Sclose,(filespace)); 587 PetscStackCallHDF5(H5Sclose,(memspace)); 588 PetscStackCallHDF5(H5Dclose,(dset_id)); 589 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 #endif 593 594 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 595 596 #if defined(PETSC_HAVE_MPIIO) 597 #undef __FUNCT__ 598 #define __FUNCT__ "DMDAArrayMPIIO" 599 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 600 { 601 PetscErrorCode ierr; 602 MPI_File mfdes; 603 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 604 MPI_Datatype view; 605 const PetscScalar *array; 606 MPI_Offset off; 607 MPI_Aint ub,ul; 608 PetscInt type,rows,vecrows,tr[2]; 609 DM_DA *dd = (DM_DA*)da->data; 610 611 PetscFunctionBegin; 612 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 613 if (!write) { 614 /* Read vector header. */ 615 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 616 type = tr[0]; 617 rows = tr[1]; 618 if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file"); 619 if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 620 } else { 621 tr[0] = VEC_FILE_CLASSID; 622 tr[1] = vecrows; 623 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 624 } 625 626 ierr = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr); 627 gsizes[0] = dof; 628 ierr = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr); 629 ierr = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr); 630 ierr = PetscMPIIntCast(dd->P,gsizes+3);CHKERRQ(ierr); 631 lsizes[0] = dof; 632 ierr = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr); 633 ierr = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr); 634 ierr = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr); 635 lstarts[0] = 0; 636 ierr = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr); 637 ierr = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr); 638 ierr = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr); 639 ierr = MPI_Type_create_subarray(da->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 640 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 641 642 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 643 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 644 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 645 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 646 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 647 if (write) { 648 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 649 } else { 650 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 651 } 652 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 653 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 654 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 655 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 #endif 659 660 #undef __FUNCT__ 661 #define __FUNCT__ "VecView_MPI_DA" 662 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 663 { 664 DM da; 665 PetscErrorCode ierr; 666 PetscInt dim; 667 Vec natural; 668 PetscBool isdraw,isvtk; 669 #if defined(PETSC_HAVE_HDF5) 670 PetscBool ishdf5; 671 #endif 672 const char *prefix,*name; 673 PetscViewerFormat format; 674 675 PetscFunctionBegin; 676 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 677 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 678 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 679 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 680 #if defined(PETSC_HAVE_HDF5) 681 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 682 #endif 683 if (isdraw) { 684 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 685 if (dim == 1) { 686 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 687 } else if (dim == 2) { 688 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 689 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 690 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 691 Vec Y; 692 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 693 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 694 if (((PetscObject)xin)->name) { 695 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 696 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 697 } 698 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 699 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 700 #if defined(PETSC_HAVE_HDF5) 701 } else if (ishdf5) { 702 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 703 #endif 704 } else { 705 #if defined(PETSC_HAVE_MPIIO) 706 PetscBool isbinary,isMPIIO; 707 708 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 709 if (isbinary) { 710 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 711 if (isMPIIO) { 712 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 713 PetscFunctionReturn(0); 714 } 715 } 716 #endif 717 718 /* call viewer on natural ordering */ 719 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 720 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 721 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 722 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 723 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 724 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 725 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 726 727 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 728 if (format == PETSC_VIEWER_BINARY_MATLAB) { 729 /* temporarily remove viewer format so it won't trigger in the VecView() */ 730 ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 731 } 732 733 ierr = VecView(natural,viewer);CHKERRQ(ierr); 734 735 if (format == PETSC_VIEWER_BINARY_MATLAB) { 736 MPI_Comm comm; 737 FILE *info; 738 const char *fieldname; 739 char fieldbuf[256]; 740 PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n; 741 742 /* set the viewer format back into the viewer */ 743 ierr = PetscViewerSetFormat(viewer,format);CHKERRQ(ierr); 744 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 745 ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr); 746 ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr); 747 ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr); 748 ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr); 749 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); } 750 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); } 751 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); } 752 753 for (n=0; n<dof; n++) { 754 ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr); 755 if (!fieldname || !fieldname[0]) { 756 ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr); 757 fieldname = fieldbuf; 758 } 759 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 760 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 761 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);} 762 } 763 ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr); 764 ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr); 765 } 766 767 ierr = VecDestroy(&natural);CHKERRQ(ierr); 768 } 769 PetscFunctionReturn(0); 770 } 771 772 #if defined(PETSC_HAVE_HDF5) 773 #undef __FUNCT__ 774 #define __FUNCT__ "VecLoad_HDF5_DA" 775 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 776 { 777 DM da; 778 PetscErrorCode ierr; 779 hsize_t dim; 780 hsize_t count[5]; 781 hsize_t offset[5]; 782 PetscInt cnt = 0, dimension; 783 PetscScalar *x; 784 const char *vecname; 785 hid_t filespace; /* file dataspace identifier */ 786 hid_t plist_id; /* property list identifier */ 787 hid_t dset_id; /* dataset identifier */ 788 hid_t memspace; /* memory dataspace identifier */ 789 hid_t file_id,group; 790 DM_DA *dd; 791 792 PetscFunctionBegin; 793 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 794 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 795 dd = (DM_DA*)da->data; 796 ierr = DMGetDimension(da, &dimension);CHKERRQ(ierr); 797 798 /* Create the dataspace for the dataset */ 799 ierr = PetscHDF5IntCast(dimension + ((dd->w == 1) ? 0 : 1),&dim);CHKERRQ(ierr); 800 #if defined(PETSC_USE_COMPLEX) 801 dim++; 802 #endif 803 804 /* Create the dataset with default properties and close filespace */ 805 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 806 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 807 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 808 #else 809 PetscStackCallHDF5Return(dset_id,H5Dopen,(group, vecname)); 810 #endif 811 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 812 813 /* Each process defines a dataset and reads it from the hyperslab in the file */ 814 cnt = 0; 815 if (dimension == 3) {ierr = PetscHDF5IntCast(dd->zs,offset + cnt++);CHKERRQ(ierr);} 816 if (dimension > 1) {ierr = PetscHDF5IntCast(dd->ys,offset + cnt++);CHKERRQ(ierr);} 817 ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + cnt++);CHKERRQ(ierr); 818 if (dd->w > 1) offset[cnt++] = 0; 819 #if defined(PETSC_USE_COMPLEX) 820 offset[cnt++] = 0; 821 #endif 822 cnt = 0; 823 if (dimension == 3) {ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + cnt++);CHKERRQ(ierr);} 824 if (dimension > 1) {ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + cnt++);CHKERRQ(ierr);} 825 ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + cnt++);CHKERRQ(ierr); 826 if (dd->w > 1) {ierr = PetscHDF5IntCast(dd->w,count + cnt++);CHKERRQ(ierr);} 827 #if defined(PETSC_USE_COMPLEX) 828 count[cnt++] = 2; 829 #endif 830 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 831 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 832 833 /* Create property list for collective dataset write */ 834 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 835 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 836 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 837 #endif 838 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 839 840 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 841 PetscStackCallHDF5(H5Dread,(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x)); 842 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 843 844 /* Close/release resources */ 845 if (group != file_id) { 846 PetscStackCallHDF5(H5Gclose,(group)); 847 } 848 PetscStackCallHDF5(H5Pclose,(plist_id)); 849 PetscStackCallHDF5(H5Sclose,(filespace)); 850 PetscStackCallHDF5(H5Sclose,(memspace)); 851 PetscStackCallHDF5(H5Dclose,(dset_id)); 852 PetscFunctionReturn(0); 853 } 854 #endif 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "VecLoad_Binary_DA" 858 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 859 { 860 DM da; 861 PetscErrorCode ierr; 862 Vec natural; 863 const char *prefix; 864 PetscInt bs; 865 PetscBool flag; 866 DM_DA *dd; 867 #if defined(PETSC_HAVE_MPIIO) 868 PetscBool isMPIIO; 869 #endif 870 871 PetscFunctionBegin; 872 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 873 dd = (DM_DA*)da->data; 874 #if defined(PETSC_HAVE_MPIIO) 875 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 876 if (isMPIIO) { 877 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 878 PetscFunctionReturn(0); 879 } 880 #endif 881 882 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 883 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 884 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 885 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 886 ierr = VecLoad(natural,viewer);CHKERRQ(ierr); 887 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 888 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 889 ierr = VecDestroy(&natural);CHKERRQ(ierr); 890 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 891 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 892 if (flag && bs != dd->w) { 893 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 894 } 895 PetscFunctionReturn(0); 896 } 897 898 #undef __FUNCT__ 899 #define __FUNCT__ "VecLoad_Default_DA" 900 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 901 { 902 PetscErrorCode ierr; 903 DM da; 904 PetscBool isbinary; 905 #if defined(PETSC_HAVE_HDF5) 906 PetscBool ishdf5; 907 #endif 908 909 PetscFunctionBegin; 910 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 911 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 912 913 #if defined(PETSC_HAVE_HDF5) 914 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 915 #endif 916 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 917 918 if (isbinary) { 919 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 920 #if defined(PETSC_HAVE_HDF5) 921 } else if (ishdf5) { 922 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 923 #endif 924 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 925 PetscFunctionReturn(0); 926 } 927