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