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,dims[5]; 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 (dd->dim == 3) dims[cnt++] = PetscHDF5IntCast(dd->P); 501 if (dd->dim > 1) dims[cnt++] = PetscHDF5IntCast(dd->N); 502 dims[cnt++] = PetscHDF5IntCast(dd->M); 503 if (dd->w > 1) PetscHDF5IntCast(dims[cnt++] = dd->w); 504 #if defined(PETSC_USE_COMPLEX) 505 dim++; 506 dims[cnt++] = 2; 507 #endif 508 509 /* Create the dataset with default properties and close filespace */ 510 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 511 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 512 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 513 #else 514 dset_id = H5Dopen(file_id, vecname); 515 #endif 516 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 517 filespace = H5Dget_space(dset_id); 518 519 /* Each process defines a dataset and reads it from the hyperslab in the file */ 520 cnt = 0; 521 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 522 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 523 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 524 if (dd->w > 1) offset[cnt++] = 0; 525 #if defined(PETSC_USE_COMPLEX) 526 offset[cnt++] = 0; 527 #endif 528 cnt = 0; 529 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 530 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 531 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 532 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 533 #if defined(PETSC_USE_COMPLEX) 534 count[cnt++] = 2; 535 #endif 536 memspace = H5Screate_simple(dim, count, NULL); 537 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 538 539 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 540 541 /* Create property list for collective dataset write */ 542 plist_id = H5Pcreate(H5P_DATASET_XFER); 543 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 544 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 545 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 546 #endif 547 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 548 549 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 550 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 551 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 552 553 /* Close/release resources */ 554 status = H5Pclose(plist_id);CHKERRQ(status); 555 status = H5Sclose(filespace);CHKERRQ(status); 556 status = H5Sclose(memspace);CHKERRQ(status); 557 status = H5Dclose(dset_id);CHKERRQ(status); 558 PetscFunctionReturn(0); 559 } 560 #endif 561 562 #undef __FUNCT__ 563 #define __FUNCT__ "VecLoad_Binary_DA" 564 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 565 { 566 DM da; 567 PetscErrorCode ierr; 568 Vec natural; 569 const char *prefix; 570 PetscInt bs; 571 PetscBool flag; 572 DM_DA *dd; 573 #if defined(PETSC_HAVE_MPIIO) 574 PetscBool isMPIIO; 575 #endif 576 577 PetscFunctionBegin; 578 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 579 dd = (DM_DA*)da->data; 580 #if defined(PETSC_HAVE_MPIIO) 581 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 582 if (isMPIIO) { 583 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 584 PetscFunctionReturn(0); 585 } 586 #endif 587 588 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 589 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 590 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 591 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 592 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 593 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 594 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 595 ierr = VecDestroy(&natural);CHKERRQ(ierr); 596 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 597 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 598 if (flag && bs != dd->w) { 599 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 600 } 601 PetscFunctionReturn(0); 602 } 603 604 EXTERN_C_BEGIN 605 #undef __FUNCT__ 606 #define __FUNCT__ "VecLoad_Default_DA" 607 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 608 { 609 PetscErrorCode ierr; 610 DM da; 611 PetscBool isbinary; 612 #if defined(PETSC_HAVE_HDF5) 613 PetscBool ishdf5; 614 #endif 615 616 PetscFunctionBegin; 617 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 618 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 619 620 #if defined(PETSC_HAVE_HDF5) 621 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 622 #endif 623 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 624 625 if (isbinary) { 626 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 627 #if defined(PETSC_HAVE_HDF5) 628 } else if (ishdf5) { 629 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 630 #endif 631 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 632 PetscFunctionReturn(0); 633 } 634 EXTERN_C_END 635