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