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