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; 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 ierr = PetscDrawViewPortsDestroy(&ports);CHKERRQ(ierr); 243 244 ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 245 ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249 250 #if defined(PETSC_HAVE_HDF5) 251 #undef __FUNCT__ 252 #define __FUNCT__ "VecView_MPI_HDF5_DA" 253 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 254 { 255 PetscErrorCode ierr; 256 DM dm; 257 DM_DA *da; 258 hsize_t dim,dims[5]; 259 hsize_t count[5]; 260 hsize_t offset[5]; 261 PetscInt cnt = 0; 262 PetscScalar *x; 263 const char *vecname; 264 hid_t filespace; /* file dataspace identifier */ 265 hid_t plist_id; /* property list identifier */ 266 hid_t dset_id; /* dataset identifier */ 267 hid_t memspace; /* memory dataspace identifier */ 268 hid_t file_id; 269 herr_t status; 270 271 PetscFunctionBegin; 272 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 273 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&dm);CHKERRQ(ierr); 274 if (!dm) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 275 da = (DM_DA*)dm->data; 276 277 /* Create the dataspace for the dataset */ 278 dim = PetscHDF5IntCast(da->dim + ((da->w == 1) ? 0 : 1)); 279 if (da->dim == 3) dims[cnt++] = PetscHDF5IntCast(da->P); 280 if (da->dim > 1) dims[cnt++] = PetscHDF5IntCast(da->N); 281 dims[cnt++] = PetscHDF5IntCast(da->M); 282 if (da->w > 1) dims[cnt++] = PetscHDF5IntCast(da->w); 283 #if defined(PETSC_USE_COMPLEX) 284 dim++; 285 dims[cnt++] = 2; 286 #endif 287 filespace = H5Screate_simple(dim, dims, NULL); 288 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 289 290 /* Create the dataset with default properties and close filespace */ 291 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 292 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 293 dset_id = H5Dcreate2(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 294 #else 295 dset_id = H5Dcreate(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT); 296 #endif 297 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 298 status = H5Sclose(filespace);CHKERRQ(status); 299 300 /* Each process defines a dataset and writes it to the hyperslab in the file */ 301 cnt = 0; 302 if (da->dim == 3) offset[cnt++] = PetscHDF5IntCast(da->zs); 303 if (da->dim > 1) offset[cnt++] = PetscHDF5IntCast(da->ys); 304 offset[cnt++] = PetscHDF5IntCast(da->xs/da->w); 305 if (da->w > 1) offset[cnt++] = 0; 306 #if defined(PETSC_USE_COMPLEX) 307 offset[cnt++] = 0; 308 #endif 309 cnt = 0; 310 if (da->dim == 3) count[cnt++] = PetscHDF5IntCast(da->ze - da->zs); 311 if (da->dim > 1) count[cnt++] = PetscHDF5IntCast(da->ye - da->ys); 312 count[cnt++] = PetscHDF5IntCast((da->xe - da->xs)/da->w); 313 if (da->w > 1) count[cnt++] = PetscHDF5IntCast(da->w); 314 #if defined(PETSC_USE_COMPLEX) 315 count[cnt++] = 2; 316 #endif 317 memspace = H5Screate_simple(dim, count, NULL); 318 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 319 320 321 filespace = H5Dget_space(dset_id); 322 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 323 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 324 325 /* Create property list for collective dataset write */ 326 plist_id = H5Pcreate(H5P_DATASET_XFER); 327 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 328 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 329 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 330 #endif 331 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 332 333 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 334 status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 335 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 336 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 337 338 /* Close/release resources */ 339 status = H5Pclose(plist_id);CHKERRQ(status); 340 status = H5Sclose(filespace);CHKERRQ(status); 341 status = H5Sclose(memspace);CHKERRQ(status); 342 status = H5Dclose(dset_id);CHKERRQ(status); 343 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 344 PetscFunctionReturn(0); 345 } 346 #endif 347 348 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 349 350 #if defined(PETSC_HAVE_MPIIO) 351 #undef __FUNCT__ 352 #define __FUNCT__ "DMDAArrayMPIIO" 353 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 354 { 355 PetscErrorCode ierr; 356 MPI_File mfdes; 357 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 358 MPI_Datatype view; 359 const PetscScalar *array; 360 MPI_Offset off; 361 MPI_Aint ub,ul; 362 PetscInt type,rows,vecrows,tr[2]; 363 DM_DA *dd = (DM_DA*)da->data; 364 365 PetscFunctionBegin; 366 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 367 if (!write) { 368 /* Read vector header. */ 369 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 370 type = tr[0]; 371 rows = tr[1]; 372 if (type != VEC_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not vector next in file"); 373 if (rows != vecrows) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 374 } else { 375 tr[0] = VEC_FILE_CLASSID; 376 tr[1] = vecrows; 377 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 378 } 379 380 dof = PetscMPIIntCast(dd->w); 381 gsizes[0] = dof; gsizes[1] = PetscMPIIntCast(dd->M); gsizes[2] = PetscMPIIntCast(dd->N); gsizes[3] = PetscMPIIntCast(dd->P); 382 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); 383 lstarts[0] = 0; lstarts[1] = PetscMPIIntCast(dd->xs/dof); lstarts[2] = PetscMPIIntCast(dd->ys); lstarts[3] = PetscMPIIntCast(dd->zs); 384 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 385 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 386 387 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 388 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 389 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char *)"native",MPI_INFO_NULL);CHKERRQ(ierr); 390 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 391 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 392 if (write) { 393 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 394 } else { 395 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 396 } 397 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 398 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 399 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 400 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 #endif 404 405 EXTERN_C_BEGIN 406 #undef __FUNCT__ 407 #define __FUNCT__ "VecView_MPI_DA" 408 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 409 { 410 DM da; 411 PetscErrorCode ierr; 412 PetscInt dim; 413 Vec natural; 414 PetscBool isdraw; 415 #if defined(PETSC_HAVE_HDF5) 416 PetscBool ishdf5; 417 #endif 418 const char *prefix; 419 420 PetscFunctionBegin; 421 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 422 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 423 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 424 #if defined(PETSC_HAVE_HDF5) 425 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 426 #endif 427 if (isdraw) { 428 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 429 if (dim == 1) { 430 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 431 } else if (dim == 2) { 432 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 433 } else { 434 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 435 } 436 #if defined(PETSC_HAVE_HDF5) 437 } else if (ishdf5) { 438 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 439 #endif 440 } else { 441 #if defined(PETSC_HAVE_MPIIO) 442 PetscBool isbinary,isMPIIO; 443 444 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 445 if (isbinary) { 446 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 447 if (isMPIIO) { 448 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 449 PetscFunctionReturn(0); 450 } 451 } 452 #endif 453 454 /* call viewer on natural ordering */ 455 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 456 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 457 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 458 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 459 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 460 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 461 ierr = VecView(natural,viewer);CHKERRQ(ierr); 462 ierr = VecDestroy(&natural);CHKERRQ(ierr); 463 } 464 PetscFunctionReturn(0); 465 } 466 EXTERN_C_END 467 468 #if defined(PETSC_HAVE_HDF5) 469 #undef __FUNCT__ 470 #define __FUNCT__ "VecLoad_HDF5_DA" 471 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 472 { 473 DM da; 474 PetscErrorCode ierr; 475 hsize_t dim,dims[5]; 476 hsize_t count[5]; 477 hsize_t offset[5]; 478 PetscInt cnt = 0; 479 PetscScalar *x; 480 const char *vecname; 481 hid_t filespace; /* file dataspace identifier */ 482 hid_t plist_id; /* property list identifier */ 483 hid_t dset_id; /* dataset identifier */ 484 hid_t memspace; /* memory dataspace identifier */ 485 hid_t file_id; 486 herr_t status; 487 DM_DA *dd; 488 489 PetscFunctionBegin; 490 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 491 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 492 dd = (DM_DA*)da->data; 493 494 /* Create the dataspace for the dataset */ 495 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 496 if (dd->dim == 3) dims[cnt++] = PetscHDF5IntCast(dd->P); 497 if (dd->dim > 1) dims[cnt++] = PetscHDF5IntCast(dd->N); 498 dims[cnt++] = PetscHDF5IntCast(dd->M); 499 if (dd->w > 1) PetscHDF5IntCast(dims[cnt++] = dd->w); 500 #if defined(PETSC_USE_COMPLEX) 501 dim++; 502 dims[cnt++] = 2; 503 #endif 504 505 /* Create the dataset with default properties and close filespace */ 506 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 507 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 508 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 509 #else 510 dset_id = H5Dopen(file_id, vecname); 511 #endif 512 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 513 filespace = H5Dget_space(dset_id); 514 515 /* Each process defines a dataset and reads it from the hyperslab in the file */ 516 cnt = 0; 517 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 518 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 519 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 520 if (dd->w > 1) offset[cnt++] = 0; 521 #if defined(PETSC_USE_COMPLEX) 522 offset[cnt++] = 0; 523 #endif 524 cnt = 0; 525 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 526 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 527 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 528 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 529 #if defined(PETSC_USE_COMPLEX) 530 count[cnt++] = 2; 531 #endif 532 memspace = H5Screate_simple(dim, count, NULL); 533 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 534 535 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 536 537 /* Create property list for collective dataset write */ 538 plist_id = H5Pcreate(H5P_DATASET_XFER); 539 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 540 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 541 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 542 #endif 543 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 544 545 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 546 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 547 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 548 549 /* Close/release resources */ 550 status = H5Pclose(plist_id);CHKERRQ(status); 551 status = H5Sclose(filespace);CHKERRQ(status); 552 status = H5Sclose(memspace);CHKERRQ(status); 553 status = H5Dclose(dset_id);CHKERRQ(status); 554 PetscFunctionReturn(0); 555 } 556 #endif 557 558 #undef __FUNCT__ 559 #define __FUNCT__ "VecLoad_Binary_DA" 560 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 561 { 562 DM da; 563 PetscErrorCode ierr; 564 Vec natural; 565 const char *prefix; 566 PetscInt bs; 567 PetscBool flag; 568 DM_DA *dd; 569 #if defined(PETSC_HAVE_MPIIO) 570 PetscBool isMPIIO; 571 #endif 572 573 PetscFunctionBegin; 574 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 575 dd = (DM_DA*)da->data; 576 #if defined(PETSC_HAVE_MPIIO) 577 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 578 if (isMPIIO) { 579 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 #endif 583 584 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 585 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 586 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 587 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 588 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 589 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 590 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 591 ierr = VecDestroy(&natural);CHKERRQ(ierr); 592 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 593 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 594 if (flag && bs != dd->w) { 595 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 596 } 597 PetscFunctionReturn(0); 598 } 599 600 EXTERN_C_BEGIN 601 #undef __FUNCT__ 602 #define __FUNCT__ "VecLoad_Default_DA" 603 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 604 { 605 PetscErrorCode ierr; 606 DM da; 607 PetscBool isbinary; 608 #if defined(PETSC_HAVE_HDF5) 609 PetscBool ishdf5; 610 #endif 611 612 PetscFunctionBegin; 613 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 614 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 615 616 #if defined(PETSC_HAVE_HDF5) 617 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 618 #endif 619 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 620 621 if (isbinary) { 622 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 623 #if defined(PETSC_HAVE_HDF5) 624 } else if (ishdf5) { 625 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 626 #endif 627 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 628 PetscFunctionReturn(0); 629 } 630 EXTERN_C_END 631