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