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