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 DMDABoundaryType 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 != DMDA_BOUNDARY_NONE || by != DMDA_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,DMDA_BOUNDARY_NONE,DMDA_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 != DMDA_BOUNDARY_NONE || by != DMDA_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",coors[0],coors[1],coors[2],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,DMDA_BOUNDARY_NONE,DMDA_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",zctx.min,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 PetscErrorCode VecGetHDF5ChunkSize(DM_DA *da, Vec xin, PetscInt timestep, hsize_t *chunkDims) 324 { 325 int comm_size; 326 hsize_t chunk_size, target_size, dim, 327 vec_size = sizeof(PetscScalar)*da->M*da->N*da->P*da->w, 328 avg_local_vec_size, 329 KiB = 1024, 330 MiB = KiB*KiB, 331 GiB = MiB*KiB, 332 min_size = MiB, 333 max_chunks = 64*KiB, /* HDF5 internal limitation */ 334 max_chunk_size = 4*GiB, /* HDF5 internal limitation */ 335 zslices=da->p, yslices=da->n, xslices=da->m; 336 337 MPI_Comm_size(PetscObjectComm(xin), &comm_size); 338 avg_local_vec_size = (hsize_t) ceil(vec_size*1.0/comm_size); /* we will attempt to use this as the chunk size */ 339 340 target_size = (hsize_t) ceil(PetscMin(vec_size, 341 PetscMin(max_chunk_size, 342 PetscMax(avg_local_vec_size, 343 PetscMax(ceil(vec_size*1.0/max_chunks), 344 min_size) 345 ) 346 ) 347 ) 348 ); 349 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); 350 351 // if size/rank > max_chunk_size, we need radical measures: even going down to 352 // avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter 353 // what, composed in the most efficient way possible. 354 // N.B. this minimises the number of chunks, which may or may not be the optimal 355 // solution. In a BG, for example, the optimal solution is probably to make # chunks = # 356 // IO nodes involved, but this author has no access to a BG to figure out how to 357 // reliably find the right number. And even then it may or may not be enough. 358 if (avg_local_vec_size > max_chunk_size) { 359 // check if we can just split local z-axis: is that enough? 360 zslices=(hsize_t) ceil(vec_size*1.0/(da->p*max_chunk_size))*zslices; 361 if (zslices > da->P) { 362 // lattice is too large in xy-directions, splitting z only is not enough 363 zslices = da->P; 364 yslices=(hsize_t) ceil(vec_size*1.0/(zslices*da->n*max_chunk_size))*yslices; 365 if (yslices > da->N) { 366 // lattice is too large in x-direction, splitting along z, y is not enough 367 yslices = da->N; 368 xslices=(hsize_t) ceil(vec_size*1.0/(zslices*yslices*da->m*max_chunk_size))*xslices; 369 } 370 } 371 dim = 0; 372 if (timestep >= 0) { 373 ++dim; 374 } 375 // prefer to split z-axis, even down to planar slices 376 if (da->dim == 3) { 377 chunkDims[dim++] = (hsize_t) floor(da->P*1.0/zslices); 378 chunkDims[dim++] = (hsize_t) floor(da->N*1.0/yslices); 379 chunkDims[dim++] = (hsize_t) floor(da->M*1.0/xslices); 380 } else { 381 // This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself 382 chunkDims[dim++] = (hsize_t) floor(da->N*1.0/yslices); 383 chunkDims[dim++] = (hsize_t) floor(da->M*1.0/xslices); 384 } 385 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); 386 } else { 387 if (target_size < chunk_size) { 388 // only change the defaults if target_size < chunk_size 389 dim = 0; 390 if (timestep >= 0) { 391 ++dim; 392 } 393 // prefer to split z-axis, even down to planar slices 394 if (da->dim == 3) { 395 // try splitting the z-axis to core-size bits, i.e. divide chunk size by # comm_size in z-direction 396 if (target_size >= chunk_size/da->p) { 397 // just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> 398 chunkDims[dim] = (hsize_t) ceil(da->P*1.0/da->p); 399 } else { 400 // oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 401 // radical and let everyone write all they've got 402 chunkDims[dim++] = (hsize_t) ceil(da->P*1.0/da->p); 403 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 404 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 405 } 406 } else { 407 // This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself 408 if (target_size >= chunk_size/da->n) { 409 // just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> 410 chunkDims[dim] = (hsize_t) ceil(da->N*1.0/da->n); 411 } else { 412 // oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be 413 // radical and let everyone write all they've got 414 chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n); 415 chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m); 416 } 417 418 } 419 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); 420 } else { 421 // precomputed chunks are fine, we don't need to do anything 422 } 423 } 424 PetscFunctionReturn(0); 425 } 426 #endif 427 428 #if defined(PETSC_HAVE_HDF5) 429 #undef __FUNCT__ 430 #define __FUNCT__ "VecView_MPI_HDF5_DA" 431 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 432 { 433 DM dm; 434 DM_DA *da; 435 hid_t filespace; /* file dataspace identifier */ 436 hid_t chunkspace; /* chunk dataset property identifier */ 437 hid_t plist_id; /* property list identifier */ 438 hid_t dset_id; /* dataset identifier */ 439 hid_t memspace; /* memory dataspace identifier */ 440 hid_t file_id; 441 hid_t group; 442 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */ 443 herr_t status; 444 hsize_t dim; 445 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 */ 446 PetscInt timestep; 447 PetscScalar *x; 448 const char *vecname; 449 PetscErrorCode ierr; 450 451 PetscFunctionBegin; 452 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 453 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 454 455 ierr = VecGetDM(xin,&dm);CHKERRQ(ierr); 456 if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 457 da = (DM_DA*)dm->data; 458 459 /* Create the dataspace for the dataset. 460 * 461 * dims - holds the current dimensions of the dataset 462 * 463 * maxDims - holds the maximum dimensions of the dataset (unlimited 464 * for the number of time steps with the current dimensions for the 465 * other dimensions; so only additional time steps can be added). 466 * 467 * chunkDims - holds the size of a single time step (required to 468 * permit extending dataset). 469 */ 470 dim = 0; 471 if (timestep >= 0) { 472 dims[dim] = timestep+1; 473 maxDims[dim] = H5S_UNLIMITED; 474 chunkDims[dim] = 1; 475 ++dim; 476 } 477 if (da->dim == 3) { 478 ierr = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr); 479 maxDims[dim] = dims[dim]; 480 chunkDims[dim] = dims[dim]; 481 ++dim; 482 } 483 if (da->dim > 1) { 484 ierr = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr); 485 maxDims[dim] = dims[dim]; 486 chunkDims[dim] = dims[dim]; 487 ++dim; 488 } 489 ierr = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr); 490 maxDims[dim] = dims[dim]; 491 chunkDims[dim] = dims[dim]; 492 ++dim; 493 if (da->w > 1) { 494 ierr = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr); 495 maxDims[dim] = dims[dim]; 496 chunkDims[dim] = dims[dim]; 497 ++dim; 498 } 499 #if defined(PETSC_USE_COMPLEX) 500 dims[dim] = 2; 501 maxDims[dim] = dims[dim]; 502 chunkDims[dim] = dims[dim]; 503 ++dim; 504 #endif 505 506 VecGetHDF5ChunkSize(da, xin, timestep, chunkDims); 507 508 filespace = H5Screate_simple(dim, dims, maxDims); 509 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 510 511 #if defined(PETSC_USE_REAL_SINGLE) 512 scalartype = H5T_NATIVE_FLOAT; 513 #elif defined(PETSC_USE_REAL___FLOAT128) 514 #error "HDF5 output with 128 bit floats not supported." 515 #else 516 scalartype = H5T_NATIVE_DOUBLE; 517 #endif 518 519 /* Create the dataset with default properties and close filespace */ 520 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 521 if (!H5Lexists(group, vecname, H5P_DEFAULT)) { 522 /* Create chunk */ 523 chunkspace = H5Pcreate(H5P_DATASET_CREATE); 524 if (chunkspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 525 status = H5Pset_chunk(chunkspace, dim, chunkDims);CHKERRQ(status); 526 527 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 528 dset_id = H5Dcreate2(group, vecname, scalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT); 529 #else 530 dset_id = H5Dcreate(group, vecname, scalartype, filespace, H5P_DEFAULT); 531 #endif 532 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 533 } else { 534 dset_id = H5Dopen2(group, vecname, H5P_DEFAULT); 535 status = H5Dset_extent(dset_id, dims);CHKERRQ(status); 536 } 537 status = H5Sclose(filespace);CHKERRQ(status); 538 539 /* Each process defines a dataset and writes it to the hyperslab in the file */ 540 dim = 0; 541 if (timestep >= 0) { 542 offset[dim] = timestep; 543 ++dim; 544 } 545 if (da->dim == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);} 546 if (da->dim > 1) {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);} 547 ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr); 548 if (da->w > 1) offset[dim++] = 0; 549 #if defined(PETSC_USE_COMPLEX) 550 offset[dim++] = 0; 551 #endif 552 dim = 0; 553 if (timestep >= 0) { 554 count[dim] = 1; 555 ++dim; 556 } 557 if (da->dim == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);} 558 if (da->dim > 1) {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);} 559 ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr); 560 if (da->w > 1) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);} 561 #if defined(PETSC_USE_COMPLEX) 562 count[dim++] = 2; 563 #endif 564 memspace = H5Screate_simple(dim, count, NULL); 565 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 566 567 filespace = H5Dget_space(dset_id); 568 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 569 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 570 571 /* Create property list for collective dataset write */ 572 plist_id = H5Pcreate(H5P_DATASET_XFER); 573 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 574 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 575 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 576 #endif 577 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 578 579 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 580 status = H5Dwrite(dset_id, scalartype, memspace, filespace, plist_id, x);CHKERRQ(status); 581 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 582 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 583 584 /* Close/release resources */ 585 if (group != file_id) { 586 status = H5Gclose(group);CHKERRQ(status); 587 } 588 status = H5Pclose(plist_id);CHKERRQ(status); 589 status = H5Sclose(filespace);CHKERRQ(status); 590 status = H5Sclose(memspace);CHKERRQ(status); 591 status = H5Dclose(dset_id);CHKERRQ(status); 592 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 593 PetscFunctionReturn(0); 594 } 595 #endif 596 597 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 598 599 #if defined(PETSC_HAVE_MPIIO) 600 #undef __FUNCT__ 601 #define __FUNCT__ "DMDAArrayMPIIO" 602 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 603 { 604 PetscErrorCode ierr; 605 MPI_File mfdes; 606 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 607 MPI_Datatype view; 608 const PetscScalar *array; 609 MPI_Offset off; 610 MPI_Aint ub,ul; 611 PetscInt type,rows,vecrows,tr[2]; 612 DM_DA *dd = (DM_DA*)da->data; 613 614 PetscFunctionBegin; 615 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 616 if (!write) { 617 /* Read vector header. */ 618 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 619 type = tr[0]; 620 rows = tr[1]; 621 if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file"); 622 if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 623 } else { 624 tr[0] = VEC_FILE_CLASSID; 625 tr[1] = vecrows; 626 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 627 } 628 629 ierr = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr); 630 gsizes[0] = dof; 631 ierr = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr); 632 ierr = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr); 633 ierr = PetscMPIIntCast(dd->P,gsizes+1);CHKERRQ(ierr); 634 lsizes[0] = dof; 635 ierr = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr); 636 ierr = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr); 637 ierr = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr); 638 lstarts[0] = 0; 639 ierr = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr); 640 ierr = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr); 641 ierr = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr); 642 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 643 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 644 645 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 646 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 647 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 648 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 649 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 650 if (write) { 651 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 652 } else { 653 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 654 } 655 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 656 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 657 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 658 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 #endif 662 663 #undef __FUNCT__ 664 #define __FUNCT__ "VecView_MPI_DA" 665 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 666 { 667 DM da; 668 PetscErrorCode ierr; 669 PetscInt dim; 670 Vec natural; 671 PetscBool isdraw,isvtk; 672 #if defined(PETSC_HAVE_HDF5) 673 PetscBool ishdf5; 674 #endif 675 const char *prefix,*name; 676 PetscViewerFormat format; 677 678 PetscFunctionBegin; 679 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 680 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 681 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 682 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 683 #if defined(PETSC_HAVE_HDF5) 684 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 685 #endif 686 if (isdraw) { 687 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 688 if (dim == 1) { 689 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 690 } else if (dim == 2) { 691 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 692 } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 693 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 694 Vec Y; 695 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 696 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 697 if (((PetscObject)xin)->name) { 698 /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */ 699 ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr); 700 } 701 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 702 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr); 703 #if defined(PETSC_HAVE_HDF5) 704 } else if (ishdf5) { 705 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 706 #endif 707 } else { 708 #if defined(PETSC_HAVE_MPIIO) 709 PetscBool isbinary,isMPIIO; 710 711 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 712 if (isbinary) { 713 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 714 if (isMPIIO) { 715 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 716 PetscFunctionReturn(0); 717 } 718 } 719 #endif 720 721 /* call viewer on natural ordering */ 722 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 723 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 724 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 725 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 726 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 727 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 728 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 729 730 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 731 if (format == PETSC_VIEWER_BINARY_MATLAB) { 732 /* temporarily remove viewer format so it won't trigger in the VecView() */ 733 ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); 734 } 735 736 ierr = VecView(natural,viewer);CHKERRQ(ierr); 737 738 if (format == PETSC_VIEWER_BINARY_MATLAB) { 739 MPI_Comm comm; 740 FILE *info; 741 const char *fieldname; 742 char fieldbuf[256]; 743 PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n; 744 745 /* set the viewer format back into the viewer */ 746 ierr = PetscViewerSetFormat(viewer,format);CHKERRQ(ierr); 747 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 748 ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr); 749 ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr); 750 ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr); 751 ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr); 752 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); } 753 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); } 754 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); } 755 756 for (n=0; n<dof; n++) { 757 ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr); 758 if (!fieldname || !fieldname[0]) { 759 ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr); 760 fieldname = fieldbuf; 761 } 762 if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 763 if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); } 764 if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);} 765 } 766 ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr); 767 ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr); 768 } 769 770 ierr = VecDestroy(&natural);CHKERRQ(ierr); 771 } 772 PetscFunctionReturn(0); 773 } 774 775 #if defined(PETSC_HAVE_HDF5) 776 #undef __FUNCT__ 777 #define __FUNCT__ "VecLoad_HDF5_DA" 778 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 779 { 780 DM da; 781 PetscErrorCode ierr; 782 hsize_t dim; 783 hsize_t count[5]; 784 hsize_t offset[5]; 785 PetscInt cnt = 0; 786 PetscScalar *x; 787 const char *vecname; 788 hid_t filespace; /* file dataspace identifier */ 789 hid_t plist_id; /* property list identifier */ 790 hid_t dset_id; /* dataset identifier */ 791 hid_t memspace; /* memory dataspace identifier */ 792 hid_t file_id; 793 herr_t status; 794 DM_DA *dd; 795 796 PetscFunctionBegin; 797 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 798 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 799 dd = (DM_DA*)da->data; 800 801 /* Create the dataspace for the dataset */ 802 ierr = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1),&dim);CHKERRQ(ierr); 803 #if defined(PETSC_USE_COMPLEX) 804 dim++; 805 #endif 806 807 /* Create the dataset with default properties and close filespace */ 808 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 809 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 810 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 811 #else 812 dset_id = H5Dopen(file_id, vecname); 813 #endif 814 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 815 filespace = H5Dget_space(dset_id); 816 817 /* Each process defines a dataset and reads it from the hyperslab in the file */ 818 cnt = 0; 819 if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->zs,offset + cnt++);CHKERRQ(ierr);} 820 if (dd->dim > 1) {ierr = PetscHDF5IntCast(dd->ys,offset + cnt++);CHKERRQ(ierr);} 821 ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + cnt++);CHKERRQ(ierr); 822 if (dd->w > 1) offset[cnt++] = 0; 823 #if defined(PETSC_USE_COMPLEX) 824 offset[cnt++] = 0; 825 #endif 826 cnt = 0; 827 if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + cnt++);CHKERRQ(ierr);} 828 if (dd->dim > 1) {ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + cnt++);CHKERRQ(ierr);} 829 ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + cnt++);CHKERRQ(ierr); 830 if (dd->w > 1) {ierr = PetscHDF5IntCast(dd->w,count + cnt++);CHKERRQ(ierr);} 831 #if defined(PETSC_USE_COMPLEX) 832 count[cnt++] = 2; 833 #endif 834 memspace = H5Screate_simple(dim, count, NULL); 835 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 836 837 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 838 839 /* Create property list for collective dataset write */ 840 plist_id = H5Pcreate(H5P_DATASET_XFER); 841 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 842 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 843 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 844 #endif 845 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 846 847 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 848 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 849 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 850 851 /* Close/release resources */ 852 status = H5Pclose(plist_id);CHKERRQ(status); 853 status = H5Sclose(filespace);CHKERRQ(status); 854 status = H5Sclose(memspace);CHKERRQ(status); 855 status = H5Dclose(dset_id);CHKERRQ(status); 856 PetscFunctionReturn(0); 857 } 858 #endif 859 860 #undef __FUNCT__ 861 #define __FUNCT__ "VecLoad_Binary_DA" 862 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 863 { 864 DM da; 865 PetscErrorCode ierr; 866 Vec natural; 867 const char *prefix; 868 PetscInt bs; 869 PetscBool flag; 870 DM_DA *dd; 871 #if defined(PETSC_HAVE_MPIIO) 872 PetscBool isMPIIO; 873 #endif 874 875 PetscFunctionBegin; 876 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 877 dd = (DM_DA*)da->data; 878 #if defined(PETSC_HAVE_MPIIO) 879 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 880 if (isMPIIO) { 881 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 882 PetscFunctionReturn(0); 883 } 884 #endif 885 886 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 887 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 888 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 889 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 890 ierr = VecLoad(natural,viewer);CHKERRQ(ierr); 891 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 892 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 893 ierr = VecDestroy(&natural);CHKERRQ(ierr); 894 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 895 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 896 if (flag && bs != dd->w) { 897 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 898 } 899 PetscFunctionReturn(0); 900 } 901 902 #undef __FUNCT__ 903 #define __FUNCT__ "VecLoad_Default_DA" 904 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 905 { 906 PetscErrorCode ierr; 907 DM da; 908 PetscBool isbinary; 909 #if defined(PETSC_HAVE_HDF5) 910 PetscBool ishdf5; 911 #endif 912 913 PetscFunctionBegin; 914 ierr = VecGetDM(xin,&da);CHKERRQ(ierr); 915 if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 916 917 #if defined(PETSC_HAVE_HDF5) 918 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 919 #endif 920 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 921 922 if (isbinary) { 923 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 924 #if defined(PETSC_HAVE_HDF5) 925 } else if (ishdf5) { 926 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 927 #endif 928 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 929 PetscFunctionReturn(0); 930 } 931