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