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