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