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);CHKERRQ(ierr); 151 152 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&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);CHKERRQ(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 PetscInt timestep, dimension; 430 const PetscScalar *x; 431 const char *vecname; 432 PetscErrorCode ierr; 433 PetscBool dim2; 434 PetscBool spoutput; 435 436 PetscFunctionBegin; 437 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 438 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 439 ierr = PetscViewerHDF5GetBaseDimension2(viewer,&dim2);CHKERRQ(ierr); 440 ierr = PetscViewerHDF5GetSPOutput(viewer,&spoutput);CHKERRQ(ierr); 441 442 ierr = VecGetDM(xin,&dm);CHKERRQ(ierr); 443 if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 444 da = (DM_DA*)dm->data; 445 ierr = DMGetDimension(dm, &dimension);CHKERRQ(ierr); 446 447 /* Create the dataspace for the dataset. 448 * 449 * dims - holds the current dimensions of the dataset 450 * 451 * maxDims - holds the maximum dimensions of the dataset (unlimited 452 * for the number of time steps with the current dimensions for the 453 * other dimensions; so only additional time steps can be added). 454 * 455 * chunkDims - holds the size of a single time step (required to 456 * permit extending dataset). 457 */ 458 dim = 0; 459 if (timestep >= 0) { 460 dims[dim] = timestep+1; 461 maxDims[dim] = H5S_UNLIMITED; 462 chunkDims[dim] = 1; 463 ++dim; 464 } 465 if (dimension == 3) { 466 ierr = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr); 467 maxDims[dim] = dims[dim]; 468 chunkDims[dim] = dims[dim]; 469 ++dim; 470 } 471 if (dimension > 1) { 472 ierr = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr); 473 maxDims[dim] = dims[dim]; 474 chunkDims[dim] = dims[dim]; 475 ++dim; 476 } 477 ierr = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr); 478 maxDims[dim] = dims[dim]; 479 chunkDims[dim] = dims[dim]; 480 ++dim; 481 if (da->w > 1 || dim2) { 482 ierr = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr); 483 maxDims[dim] = dims[dim]; 484 chunkDims[dim] = dims[dim]; 485 ++dim; 486 } 487 #if defined(PETSC_USE_COMPLEX) 488 dims[dim] = 2; 489 maxDims[dim] = dims[dim]; 490 chunkDims[dim] = dims[dim]; 491 ++dim; 492 #endif 493 494 ierr = VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims); CHKERRQ(ierr); 495 496 PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims)); 497 498 #if defined(PETSC_USE_REAL_SINGLE) 499 memscalartype = H5T_NATIVE_FLOAT; 500 filescalartype = H5T_NATIVE_FLOAT; 501 #elif defined(PETSC_USE_REAL___FLOAT128) 502 #error "HDF5 output with 128 bit floats not supported." 503 #elif defined(PETSC_USE_REAL___FP16) 504 #error "HDF5 output with 16 bit floats not supported." 505 #else 506 memscalartype = H5T_NATIVE_DOUBLE; 507 if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT; 508 else filescalartype = H5T_NATIVE_DOUBLE; 509 #endif 510 511 /* Create the dataset with default properties and close filespace */ 512 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 513 if (!H5Lexists(group, vecname, H5P_DEFAULT)) { 514 /* Create chunk */ 515 PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE)); 516 PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims)); 517 518 PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT)); 519 } else { 520 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 521 PetscStackCallHDF5(H5Dset_extent,(dset_id, dims)); 522 } 523 PetscStackCallHDF5(H5Sclose,(filespace)); 524 525 /* Each process defines a dataset and writes it to the hyperslab in the file */ 526 dim = 0; 527 if (timestep >= 0) { 528 offset[dim] = timestep; 529 ++dim; 530 } 531 if (dimension == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);} 532 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);} 533 ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr); 534 if (da->w > 1 || dim2) offset[dim++] = 0; 535 #if defined(PETSC_USE_COMPLEX) 536 offset[dim++] = 0; 537 #endif 538 dim = 0; 539 if (timestep >= 0) { 540 count[dim] = 1; 541 ++dim; 542 } 543 if (dimension == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);} 544 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);} 545 ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr); 546 if (da->w > 1 || dim2) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);} 547 #if defined(PETSC_USE_COMPLEX) 548 count[dim++] = 2; 549 #endif 550 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 551 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 552 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 553 554 ierr = VecGetArrayRead(xin, &x);CHKERRQ(ierr); 555 PetscStackCallHDF5(H5Dwrite,(dset_id, memscalartype, memspace, filespace, hdf5->dxpl_id, x)); 556 PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL)); 557 ierr = VecRestoreArrayRead(xin, &x);CHKERRQ(ierr); 558 559 #if defined(PETSC_USE_COMPLEX) 560 { 561 PetscBool tru = PETSC_TRUE; 562 ierr = PetscViewerHDF5WriteObjectAttribute(viewer,(PetscObject)xin,"complex",PETSC_BOOL,&tru);CHKERRQ(ierr); 563 } 564 #endif 565 566 /* Close/release resources */ 567 if (group != file_id) { 568 PetscStackCallHDF5(H5Gclose,(group)); 569 } 570 PetscStackCallHDF5(H5Sclose,(filespace)); 571 PetscStackCallHDF5(H5Sclose,(memspace)); 572 PetscStackCallHDF5(H5Dclose,(dset_id)); 573 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 574 PetscFunctionReturn(0); 575 } 576 #endif 577 578 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 579 580 #if defined(PETSC_HAVE_MPIIO) 581 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 582 { 583 PetscErrorCode ierr; 584 MPI_File mfdes; 585 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 586 MPI_Datatype view; 587 const PetscScalar *array; 588 MPI_Offset off; 589 MPI_Aint ub,ul; 590 PetscInt type,rows,vecrows,tr[2]; 591 DM_DA *dd = (DM_DA*)da->data; 592 PetscBool skipheader; 593 594 PetscFunctionBegin; 595 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 596 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr); 597 if (!write) { 598 /* Read vector header. */ 599 if (!skipheader) { 600 ierr = PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);CHKERRQ(ierr); 601 type = tr[0]; 602 rows = tr[1]; 603 if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file"); 604 if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 605 } 606 } else { 607 tr[0] = VEC_FILE_CLASSID; 608 tr[1] = vecrows; 609 if (!skipheader) { 610 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 611 } 612 } 613 614 ierr = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr); 615 gsizes[0] = dof; 616 ierr = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr); 617 ierr = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr); 618 ierr = PetscMPIIntCast(dd->P,gsizes+3);CHKERRQ(ierr); 619 lsizes[0] = dof; 620 ierr = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr); 621 ierr = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr); 622 ierr = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr); 623 lstarts[0] = 0; 624 ierr = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr); 625 ierr = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr); 626 ierr = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr); 627 ierr = MPI_Type_create_subarray(da->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 628 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 629 630 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 631 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 632 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 633 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 634 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 635 if (write) { 636 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 637 } else { 638 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 639 } 640 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 641 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 642 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 643 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 644 PetscFunctionReturn(0); 645 } 646 #endif 647 648 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 649 { 650 DM da; 651 PetscErrorCode ierr; 652 PetscInt dim; 653 Vec natural; 654 PetscBool isdraw,isvtk,isglvis; 655 #if defined(PETSC_HAVE_HDF5) 656 PetscBool ishdf5; 657 #endif 658 const char *prefix,*name; 659 PetscViewerFormat format; 660 661 PetscFunctionBegin; 662 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 663 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 664 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 665 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 666 #if defined(PETSC_HAVE_HDF5) 667 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 668 #endif 669 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);CHKERRQ(ierr); 670 if (isdraw) { 671 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 672 if (dim == 1) { 673 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 674 } else if (dim == 2) { 675 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 676 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 677 } else if (isvtk) { /* Duplicate the Vec */ 678 Vec Y; 679 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 680 if (((PetscObject)xin)->name) { 681 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 682 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 683 } 684 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 685 { 686 PetscObject dmvtk; 687 PetscBool compatible,compatibleSet; 688 ierr = PetscViewerVTKGetDM(viewer,&dmvtk);CHKERRQ(ierr); 689 if (dmvtk) { 690 PetscValidHeaderSpecific((DM)dmvtk,DM_CLASSID,-1); 691 ierr = DMGetCompatibility(da,(DM)dmvtk,&compatible,&compatibleSet);CHKERRQ(ierr); 692 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."); 693 } 694 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,(PetscObject)Y);CHKERRQ(ierr); 695 } 696 #if defined(PETSC_HAVE_HDF5) 697 } else if (ishdf5) { 698 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 699 #endif 700 } else if (isglvis) { 701 ierr = VecView_GLVis(xin,viewer);CHKERRQ(ierr); 702 } else { 703 #if defined(PETSC_HAVE_MPIIO) 704 PetscBool isbinary,isMPIIO; 705 706 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 707 if (isbinary) { 708 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 709 if (isMPIIO) { 710 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 711 PetscFunctionReturn(0); 712 } 713 } 714 #endif 715 716 /* call viewer on natural ordering */ 717 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 718 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 719 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 720 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 721 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 722 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 723 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 724 725 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 726 if (format == PETSC_VIEWER_BINARY_MATLAB) { 727 /* temporarily remove viewer format so it won't trigger in the VecView() */ 728 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 729 } 730 731 ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE; 732 ierr = VecView(natural,viewer);CHKERRQ(ierr); 733 ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE; 734 735 if (format == PETSC_VIEWER_BINARY_MATLAB) { 736 MPI_Comm comm; 737 FILE *info; 738 const char *fieldname; 739 char fieldbuf[256]; 740 PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n; 741 742 /* set the viewer format back into the viewer */ 743 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 744 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 745 ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr); 746 ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr); 747 ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr); 748 ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr); 749 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); } 750 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); } 751 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); } 752 753 for (n=0; n<dof; n++) { 754 ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr); 755 if (!fieldname || !fieldname[0]) { 756 ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr); 757 fieldname = fieldbuf; 758 } 759 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 760 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 761 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);} 762 } 763 ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr); 764 ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr); 765 } 766 767 ierr = VecDestroy(&natural);CHKERRQ(ierr); 768 } 769 PetscFunctionReturn(0); 770 } 771 772 #if defined(PETSC_HAVE_HDF5) 773 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 774 { 775 PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data; 776 DM da; 777 PetscErrorCode ierr; 778 int dim,rdim; 779 hsize_t dims[6]={0},count[6]={0},offset[6]={0}; 780 PetscInt dimension,timestep,dofInd; 781 PetscScalar *x; 782 const char *vecname; 783 hid_t filespace; /* file dataspace identifier */ 784 hid_t dset_id; /* dataset identifier */ 785 hid_t memspace; /* memory dataspace identifier */ 786 hid_t file_id,group; 787 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 788 DM_DA *dd; 789 PetscBool dim2 = PETSC_FALSE; 790 791 PetscFunctionBegin; 792 #if defined(PETSC_USE_REAL_SINGLE) 793 scalartype = H5T_NATIVE_FLOAT; 794 #elif defined(PETSC_USE_REAL___FLOAT128) 795 #error "HDF5 output with 128 bit floats not supported." 796 #elif defined(PETSC_USE_REAL___FP16) 797 #error "HDF5 output with 16 bit floats not supported." 798 #else 799 scalartype = H5T_NATIVE_DOUBLE; 800 #endif 801 802 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 803 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 804 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 805 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 806 dd = (DM_DA*)da->data; 807 ierr = DMGetDimension(da, &dimension);CHKERRQ(ierr); 808 809 /* Open dataset */ 810 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 811 812 /* Retrieve the dataspace for the dataset */ 813 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 814 PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL)); 815 816 /* Expected dimension for holding the dof's */ 817 #if defined(PETSC_USE_COMPLEX) 818 dofInd = rdim-2; 819 #else 820 dofInd = rdim-1; 821 #endif 822 823 /* The expected number of dimensions, assuming basedimension2 = false */ 824 dim = dimension; 825 if (dd->w > 1) ++dim; 826 if (timestep >= 0) ++dim; 827 #if defined(PETSC_USE_COMPLEX) 828 ++dim; 829 #endif 830 831 /* In this case the input dataset have one extra, unexpected dimension. */ 832 if (rdim == dim+1) { 833 /* In this case the block size unity */ 834 if (dd->w == 1 && dims[dofInd] == 1) dim2 = PETSC_TRUE; 835 836 /* Special error message for the case where dof does not match the input file */ 837 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); 838 839 /* Other cases where rdim != dim cannot be handled currently */ 840 } 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); 841 842 /* Set up the hyperslab size */ 843 dim = 0; 844 if (timestep >= 0) { 845 offset[dim] = timestep; 846 count[dim] = 1; 847 ++dim; 848 } 849 if (dimension == 3) { 850 ierr = PetscHDF5IntCast(dd->zs,offset + dim);CHKERRQ(ierr); 851 ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + dim);CHKERRQ(ierr); 852 ++dim; 853 } 854 if (dimension > 1) { 855 ierr = PetscHDF5IntCast(dd->ys,offset + dim);CHKERRQ(ierr); 856 ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + dim);CHKERRQ(ierr); 857 ++dim; 858 } 859 ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + dim);CHKERRQ(ierr); 860 ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + dim);CHKERRQ(ierr); 861 ++dim; 862 if (dd->w > 1 || dim2) { 863 offset[dim] = 0; 864 ierr = PetscHDF5IntCast(dd->w,count + dim);CHKERRQ(ierr); 865 ++dim; 866 } 867 #if defined(PETSC_USE_COMPLEX) 868 offset[dim] = 0; 869 count[dim] = 2; 870 ++dim; 871 #endif 872 873 /* Create the memory and filespace */ 874 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 875 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 876 877 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 878 PetscStackCallHDF5(H5Dread,(dset_id, scalartype, memspace, filespace, hdf5->dxpl_id, x)); 879 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 880 881 /* Close/release resources */ 882 if (group != file_id) { 883 PetscStackCallHDF5(H5Gclose,(group)); 884 } 885 PetscStackCallHDF5(H5Sclose,(filespace)); 886 PetscStackCallHDF5(H5Sclose,(memspace)); 887 PetscStackCallHDF5(H5Dclose,(dset_id)); 888 PetscFunctionReturn(0); 889 } 890 #endif 891 892 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 893 { 894 DM da; 895 PetscErrorCode ierr; 896 Vec natural; 897 const char *prefix; 898 PetscInt bs; 899 PetscBool flag; 900 DM_DA *dd; 901 #if defined(PETSC_HAVE_MPIIO) 902 PetscBool isMPIIO; 903 #endif 904 905 PetscFunctionBegin; 906 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 907 dd = (DM_DA*)da->data; 908 #if defined(PETSC_HAVE_MPIIO) 909 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 910 if (isMPIIO) { 911 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 912 PetscFunctionReturn(0); 913 } 914 #endif 915 916 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 917 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 918 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 919 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 920 ierr = VecLoad(natural,viewer);CHKERRQ(ierr); 921 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 922 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 923 ierr = VecDestroy(&natural);CHKERRQ(ierr); 924 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 925 ierr = PetscOptionsGetInt(NULL,((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 926 if (flag && bs != dd->w) { 927 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 928 } 929 PetscFunctionReturn(0); 930 } 931 932 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 933 { 934 PetscErrorCode ierr; 935 DM da; 936 PetscBool isbinary; 937 #if defined(PETSC_HAVE_HDF5) 938 PetscBool ishdf5; 939 #endif 940 941 PetscFunctionBegin; 942 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 943 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 944 945 #if defined(PETSC_HAVE_HDF5) 946 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 947 #endif 948 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 949 950 if (isbinary) { 951 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 952 #if defined(PETSC_HAVE_HDF5) 953 } else if (ishdf5) { 954 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 955 #endif 956 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 957 PetscFunctionReturn(0); 958 } 959