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