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 DM da2D; 22 PetscInt i,j,ixs, ixm, iys, iym; 23 PetscViewer H5viewer; 24 PetscScalar xm = -1.0, xp=1.0; 25 PetscScalar ym = -1.0, yp=1.0; 26 PetscScalar value = 1.0,dx,dy; 27 PetscInt Nx = 40, Ny=40; 28 Vec gauss,input; 29 PetscScalar **gauss_ptr; 30 PetscReal norm; 31 const char *vecname; 32 33 dx=(xp-xm)/(Nx-1); 34 dy=(yp-ym)/(Ny-1); 35 36 /* Initialize the Petsc context */ 37 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 38 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,Nx,Ny,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da2D)); 39 PetscCall(DMSetFromOptions(da2D)); 40 PetscCall(DMSetUp(da2D)); 41 42 /* Set the coordinates */ 43 PetscCall(DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0)); 44 45 /* Declare gauss as a DMDA component */ 46 PetscCall(DMCreateGlobalVector(da2D,&gauss)); 47 PetscCall(PetscObjectSetName((PetscObject) gauss, "pressure")); 48 49 /* Initialize vector gauss with a constant value (=1) */ 50 PetscCall(VecSet(gauss,value)); 51 52 /* Get the coordinates of the corners for each process */ 53 PetscCall(DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0)); 54 55 /* Build the gaussian profile (exp(-x^2-y^2)) */ 56 PetscCall(DMDAVecGetArray(da2D,gauss,&gauss_ptr)); 57 for (j=iys; j<iys+iym; j++) { 58 for (i=ixs; i<ixs+ixm; i++) { 59 gauss_ptr[j][i]=PetscExpScalar(-(xm+i*dx)*(xm+i*dx)-(ym+j*dy)*(ym+j*dy)); 60 } 61 } 62 PetscCall(DMDAVecRestoreArray(da2D,gauss,&gauss_ptr)); 63 64 /* Create the HDF5 viewer */ 65 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_WRITE,&H5viewer)); 66 PetscCall(PetscViewerSetFromOptions(H5viewer)); 67 68 /* Write the H5 file */ 69 PetscCall(VecView(gauss,H5viewer)); 70 71 /* Close the viewer */ 72 PetscCall(PetscViewerDestroy(&H5viewer)); 73 74 PetscCall(VecDuplicate(gauss,&input)); 75 PetscCall(PetscObjectGetName((PetscObject)gauss,&vecname)); 76 PetscCall(PetscObjectSetName((PetscObject)input,vecname)); 77 78 /* Create the HDF5 viewer for reading */ 79 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_READ,&H5viewer)); 80 PetscCall(PetscViewerSetFromOptions(H5viewer)); 81 PetscCall(VecLoad(input,H5viewer)); 82 PetscCall(PetscViewerDestroy(&H5viewer)); 83 84 PetscCall(VecAXPY(input,-1.0,gauss)); 85 PetscCall(VecNorm(input,NORM_2,&norm)); 86 PetscCheck(norm <= 1.e-6,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Vec read in does not match vector written out"); 87 88 PetscCall(VecDestroy(&input)); 89 PetscCall(VecDestroy(&gauss)); 90 PetscCall(DMDestroy(&da2D)); 91 PetscCall(PetscFinalize()); 92 return 0; 93 } 94 95 /*TEST 96 97 build: 98 requires: hdf5 !defined(PETSC_USE_CXXCOMPLEX) 99 100 test: 101 nsize: 4 102 103 test: 104 nsize: 4 105 suffix: 2 106 args: -viewer_hdf5_base_dimension2 107 output_file: output/ex10_1.out 108 109 test: 110 nsize: 4 111 suffix: 3 112 args: -viewer_hdf5_sp_output 113 output_file: output/ex10_1.out 114 115 TEST*/ 116