1 /* 2 Demonstrates using the HDF5 viewer with a DMDA Vec 3 - create a global vector containing a gauss profile (exp(-x^2-y^2)) 4 - write the global vector in a hdf5 file 5 6 The resulting file gauss.h5 can be viewed with Visit (an open source visualization package) 7 Or with some versions of MATLAB with data=hdfread('gauss.h5','pressure'); mesh(data); 8 9 The file storage of the vector is independent of the number of processes used. 10 */ 11 12 #include <petscdm.h> 13 #include <petscdmda.h> 14 #include <petscsys.h> 15 #include <petscviewerhdf5.h> 16 17 static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n"; 18 19 int main(int argc,char **argv) 20 { 21 PetscErrorCode ierr; 22 DM da2D; 23 PetscInt i,j,ixs, ixm, iys, iym; 24 PetscViewer H5viewer; 25 PetscScalar xm = -1.0, xp=1.0; 26 PetscScalar ym = -1.0, yp=1.0; 27 PetscScalar value = 1.0,dx,dy; 28 PetscInt Nx = 40, Ny=40; 29 Vec gauss,input; 30 PetscScalar **gauss_ptr; 31 PetscReal norm; 32 const char *vecname; 33 34 dx=(xp-xm)/(Nx-1); 35 dy=(yp-ym)/(Ny-1); 36 37 /* Initialize the Petsc context */ 38 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 39 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,Nx,Ny,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da2D);CHKERRQ(ierr); 40 ierr = DMSetFromOptions(da2D);CHKERRQ(ierr); 41 ierr = DMSetUp(da2D);CHKERRQ(ierr); 42 43 /* Set the coordinates */ 44 ierr = DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);CHKERRQ(ierr); 45 46 /* Declare gauss as a DMDA component */ 47 ierr = DMCreateGlobalVector(da2D,&gauss);CHKERRQ(ierr); 48 ierr = PetscObjectSetName((PetscObject) gauss, "pressure");CHKERRQ(ierr); 49 50 /* Initialize vector gauss with a constant value (=1) */ 51 ierr = VecSet(gauss,value);CHKERRQ(ierr); 52 53 /* Get the coordinates of the corners for each process */ 54 ierr = DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0);CHKERRQ(ierr); 55 56 /* Build the gaussian profile (exp(-x^2-y^2)) */ 57 ierr = DMDAVecGetArray(da2D,gauss,&gauss_ptr);CHKERRQ(ierr); 58 for (j=iys; j<iys+iym; j++) { 59 for (i=ixs; i<ixs+ixm; i++) { 60 gauss_ptr[j][i]=PetscExpScalar(-(xm+i*dx)*(xm+i*dx)-(ym+j*dy)*(ym+j*dy)); 61 } 62 } 63 ierr = DMDAVecRestoreArray(da2D,gauss,&gauss_ptr);CHKERRQ(ierr); 64 65 /* Create the HDF5 viewer */ 66 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_WRITE,&H5viewer);CHKERRQ(ierr); 67 ierr = PetscViewerSetFromOptions(H5viewer);CHKERRQ(ierr); 68 69 /* Write the H5 file */ 70 ierr = VecView(gauss,H5viewer);CHKERRQ(ierr); 71 72 /* Close the viewer */ 73 ierr = PetscViewerDestroy(&H5viewer);CHKERRQ(ierr); 74 75 ierr = VecDuplicate(gauss,&input);CHKERRQ(ierr); 76 ierr = PetscObjectGetName((PetscObject)gauss,&vecname);CHKERRQ(ierr); 77 ierr = PetscObjectSetName((PetscObject)input,vecname);CHKERRQ(ierr); 78 79 /* Create the HDF5 viewer for reading */ 80 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_READ,&H5viewer);CHKERRQ(ierr); 81 ierr = PetscViewerSetFromOptions(H5viewer);CHKERRQ(ierr); 82 ierr = VecLoad(input,H5viewer);CHKERRQ(ierr); 83 ierr = PetscViewerDestroy(&H5viewer);CHKERRQ(ierr); 84 85 ierr = VecAXPY(input,-1.0,gauss);CHKERRQ(ierr); 86 ierr = VecNorm(input,NORM_2,&norm);CHKERRQ(ierr); 87 if (norm > 1.e-6) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Vec read in does not match vector written out"); 88 89 ierr = VecDestroy(&input);CHKERRQ(ierr); 90 ierr = VecDestroy(&gauss);CHKERRQ(ierr); 91 ierr = DMDestroy(&da2D);CHKERRQ(ierr); 92 ierr = PetscFinalize(); 93 return ierr; 94 } 95 96 97 /*TEST 98 99 build: 100 requires: hdf5 !define(PETSC_USE_CXXCOMPLEX) 101 102 test: 103 nsize: 4 104 105 test: 106 nsize: 4 107 suffix: 2 108 args: -viewer_hdf5_base_dimension2 109 output_file: output/ex10_1.out 110 111 test: 112 nsize: 4 113 suffix: 3 114 args: -viewer_hdf5_sp_output 115 output_file: output/ex10_1.out 116 117 TEST*/ 118