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