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