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