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