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