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