1 2 /* 3 Plots vectors obtained with DMDACreate2d() 4 */ 5 6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 #include <petscdraw.h> 8 #include <petscviewerhdf5.h> 9 10 /* 11 The data that is passed into the graphics callback 12 */ 13 typedef struct { 14 PetscMPIInt rank; 15 PetscInt m,n,dof,k; 16 PetscReal xmin,xmax,ymin,ymax,min,max; 17 const PetscScalar *xy,*v; 18 PetscBool showaxis,showgrid; 19 const char *name0,*name1; 20 } ZoomCtx; 21 22 /* 23 This does the drawing for one particular field 24 in one particular set of coordinates. It is a callback 25 called from PetscDrawZoom() 26 */ 27 #undef __FUNCT__ 28 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom" 29 PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx) 30 { 31 ZoomCtx *zctx = (ZoomCtx*)ctx; 32 PetscErrorCode ierr; 33 PetscInt m,n,i,j,k,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 #undef __FUNCT__ 117 #define __FUNCT__ "VecView_MPI_Draw_DA2d" 118 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer) 119 { 120 DM da,dac,dag; 121 PetscErrorCode ierr; 122 PetscInt N,s,M,w,ncoors = 4; 123 const PetscInt *lx,*ly; 124 PetscReal coors[4]; 125 PetscDraw draw,popup; 126 PetscBool isnull,useports = PETSC_FALSE; 127 MPI_Comm comm; 128 Vec xlocal,xcoor,xcoorl; 129 DMBoundaryType bx,by; 130 DMDAStencilType st; 131 ZoomCtx zctx; 132 PetscDrawViewPorts *ports = NULL; 133 PetscViewerFormat format; 134 PetscInt *displayfields; 135 PetscInt ndisplayfields,i,nbounds; 136 const PetscReal *bounds; 137 138 PetscFunctionBegin; 139 zctx.showgrid = PETSC_FALSE; 140 zctx.showaxis = PETSC_TRUE; 141 142 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 143 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 144 if (isnull) PetscFunctionReturn(0); 145 146 ierr = PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);CHKERRQ(ierr); 147 148 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 149 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 150 151 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 152 ierr = MPI_Comm_rank(comm,&zctx.rank);CHKERRQ(ierr); 153 154 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr); 155 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,NULL);CHKERRQ(ierr); 156 157 /* 158 Obtain a sequential vector that is going to contain the local values plus ONE layer of 159 ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will 160 update the local values pluse ONE layer of ghost values. 161 */ 162 ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr); 163 if (!xlocal) { 164 if (bx != DM_BOUNDARY_NONE || by != DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 165 /* 166 if original da is not of stencil width one, or periodic or not a box stencil then 167 create a special DMDA to handle one level of ghost points for graphics 168 */ 169 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); 170 ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr); 171 } else { 172 /* otherwise we can use the da we already have */ 173 dac = da; 174 } 175 /* create local vector for holding ghosted values used in graphics */ 176 ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 177 if (dac != da) { 178 /* don't keep any public reference of this DMDA, is is only available through xlocal */ 179 ierr = PetscObjectDereference((PetscObject)dac);CHKERRQ(ierr); 180 } else { 181 /* remove association between xlocal and da, because below we compose in the opposite 182 direction and if we left this connect we'd get a loop, so the objects could 183 never be destroyed */ 184 ierr = PetscObjectRemoveReference((PetscObject)xlocal,"__PETSc_dm");CHKERRQ(ierr); 185 } 186 ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr); 187 ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr); 188 } else { 189 if (bx != DM_BOUNDARY_NONE || by != DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 190 ierr = VecGetDM(xlocal, &dac);CHKERRQ(ierr); 191 } else { 192 dac = da; 193 } 194 } 195 196 /* 197 Get local (ghosted) values of vector 198 */ 199 ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 200 ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 201 ierr = VecGetArrayRead(xlocal,&zctx.v);CHKERRQ(ierr); 202 203 /* 204 Get coordinates of nodes 205 */ 206 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 207 if (!xcoor) { 208 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 209 ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr); 210 } 211 212 /* 213 Determine the min and max coordinates in plot 214 */ 215 ierr = VecStrideMin(xcoor,0,NULL,&zctx.xmin);CHKERRQ(ierr); 216 ierr = VecStrideMax(xcoor,0,NULL,&zctx.xmax);CHKERRQ(ierr); 217 ierr = VecStrideMin(xcoor,1,NULL,&zctx.ymin);CHKERRQ(ierr); 218 ierr = VecStrideMax(xcoor,1,NULL,&zctx.ymax);CHKERRQ(ierr); 219 ierr = PetscOptionsGetBool(NULL,NULL,"-draw_contour_axis",&zctx.showaxis,NULL);CHKERRQ(ierr); 220 if (zctx.showaxis) { 221 coors[0] = zctx.xmin - .05*(zctx.xmax - zctx.xmin); coors[1] = zctx.ymin - .05*(zctx.ymax - zctx.ymin); 222 coors[2] = zctx.xmax + .05*(zctx.xmax - zctx.xmin); coors[3] = zctx.ymax + .05*(zctx.ymax - zctx.ymin); 223 } else { 224 coors[0] = zctx.xmin; coors[1] = zctx.ymin; coors[2] = zctx.xmax; coors[3] = zctx.ymax; 225 } 226 ierr = PetscOptionsGetRealArray(NULL,NULL,"-draw_coordinates",coors,&ncoors,NULL);CHKERRQ(ierr); 227 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); 228 229 /* 230 Get local ghosted version of coordinates 231 */ 232 ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr); 233 if (!xcoorl) { 234 /* create DMDA to get local version of graphics */ 235 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); 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 #undef __FUNCT__ 316 #define __FUNCT__ "VecGetHDF5ChunkSize" 317 static PetscErrorCode VecGetHDF5ChunkSize(DM_DA *da, Vec xin, PetscInt dimension, PetscInt timestep, hsize_t *chunkDims) 318 { 319 PetscMPIInt comm_size; 320 PetscErrorCode ierr; 321 hsize_t chunk_size, target_size, dim; 322 hsize_t vec_size = sizeof(PetscScalar)*da->M*da->N*da->P*da->w; 323 hsize_t avg_local_vec_size,KiB = 1024,MiB = KiB*KiB,GiB = MiB*KiB,min_size = MiB; 324 hsize_t max_chunks = 64*KiB; /* HDF5 internal limitation */ 325 hsize_t max_chunk_size = 4*GiB; /* HDF5 internal limitation */ 326 PetscInt zslices=da->p, yslices=da->n, xslices=da->m; 327 328 PetscFunctionBegin; 329 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)xin), &comm_size);CHKERRQ(ierr); 330 avg_local_vec_size = (hsize_t) ceil(vec_size*1.0/comm_size); /* we will attempt to use this as the chunk size */ 331 332 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))))); 333 /* following line uses sizeof(PetscReal) instead of sizeof(PetscScalar) because the last dimension of chunkDims[] captures the 2* when complex numbers are being used */ 334 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); 335 336 /* 337 if size/rank > max_chunk_size, we need radical measures: even going down to 338 avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter 339 what, composed in the most efficient way possible. 340 N.B. this minimises the number of chunks, which may or may not be the optimal 341 solution. In a BG, for example, the optimal solution is probably to make # chunks = # 342 IO nodes involved, but this author has no access to a BG to figure out how to 343 reliably find the right number. And even then it may or may not be enough. 344 */ 345 if (avg_local_vec_size > max_chunk_size) { 346 /* check if we can just split local z-axis: is that enough? */ 347 zslices = (PetscInt)ceil(vec_size*1.0/(da->p*max_chunk_size))*zslices; 348 if (zslices > da->P) { 349 /* lattice is too large in xy-directions, splitting z only is not enough */ 350 zslices = da->P; 351 yslices= (PetscInt)ceil(vec_size*1.0/(zslices*da->n*max_chunk_size))*yslices; 352 if (yslices > da->N) { 353 /* lattice is too large in x-direction, splitting along z, y is not enough */ 354 yslices = da->N; 355 xslices= (PetscInt)ceil(vec_size*1.0/(zslices*yslices*da->m*max_chunk_size))*xslices; 356 } 357 } 358 dim = 0; 359 if (timestep >= 0) { 360 ++dim; 361 } 362 /* prefer to split z-axis, even down to planar slices */ 363 if (dimension == 3) { 364 chunkDims[dim++] = (hsize_t) da->P/zslices; 365 chunkDims[dim++] = (hsize_t) da->N/yslices; 366 chunkDims[dim++] = (hsize_t) da->M/xslices; 367 } else { 368 /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */ 369 chunkDims[dim++] = (hsize_t) da->N/yslices; 370 chunkDims[dim++] = (hsize_t) da->M/xslices; 371 } 372 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); 373 } else { 374 if (target_size < chunk_size) { 375 /* only change the defaults if target_size < chunk_size */ 376 dim = 0; 377 if (timestep >= 0) { 378 ++dim; 379 } 380 /* prefer to split z-axis, even down to planar slices */ 381 if (dimension == 3) { 382 /* try splitting the z-axis to core-size bits, i.e. divide chunk size by # comm_size in z-direction */ 383 if (target_size >= chunk_size/da->p) { 384 /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */ 385 chunkDims[dim] = (hsize_t) ceil(da->P*1.0/da->p); 386 } else { 387 /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 388 radical and let everyone write all they've got */ 389 chunkDims[dim++] = (hsize_t) ceil(da->P*1.0/da->p); 390 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 391 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 392 } 393 } else { 394 /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */ 395 if (target_size >= chunk_size/da->n) { 396 /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */ 397 chunkDims[dim] = (hsize_t) ceil(da->N*1.0/da->n); 398 } else { 399 /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 400 radical and let everyone write all they've got */ 401 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 402 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 403 } 404 405 } 406 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); 407 } else { 408 /* precomputed chunks are fine, we don't need to do anything */ 409 } 410 } 411 PetscFunctionReturn(0); 412 } 413 #endif 414 415 #if defined(PETSC_HAVE_HDF5) 416 #undef __FUNCT__ 417 #define __FUNCT__ "VecView_MPI_HDF5_DA" 418 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 419 { 420 DM dm; 421 DM_DA *da; 422 hid_t filespace; /* file dataspace identifier */ 423 hid_t chunkspace; /* chunk dataset property identifier */ 424 hid_t plist_id; /* property list identifier */ 425 hid_t dset_id; /* dataset identifier */ 426 hid_t memspace; /* memory dataspace identifier */ 427 hid_t file_id; 428 hid_t group; 429 hid_t memscalartype; /* scalar type for mem (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 430 hid_t filescalartype; /* scalar type for file (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 431 hsize_t dim; 432 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 */ 433 PetscInt timestep, dimension; 434 const PetscScalar *x; 435 const char *vecname; 436 PetscErrorCode ierr; 437 PetscBool dim2; 438 PetscBool spoutput; 439 440 PetscFunctionBegin; 441 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 442 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 443 ierr = PetscViewerHDF5GetBaseDimension2(viewer,&dim2);CHKERRQ(ierr); 444 ierr = PetscViewerHDF5GetSPOutput(viewer,&spoutput);CHKERRQ(ierr); 445 446 ierr = VecGetDM(xin,&dm);CHKERRQ(ierr); 447 if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 448 da = (DM_DA*)dm->data; 449 ierr = DMGetDimension(dm, &dimension);CHKERRQ(ierr); 450 451 /* Create the dataspace for the dataset. 452 * 453 * dims - holds the current dimensions of the dataset 454 * 455 * maxDims - holds the maximum dimensions of the dataset (unlimited 456 * for the number of time steps with the current dimensions for the 457 * other dimensions; so only additional time steps can be added). 458 * 459 * chunkDims - holds the size of a single time step (required to 460 * permit extending dataset). 461 */ 462 dim = 0; 463 if (timestep >= 0) { 464 dims[dim] = timestep+1; 465 maxDims[dim] = H5S_UNLIMITED; 466 chunkDims[dim] = 1; 467 ++dim; 468 } 469 if (dimension == 3) { 470 ierr = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr); 471 maxDims[dim] = dims[dim]; 472 chunkDims[dim] = dims[dim]; 473 ++dim; 474 } 475 if (dimension > 1) { 476 ierr = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr); 477 maxDims[dim] = dims[dim]; 478 chunkDims[dim] = dims[dim]; 479 ++dim; 480 } 481 ierr = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr); 482 maxDims[dim] = dims[dim]; 483 chunkDims[dim] = dims[dim]; 484 ++dim; 485 if (da->w > 1 || dim2) { 486 ierr = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr); 487 maxDims[dim] = dims[dim]; 488 chunkDims[dim] = dims[dim]; 489 ++dim; 490 } 491 #if defined(PETSC_USE_COMPLEX) 492 dims[dim] = 2; 493 maxDims[dim] = dims[dim]; 494 chunkDims[dim] = dims[dim]; 495 ++dim; 496 #endif 497 498 ierr = VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims); CHKERRQ(ierr); 499 500 PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims)); 501 502 #if defined(PETSC_USE_REAL_SINGLE) 503 memscalartype = H5T_NATIVE_FLOAT; 504 filescalartype = H5T_NATIVE_FLOAT; 505 #elif defined(PETSC_USE_REAL___FLOAT128) 506 #error "HDF5 output with 128 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 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 521 PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT)); 522 #else 523 PetscStackCallHDF5Return(dset_id,H5Dcreate,(group, vecname, filescalartype, filespace, H5P_DEFAULT)); 524 #endif 525 } else { 526 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 527 PetscStackCallHDF5(H5Dset_extent,(dset_id, dims)); 528 } 529 PetscStackCallHDF5(H5Sclose,(filespace)); 530 531 /* Each process defines a dataset and writes it to the hyperslab in the file */ 532 dim = 0; 533 if (timestep >= 0) { 534 offset[dim] = timestep; 535 ++dim; 536 } 537 if (dimension == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);} 538 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);} 539 ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr); 540 if (da->w > 1 || dim2) offset[dim++] = 0; 541 #if defined(PETSC_USE_COMPLEX) 542 offset[dim++] = 0; 543 #endif 544 dim = 0; 545 if (timestep >= 0) { 546 count[dim] = 1; 547 ++dim; 548 } 549 if (dimension == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);} 550 if (dimension > 1) {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);} 551 ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr); 552 if (da->w > 1 || dim2) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);} 553 #if defined(PETSC_USE_COMPLEX) 554 count[dim++] = 2; 555 #endif 556 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 557 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 558 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 559 560 /* Create property list for collective dataset write */ 561 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 562 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 563 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 564 #endif 565 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 566 567 ierr = VecGetArrayRead(xin, &x);CHKERRQ(ierr); 568 PetscStackCallHDF5(H5Dwrite,(dset_id, memscalartype, memspace, filespace, plist_id, x)); 569 PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL)); 570 ierr = VecRestoreArrayRead(xin, &x);CHKERRQ(ierr); 571 572 /* Close/release resources */ 573 if (group != file_id) { 574 PetscStackCallHDF5(H5Gclose,(group)); 575 } 576 PetscStackCallHDF5(H5Pclose,(plist_id)); 577 PetscStackCallHDF5(H5Sclose,(filespace)); 578 PetscStackCallHDF5(H5Sclose,(memspace)); 579 PetscStackCallHDF5(H5Dclose,(dset_id)); 580 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 #endif 584 585 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 586 587 #if defined(PETSC_HAVE_MPIIO) 588 #undef __FUNCT__ 589 #define __FUNCT__ "DMDAArrayMPIIO" 590 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 591 { 592 PetscErrorCode ierr; 593 MPI_File mfdes; 594 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 595 MPI_Datatype view; 596 const PetscScalar *array; 597 MPI_Offset off; 598 MPI_Aint ub,ul; 599 PetscInt type,rows,vecrows,tr[2]; 600 DM_DA *dd = (DM_DA*)da->data; 601 PetscBool skipheader; 602 603 PetscFunctionBegin; 604 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 605 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr); 606 if (!write) { 607 /* Read vector header. */ 608 if (!skipheader) { 609 ierr = PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);CHKERRQ(ierr); 610 type = tr[0]; 611 rows = tr[1]; 612 if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file"); 613 if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 614 } 615 } else { 616 tr[0] = VEC_FILE_CLASSID; 617 tr[1] = vecrows; 618 if (!skipheader) { 619 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 620 } 621 } 622 623 ierr = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr); 624 gsizes[0] = dof; 625 ierr = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr); 626 ierr = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr); 627 ierr = PetscMPIIntCast(dd->P,gsizes+3);CHKERRQ(ierr); 628 lsizes[0] = dof; 629 ierr = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr); 630 ierr = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr); 631 ierr = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr); 632 lstarts[0] = 0; 633 ierr = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr); 634 ierr = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr); 635 ierr = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr); 636 ierr = MPI_Type_create_subarray(da->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 637 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 638 639 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 640 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 641 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 642 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 643 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 644 if (write) { 645 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 646 } else { 647 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 648 } 649 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 650 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 651 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 652 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 #endif 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "VecView_MPI_DA" 659 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 660 { 661 DM da; 662 PetscErrorCode ierr; 663 PetscInt dim; 664 Vec natural; 665 PetscBool isdraw,isvtk; 666 #if defined(PETSC_HAVE_HDF5) 667 PetscBool ishdf5; 668 #endif 669 const char *prefix,*name; 670 PetscViewerFormat format; 671 672 PetscFunctionBegin; 673 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 674 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 675 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 676 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 677 #if defined(PETSC_HAVE_HDF5) 678 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 679 #endif 680 if (isdraw) { 681 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 682 if (dim == 1) { 683 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 684 } else if (dim == 2) { 685 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 686 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 687 } else if (isvtk) { /* Duplicate the Vec */ 688 Vec Y; 689 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 690 if (((PetscObject)xin)->name) { 691 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 692 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 693 } 694 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 695 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 696 #if defined(PETSC_HAVE_HDF5) 697 } else if (ishdf5) { 698 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 699 #endif 700 } else { 701 #if defined(PETSC_HAVE_MPIIO) 702 PetscBool isbinary,isMPIIO; 703 704 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 705 if (isbinary) { 706 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 707 if (isMPIIO) { 708 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 709 PetscFunctionReturn(0); 710 } 711 } 712 #endif 713 714 /* call viewer on natural ordering */ 715 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 716 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 717 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 718 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 719 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 720 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 721 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 722 723 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 724 if (format == PETSC_VIEWER_BINARY_MATLAB) { 725 /* temporarily remove viewer format so it won't trigger in the VecView() */ 726 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 727 } 728 729 ierr = VecView(natural,viewer);CHKERRQ(ierr); 730 731 if (format == PETSC_VIEWER_BINARY_MATLAB) { 732 MPI_Comm comm; 733 FILE *info; 734 const char *fieldname; 735 char fieldbuf[256]; 736 PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n; 737 738 /* set the viewer format back into the viewer */ 739 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 740 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 741 ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr); 742 ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr); 743 ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr); 744 ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr); 745 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); } 746 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); } 747 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); } 748 749 for (n=0; n<dof; n++) { 750 ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr); 751 if (!fieldname || !fieldname[0]) { 752 ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr); 753 fieldname = fieldbuf; 754 } 755 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 756 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 757 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);} 758 } 759 ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr); 760 ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr); 761 } 762 763 ierr = VecDestroy(&natural);CHKERRQ(ierr); 764 } 765 PetscFunctionReturn(0); 766 } 767 768 #if defined(PETSC_HAVE_HDF5) 769 #undef __FUNCT__ 770 #define __FUNCT__ "VecLoad_HDF5_DA" 771 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 772 { 773 DM da; 774 PetscErrorCode ierr; 775 hsize_t dim,rdim; 776 hsize_t dims[6]={0},count[6]={0},offset[6]={0}; 777 PetscInt dimension,timestep,dofInd; 778 PetscScalar *x; 779 const char *vecname; 780 hid_t filespace; /* file dataspace identifier */ 781 hid_t plist_id; /* property list identifier */ 782 hid_t dset_id; /* dataset identifier */ 783 hid_t memspace; /* memory dataspace identifier */ 784 hid_t file_id,group; 785 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 786 DM_DA *dd; 787 PetscBool dim2 = PETSC_FALSE; 788 789 PetscFunctionBegin; 790 #if defined(PETSC_USE_REAL_SINGLE) 791 scalartype = H5T_NATIVE_FLOAT; 792 #elif defined(PETSC_USE_REAL___FLOAT128) 793 #error "HDF5 output with 128 bit floats not supported." 794 #else 795 scalartype = H5T_NATIVE_DOUBLE; 796 #endif 797 798 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 799 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 800 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 801 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 802 dd = (DM_DA*)da->data; 803 ierr = DMGetDimension(da, &dimension);CHKERRQ(ierr); 804 805 /* Open dataset */ 806 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 807 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 808 #else 809 PetscStackCallHDF5Return(dset_id,H5Dopen,(group, vecname)); 810 #endif 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 /* Create property list for collective dataset write */ 878 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 879 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 880 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 881 #endif 882 /* To read dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 883 884 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 885 PetscStackCallHDF5(H5Dread,(dset_id, scalartype, memspace, filespace, plist_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(H5Pclose,(plist_id)); 893 PetscStackCallHDF5(H5Sclose,(filespace)); 894 PetscStackCallHDF5(H5Sclose,(memspace)); 895 PetscStackCallHDF5(H5Dclose,(dset_id)); 896 PetscFunctionReturn(0); 897 } 898 #endif 899 900 #undef __FUNCT__ 901 #define __FUNCT__ "VecLoad_Binary_DA" 902 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 903 { 904 DM da; 905 PetscErrorCode ierr; 906 Vec natural; 907 const char *prefix; 908 PetscInt bs; 909 PetscBool flag; 910 DM_DA *dd; 911 #if defined(PETSC_HAVE_MPIIO) 912 PetscBool isMPIIO; 913 #endif 914 915 PetscFunctionBegin; 916 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 917 dd = (DM_DA*)da->data; 918 #if defined(PETSC_HAVE_MPIIO) 919 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 920 if (isMPIIO) { 921 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 922 PetscFunctionReturn(0); 923 } 924 #endif 925 926 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 927 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 928 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 929 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 930 ierr = VecLoad(natural,viewer);CHKERRQ(ierr); 931 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 932 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 933 ierr = VecDestroy(&natural);CHKERRQ(ierr); 934 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 935 ierr = PetscOptionsGetInt(NULL,((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 936 if (flag && bs != dd->w) { 937 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 938 } 939 PetscFunctionReturn(0); 940 } 941 942 #undef __FUNCT__ 943 #define __FUNCT__ "VecLoad_Default_DA" 944 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 945 { 946 PetscErrorCode ierr; 947 DM da; 948 PetscBool isbinary; 949 #if defined(PETSC_HAVE_HDF5) 950 PetscBool ishdf5; 951 #endif 952 953 PetscFunctionBegin; 954 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 955 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 956 957 #if defined(PETSC_HAVE_HDF5) 958 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 959 #endif 960 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 961 962 if (isbinary) { 963 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 964 #if defined(PETSC_HAVE_HDF5) 965 } else if (ishdf5) { 966 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 967 #endif 968 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 969 PetscFunctionReturn(0); 970 } 971