1 2 /* 3 Plots vectors obtained with DMDACreate2d() 4 */ 5 6 #include <private/daimpl.h> /*I "petscdmda.h" I*/ 7 #include <private/vecimpl.h> 8 9 /* 10 The data that is passed into the graphics callback 11 */ 12 typedef struct { 13 PetscInt m,n,step,k; 14 PetscReal min,max,scale; 15 PetscScalar *xy,*v; 16 PetscBool showgrid; 17 } ZoomCtx; 18 19 /* 20 This does the drawing for one particular field 21 in one particular set of coordinates. It is a callback 22 called from PetscDrawZoom() 23 */ 24 #undef __FUNCT__ 25 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom" 26 PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx) 27 { 28 ZoomCtx *zctx = (ZoomCtx*)ctx; 29 PetscErrorCode ierr; 30 PetscInt m,n,i,j,k,step,id,c1,c2,c3,c4; 31 PetscReal s,min,x1,x2,x3,x4,y_1,y2,y3,y4; 32 PetscScalar *v,*xy; 33 34 PetscFunctionBegin; 35 m = zctx->m; 36 n = zctx->n; 37 step = zctx->step; 38 k = zctx->k; 39 v = zctx->v; 40 xy = zctx->xy; 41 s = zctx->scale; 42 min = zctx->min; 43 44 /* PetscDraw the contour plot patch */ 45 for (j=0; j<n-1; j++) { 46 for (i=0; i<m-1; i++) { 47 #if !defined(PETSC_USE_COMPLEX) 48 id = i+j*m; x1 = xy[2*id];y_1 = xy[2*id+1];c1 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 49 id = i+j*m+1; x2 = xy[2*id];y2 = y_1; c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 50 id = i+j*m+1+m;x3 = x2; y3 = xy[2*id+1];c3 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 51 id = i+j*m+m; x4 = x1; y4 = y3; c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 52 #else 53 id = i+j*m; x1 = PetscRealPart(xy[2*id]);y_1 = PetscRealPart(xy[2*id+1]);c1 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 54 id = i+j*m+1; x2 = PetscRealPart(xy[2*id]);y2 = y_1; c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 55 id = i+j*m+1+m;x3 = x2; y3 = PetscRealPart(xy[2*id+1]);c3 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 56 id = i+j*m+m; x4 = x1; y4 = y3; c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 57 #endif 58 ierr = PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);CHKERRQ(ierr); 59 ierr = PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);CHKERRQ(ierr); 60 if (zctx->showgrid) { 61 ierr = PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);CHKERRQ(ierr); 62 ierr = PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);CHKERRQ(ierr); 63 ierr = PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);CHKERRQ(ierr); 64 ierr = PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);CHKERRQ(ierr); 65 } 66 } 67 } 68 PetscFunctionReturn(0); 69 } 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "VecView_MPI_Draw_DA2d" 73 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer) 74 { 75 DM da,dac,dag; 76 PetscErrorCode ierr; 77 PetscMPIInt rank; 78 PetscInt igstart,N,s,M,istart,isize,jgstart,w; 79 const PetscInt *lx,*ly; 80 PetscReal coors[4],ymin,ymax,xmin,xmax; 81 PetscDraw draw,popup; 82 PetscBool isnull,useports = PETSC_FALSE; 83 MPI_Comm comm; 84 Vec xlocal,xcoor,xcoorl; 85 DMDABoundaryType bx,by; 86 DMDAStencilType st; 87 ZoomCtx zctx; 88 PetscDrawViewPorts *ports = PETSC_NULL; 89 PetscViewerFormat format; 90 91 PetscFunctionBegin; 92 zctx.showgrid = PETSC_FALSE; 93 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 94 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 95 96 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 97 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 98 99 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 100 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 101 102 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr); 103 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,PETSC_NULL);CHKERRQ(ierr); 104 105 /* 106 Obtain a sequential vector that is going to contain the local values plus ONE layer of 107 ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will 108 update the local values pluse ONE layer of ghost values. 109 */ 110 ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr); 111 if (!xlocal) { 112 if (!bx || !by || s != 1 || st != DMDA_STENCIL_BOX) { 113 /* 114 if original da is not of stencil width one, or periodic or not a box stencil then 115 create a special DMDA to handle one level of ghost points for graphics 116 */ 117 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); 118 ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr); 119 } else { 120 /* otherwise we can use the da we already have */ 121 dac = da; 122 } 123 /* create local vector for holding ghosted values used in graphics */ 124 ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 125 if (dac != da) { 126 /* don't keep any public reference of this DMDA, is is only available through xlocal */ 127 DM dacd = dac; 128 ierr = DMDestroy(&dacd);CHKERRQ(ierr); 129 } else { 130 /* remove association between xlocal and da, because below we compose in the opposite 131 direction and if we left this connect we'd get a loop, so the objects could 132 never be destroyed */ 133 ierr = PetscObjectCompose((PetscObject)xlocal,"DM",0);CHKERRQ(ierr); 134 } 135 ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr); 136 ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr); 137 } else { 138 if (!bx && !by && s == 1 && st == DMDA_STENCIL_BOX) { 139 dac = da; 140 } else { 141 ierr = PetscObjectQuery((PetscObject)xlocal,"DM",(PetscObject*)&dac);CHKERRQ(ierr); 142 } 143 } 144 145 /* 146 Get local (ghosted) values of vector 147 */ 148 ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 149 ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 150 ierr = VecGetArray(xlocal,&zctx.v);CHKERRQ(ierr); 151 152 /* get coordinates of nodes */ 153 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 154 if (!xcoor) { 155 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 156 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 157 } 158 159 /* 160 Determine the min and max coordinates in plot 161 */ 162 ierr = VecStrideMin(xcoor,0,PETSC_NULL,&xmin);CHKERRQ(ierr); 163 ierr = VecStrideMax(xcoor,0,PETSC_NULL,&xmax);CHKERRQ(ierr); 164 ierr = VecStrideMin(xcoor,1,PETSC_NULL,&ymin);CHKERRQ(ierr); 165 ierr = VecStrideMax(xcoor,1,PETSC_NULL,&ymax);CHKERRQ(ierr); 166 coors[0] = xmin - .05*(xmax- xmin); coors[2] = xmax + .05*(xmax - xmin); 167 coors[1] = ymin - .05*(ymax- ymin); coors[3] = ymax + .05*(ymax - ymin); 168 ierr = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %G %G %G %G\n",coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 169 170 /* 171 get local ghosted version of coordinates 172 */ 173 ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr); 174 if (!xcoorl) { 175 DM dagd; 176 /* create DMDA to get local version of graphics */ 177 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); 178 ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr); 179 ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr); 180 ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr); 181 dagd = dag; 182 ierr = DMDestroy(&dagd);CHKERRQ(ierr);/* dereference dag */ 183 ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr); 184 } else { 185 ierr = PetscObjectQuery((PetscObject)xcoorl,"DM",(PetscObject*)&dag);CHKERRQ(ierr); 186 } 187 ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 188 ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 189 ierr = VecGetArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 190 191 /* 192 Get information about size of area each processor must do graphics for 193 */ 194 ierr = DMDAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&bx,&by,0,0);CHKERRQ(ierr); 195 ierr = DMDAGetGhostCorners(dac,&igstart,&jgstart,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr); 196 ierr = DMDAGetCorners(dac,&istart,0,0,&isize,0,0);CHKERRQ(ierr); 197 198 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_contour_grid",&zctx.showgrid,PETSC_NULL);CHKERRQ(ierr); 199 200 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 201 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_ports",&useports,PETSC_NULL);CHKERRQ(ierr); 202 if (useports || format == PETSC_VIEWER_DRAW_PORTS){ 203 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 204 ierr = PetscDrawViewPortsCreate(draw,zctx.step,&ports);CHKERRQ(ierr); 205 } 206 /* 207 Loop over each field; drawing each in a different window 208 */ 209 for (zctx.k=0; zctx.k<zctx.step; zctx.k++) { 210 if (useports) { 211 ierr = PetscDrawViewPortsSet(ports,zctx.k);CHKERRQ(ierr); 212 } else { 213 ierr = PetscViewerDrawGetDraw(viewer,zctx.k,&draw);CHKERRQ(ierr); 214 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 215 } 216 217 /* 218 Determine the min and max color in plot 219 */ 220 ierr = VecStrideMin(xin,zctx.k,PETSC_NULL,&zctx.min);CHKERRQ(ierr); 221 ierr = VecStrideMax(xin,zctx.k,PETSC_NULL,&zctx.max);CHKERRQ(ierr); 222 if (zctx.min == zctx.max) { 223 zctx.min -= 1.e-12; 224 zctx.max += 1.e-12; 225 } 226 227 if (!rank) { 228 const char *title; 229 230 ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr); 231 if (title) { 232 ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr); 233 } 234 } 235 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 236 ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr); 237 238 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 239 if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);} 240 241 zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min); 242 243 ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr); 244 } 245 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 246 247 ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 248 ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 252 253 #if defined(PETSC_HAVE_HDF5) 254 #undef __FUNCT__ 255 #define __FUNCT__ "VecView_MPI_HDF5_DA" 256 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 257 { 258 PetscErrorCode ierr; 259 DM dm; 260 DM_DA *da; 261 hsize_t dim,dims[5]; 262 hsize_t count[5]; 263 hsize_t offset[5]; 264 PetscInt cnt = 0; 265 PetscScalar *x; 266 const char *vecname; 267 hid_t filespace; /* file dataspace identifier */ 268 hid_t plist_id; /* property list identifier */ 269 hid_t dset_id; /* dataset identifier */ 270 hid_t memspace; /* memory dataspace identifier */ 271 hid_t file_id; 272 herr_t status; 273 274 PetscFunctionBegin; 275 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 276 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&dm);CHKERRQ(ierr); 277 if (!dm) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 278 da = (DM_DA*)dm->data; 279 280 /* Create the dataspace for the dataset */ 281 dim = PetscHDF5IntCast(da->dim + ((da->w == 1) ? 0 : 1)); 282 if (da->dim == 3) dims[cnt++] = PetscHDF5IntCast(da->P); 283 if (da->dim > 1) dims[cnt++] = PetscHDF5IntCast(da->N); 284 dims[cnt++] = PetscHDF5IntCast(da->M); 285 if (da->w > 1) dims[cnt++] = PetscHDF5IntCast(da->w); 286 #if defined(PETSC_USE_COMPLEX) 287 dim++; 288 dims[cnt++] = 2; 289 #endif 290 filespace = H5Screate_simple(dim, dims, NULL); 291 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 292 293 /* Create the dataset with default properties and close filespace */ 294 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 295 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 296 dset_id = H5Dcreate2(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 297 #else 298 dset_id = H5Dcreate(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT); 299 #endif 300 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 301 status = H5Sclose(filespace);CHKERRQ(status); 302 303 /* Each process defines a dataset and writes it to the hyperslab in the file */ 304 cnt = 0; 305 if (da->dim == 3) offset[cnt++] = PetscHDF5IntCast(da->zs); 306 if (da->dim > 1) offset[cnt++] = PetscHDF5IntCast(da->ys); 307 offset[cnt++] = PetscHDF5IntCast(da->xs/da->w); 308 if (da->w > 1) offset[cnt++] = 0; 309 #if defined(PETSC_USE_COMPLEX) 310 offset[cnt++] = 0; 311 #endif 312 cnt = 0; 313 if (da->dim == 3) count[cnt++] = PetscHDF5IntCast(da->ze - da->zs); 314 if (da->dim > 1) count[cnt++] = PetscHDF5IntCast(da->ye - da->ys); 315 count[cnt++] = PetscHDF5IntCast((da->xe - da->xs)/da->w); 316 if (da->w > 1) count[cnt++] = PetscHDF5IntCast(da->w); 317 #if defined(PETSC_USE_COMPLEX) 318 count[cnt++] = 2; 319 #endif 320 memspace = H5Screate_simple(dim, count, NULL); 321 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 322 323 324 filespace = H5Dget_space(dset_id); 325 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 326 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 327 328 /* Create property list for collective dataset write */ 329 plist_id = H5Pcreate(H5P_DATASET_XFER); 330 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 331 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 332 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 333 #endif 334 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 335 336 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 337 status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 338 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 339 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 340 341 /* Close/release resources */ 342 status = H5Pclose(plist_id);CHKERRQ(status); 343 status = H5Sclose(filespace);CHKERRQ(status); 344 status = H5Sclose(memspace);CHKERRQ(status); 345 status = H5Dclose(dset_id);CHKERRQ(status); 346 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 347 PetscFunctionReturn(0); 348 } 349 #endif 350 351 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 352 353 #if defined(PETSC_HAVE_MPIIO) 354 #undef __FUNCT__ 355 #define __FUNCT__ "DMDAArrayMPIIO" 356 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 357 { 358 PetscErrorCode ierr; 359 MPI_File mfdes; 360 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 361 MPI_Datatype view; 362 const PetscScalar *array; 363 MPI_Offset off; 364 MPI_Aint ub,ul; 365 PetscInt type,rows,vecrows,tr[2]; 366 DM_DA *dd = (DM_DA*)da->data; 367 368 PetscFunctionBegin; 369 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 370 if (!write) { 371 /* Read vector header. */ 372 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 373 type = tr[0]; 374 rows = tr[1]; 375 if (type != VEC_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not vector next in file"); 376 if (rows != vecrows) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 377 } else { 378 tr[0] = VEC_FILE_CLASSID; 379 tr[1] = vecrows; 380 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 381 } 382 383 dof = PetscMPIIntCast(dd->w); 384 gsizes[0] = dof; gsizes[1] = PetscMPIIntCast(dd->M); gsizes[2] = PetscMPIIntCast(dd->N); gsizes[3] = PetscMPIIntCast(dd->P); 385 lsizes[0] = dof;lsizes[1] = PetscMPIIntCast((dd->xe-dd->xs)/dof); lsizes[2] = PetscMPIIntCast(dd->ye-dd->ys); lsizes[3] = PetscMPIIntCast(dd->ze-dd->zs); 386 lstarts[0] = 0; lstarts[1] = PetscMPIIntCast(dd->xs/dof); lstarts[2] = PetscMPIIntCast(dd->ys); lstarts[3] = PetscMPIIntCast(dd->zs); 387 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 388 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 389 390 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 391 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 392 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char *)"native",MPI_INFO_NULL);CHKERRQ(ierr); 393 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 394 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 395 if (write) { 396 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 397 } else { 398 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 399 } 400 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 401 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 402 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 403 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 404 PetscFunctionReturn(0); 405 } 406 #endif 407 408 EXTERN_C_BEGIN 409 #undef __FUNCT__ 410 #define __FUNCT__ "VecView_MPI_DA" 411 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 412 { 413 DM da; 414 PetscErrorCode ierr; 415 PetscInt dim; 416 Vec natural; 417 PetscBool isdraw; 418 #if defined(PETSC_HAVE_HDF5) 419 PetscBool ishdf5; 420 #endif 421 const char *prefix,*name; 422 423 PetscFunctionBegin; 424 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 425 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 426 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 427 #if defined(PETSC_HAVE_HDF5) 428 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 429 #endif 430 if (isdraw) { 431 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 432 if (dim == 1) { 433 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 434 } else if (dim == 2) { 435 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 436 } else { 437 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 438 } 439 #if defined(PETSC_HAVE_HDF5) 440 } else if (ishdf5) { 441 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 442 #endif 443 } else { 444 #if defined(PETSC_HAVE_MPIIO) 445 PetscBool isbinary,isMPIIO; 446 447 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 448 if (isbinary) { 449 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 450 if (isMPIIO) { 451 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 } 455 #endif 456 457 /* call viewer on natural ordering */ 458 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 459 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 460 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 461 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 462 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 463 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 464 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 465 ierr = VecView(natural,viewer);CHKERRQ(ierr); 466 ierr = VecDestroy(&natural);CHKERRQ(ierr); 467 } 468 PetscFunctionReturn(0); 469 } 470 EXTERN_C_END 471 472 #if defined(PETSC_HAVE_HDF5) 473 #undef __FUNCT__ 474 #define __FUNCT__ "VecLoad_HDF5_DA" 475 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 476 { 477 DM da; 478 PetscErrorCode ierr; 479 hsize_t dim; 480 hsize_t count[5]; 481 hsize_t offset[5]; 482 PetscInt cnt = 0; 483 PetscScalar *x; 484 const char *vecname; 485 hid_t filespace; /* file dataspace identifier */ 486 hid_t plist_id; /* property list identifier */ 487 hid_t dset_id; /* dataset identifier */ 488 hid_t memspace; /* memory dataspace identifier */ 489 hid_t file_id; 490 herr_t status; 491 DM_DA *dd; 492 493 PetscFunctionBegin; 494 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 495 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 496 dd = (DM_DA*)da->data; 497 498 /* Create the dataspace for the dataset */ 499 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 500 #if defined(PETSC_USE_COMPLEX) 501 dim++; 502 #endif 503 504 /* Create the dataset with default properties and close filespace */ 505 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 506 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 507 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 508 #else 509 dset_id = H5Dopen(file_id, vecname); 510 #endif 511 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 512 filespace = H5Dget_space(dset_id); 513 514 /* Each process defines a dataset and reads it from the hyperslab in the file */ 515 cnt = 0; 516 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 517 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 518 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 519 if (dd->w > 1) offset[cnt++] = 0; 520 #if defined(PETSC_USE_COMPLEX) 521 offset[cnt++] = 0; 522 #endif 523 cnt = 0; 524 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 525 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 526 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 527 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 528 #if defined(PETSC_USE_COMPLEX) 529 count[cnt++] = 2; 530 #endif 531 memspace = H5Screate_simple(dim, count, NULL); 532 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 533 534 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 535 536 /* Create property list for collective dataset write */ 537 plist_id = H5Pcreate(H5P_DATASET_XFER); 538 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 539 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 540 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 541 #endif 542 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 543 544 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 545 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 546 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 547 548 /* Close/release resources */ 549 status = H5Pclose(plist_id);CHKERRQ(status); 550 status = H5Sclose(filespace);CHKERRQ(status); 551 status = H5Sclose(memspace);CHKERRQ(status); 552 status = H5Dclose(dset_id);CHKERRQ(status); 553 PetscFunctionReturn(0); 554 } 555 #endif 556 557 #undef __FUNCT__ 558 #define __FUNCT__ "VecLoad_Binary_DA" 559 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 560 { 561 DM da; 562 PetscErrorCode ierr; 563 Vec natural; 564 const char *prefix; 565 PetscInt bs; 566 PetscBool flag; 567 DM_DA *dd; 568 #if defined(PETSC_HAVE_MPIIO) 569 PetscBool isMPIIO; 570 #endif 571 572 PetscFunctionBegin; 573 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 574 dd = (DM_DA*)da->data; 575 #if defined(PETSC_HAVE_MPIIO) 576 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 577 if (isMPIIO) { 578 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 579 PetscFunctionReturn(0); 580 } 581 #endif 582 583 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 584 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 585 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 586 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 587 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 588 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 589 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 590 ierr = VecDestroy(&natural);CHKERRQ(ierr); 591 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 592 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 593 if (flag && bs != dd->w) { 594 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 595 } 596 PetscFunctionReturn(0); 597 } 598 599 EXTERN_C_BEGIN 600 #undef __FUNCT__ 601 #define __FUNCT__ "VecLoad_Default_DA" 602 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 603 { 604 PetscErrorCode ierr; 605 DM da; 606 PetscBool isbinary; 607 #if defined(PETSC_HAVE_HDF5) 608 PetscBool ishdf5; 609 #endif 610 611 PetscFunctionBegin; 612 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 613 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 614 615 #if defined(PETSC_HAVE_HDF5) 616 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 617 #endif 618 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 619 620 if (isbinary) { 621 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 622 #if defined(PETSC_HAVE_HDF5) 623 } else if (ishdf5) { 624 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 625 #endif 626 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 627 PetscFunctionReturn(0); 628 } 629 EXTERN_C_END 630