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