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