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