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 and hold a reference to the DM */ 688 Vec Y; 689 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 690 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 691 if (((PetscObject)xin)->name) { 692 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 693 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 694 } 695 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 696 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 697 #if defined(PETSC_HAVE_HDF5) 698 } else if (ishdf5) { 699 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 700 #endif 701 } else { 702 #if defined(PETSC_HAVE_MPIIO) 703 PetscBool isbinary,isMPIIO; 704 705 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 706 if (isbinary) { 707 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 708 if (isMPIIO) { 709 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 710 PetscFunctionReturn(0); 711 } 712 } 713 #endif 714 715 /* call viewer on natural ordering */ 716 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 717 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 718 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 719 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 720 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 721 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 722 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 723 724 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 725 if (format == PETSC_VIEWER_BINARY_MATLAB) { 726 /* temporarily remove viewer format so it won't trigger in the VecView() */ 727 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 728 } 729 730 ierr = VecView(natural,viewer);CHKERRQ(ierr); 731 732 if (format == PETSC_VIEWER_BINARY_MATLAB) { 733 MPI_Comm comm; 734 FILE *info; 735 const char *fieldname; 736 char fieldbuf[256]; 737 PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n; 738 739 /* set the viewer format back into the viewer */ 740 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 741 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 742 ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr); 743 ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr); 744 ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr); 745 ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr); 746 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); } 747 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); } 748 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); } 749 750 for (n=0; n<dof; n++) { 751 ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr); 752 if (!fieldname || !fieldname[0]) { 753 ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr); 754 fieldname = fieldbuf; 755 } 756 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 757 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 758 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);} 759 } 760 ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr); 761 ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr); 762 } 763 764 ierr = VecDestroy(&natural);CHKERRQ(ierr); 765 } 766 PetscFunctionReturn(0); 767 } 768 769 #if defined(PETSC_HAVE_HDF5) 770 #undef __FUNCT__ 771 #define __FUNCT__ "VecLoad_HDF5_DA" 772 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 773 { 774 DM da; 775 PetscErrorCode ierr; 776 hsize_t dim,rdim; 777 hsize_t dims[6]={0},count[6]={0},offset[6]={0}; 778 PetscInt dimension,timestep,dofInd; 779 PetscScalar *x; 780 const char *vecname; 781 hid_t filespace; /* file dataspace identifier */ 782 hid_t plist_id; /* property list identifier */ 783 hid_t dset_id; /* dataset identifier */ 784 hid_t memspace; /* memory dataspace identifier */ 785 hid_t file_id,group; 786 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 787 DM_DA *dd; 788 PetscBool dim2 = PETSC_FALSE; 789 790 PetscFunctionBegin; 791 #if defined(PETSC_USE_REAL_SINGLE) 792 scalartype = H5T_NATIVE_FLOAT; 793 #elif defined(PETSC_USE_REAL___FLOAT128) 794 #error "HDF5 output with 128 bit floats not supported." 795 #else 796 scalartype = H5T_NATIVE_DOUBLE; 797 #endif 798 799 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 800 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 801 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 802 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 803 dd = (DM_DA*)da->data; 804 ierr = DMGetDimension(da, &dimension);CHKERRQ(ierr); 805 806 /* Open dataset */ 807 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 808 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT)); 809 #else 810 PetscStackCallHDF5Return(dset_id,H5Dopen,(group, vecname)); 811 #endif 812 813 /* Retrieve the dataspace for the dataset */ 814 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 815 PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL)); 816 817 /* Expected dimension for holding the dof's */ 818 #if defined(PETSC_USE_COMPLEX) 819 dofInd = rdim-2; 820 #else 821 dofInd = rdim-1; 822 #endif 823 824 /* The expected number of dimensions, assuming basedimension2 = false */ 825 dim = dimension; 826 if (dd->w > 1) ++dim; 827 if (timestep >= 0) ++dim; 828 #if defined(PETSC_USE_COMPLEX) 829 ++dim; 830 #endif 831 832 /* In this case the input dataset have one extra, unexpected dimension. */ 833 if (rdim == dim+1) { 834 /* In this case the block size unity */ 835 if (dd->w == 1 && dims[dofInd] == 1) dim2 = PETSC_TRUE; 836 837 /* Special error message for the case where dof does not match the input file */ 838 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); 839 840 /* Other cases where rdim != dim cannot be handled currently */ 841 } 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); 842 843 /* Set up the hyperslab size */ 844 dim = 0; 845 if (timestep >= 0) { 846 offset[dim] = timestep; 847 count[dim] = 1; 848 ++dim; 849 } 850 if (dimension == 3) { 851 ierr = PetscHDF5IntCast(dd->zs,offset + dim);CHKERRQ(ierr); 852 ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + dim);CHKERRQ(ierr); 853 ++dim; 854 } 855 if (dimension > 1) { 856 ierr = PetscHDF5IntCast(dd->ys,offset + dim);CHKERRQ(ierr); 857 ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + dim);CHKERRQ(ierr); 858 ++dim; 859 } 860 ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + dim);CHKERRQ(ierr); 861 ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + dim);CHKERRQ(ierr); 862 ++dim; 863 if (dd->w > 1 || dim2) { 864 offset[dim] = 0; 865 ierr = PetscHDF5IntCast(dd->w,count + dim);CHKERRQ(ierr); 866 ++dim; 867 } 868 #if defined(PETSC_USE_COMPLEX) 869 offset[dim] = 0; 870 count[dim] = 2; 871 ++dim; 872 #endif 873 874 /* Create the memory and filespace */ 875 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 876 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 877 878 /* Create property list for collective dataset write */ 879 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 880 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 881 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 882 #endif 883 /* To read dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 884 885 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 886 PetscStackCallHDF5(H5Dread,(dset_id, scalartype, memspace, filespace, plist_id, x)); 887 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 888 889 /* Close/release resources */ 890 if (group != file_id) { 891 PetscStackCallHDF5(H5Gclose,(group)); 892 } 893 PetscStackCallHDF5(H5Pclose,(plist_id)); 894 PetscStackCallHDF5(H5Sclose,(filespace)); 895 PetscStackCallHDF5(H5Sclose,(memspace)); 896 PetscStackCallHDF5(H5Dclose,(dset_id)); 897 PetscFunctionReturn(0); 898 } 899 #endif 900 901 #undef __FUNCT__ 902 #define __FUNCT__ "VecLoad_Binary_DA" 903 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 904 { 905 DM da; 906 PetscErrorCode ierr; 907 Vec natural; 908 const char *prefix; 909 PetscInt bs; 910 PetscBool flag; 911 DM_DA *dd; 912 #if defined(PETSC_HAVE_MPIIO) 913 PetscBool isMPIIO; 914 #endif 915 916 PetscFunctionBegin; 917 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 918 dd = (DM_DA*)da->data; 919 #if defined(PETSC_HAVE_MPIIO) 920 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 921 if (isMPIIO) { 922 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 923 PetscFunctionReturn(0); 924 } 925 #endif 926 927 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 928 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 929 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 930 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 931 ierr = VecLoad(natural,viewer);CHKERRQ(ierr); 932 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 933 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 934 ierr = VecDestroy(&natural);CHKERRQ(ierr); 935 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 936 ierr = PetscOptionsGetInt(NULL,((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 937 if (flag && bs != dd->w) { 938 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 939 } 940 PetscFunctionReturn(0); 941 } 942 943 #undef __FUNCT__ 944 #define __FUNCT__ "VecLoad_Default_DA" 945 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 946 { 947 PetscErrorCode ierr; 948 DM da; 949 PetscBool isbinary; 950 #if defined(PETSC_HAVE_HDF5) 951 PetscBool ishdf5; 952 #endif 953 954 PetscFunctionBegin; 955 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 956 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 957 958 #if defined(PETSC_HAVE_HDF5) 959 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 960 #endif 961 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 962 963 if (isbinary) { 964 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 965 #if defined(PETSC_HAVE_HDF5) 966 } else if (ishdf5) { 967 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 968 #endif 969 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 970 PetscFunctionReturn(0); 971 } 972