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 199371c9d4SSatish Balay int main(int argc, char **argv) { 20c4762a1bSJed Brown DM da2D; 21c4762a1bSJed Brown PetscInt i, j, ixs, ixm, iys, iym; 22c4762a1bSJed Brown PetscViewer H5viewer; 23c4762a1bSJed Brown PetscScalar xm = -1.0, xp = 1.0; 24c4762a1bSJed Brown PetscScalar ym = -1.0, yp = 1.0; 25c4762a1bSJed Brown PetscScalar value = 1.0, dx, dy; 26c4762a1bSJed Brown PetscInt Nx = 40, Ny = 40; 27c4762a1bSJed Brown Vec gauss, input; 28c4762a1bSJed Brown PetscScalar **gauss_ptr; 29c4762a1bSJed Brown PetscReal norm; 30c4762a1bSJed Brown const char *vecname; 31c4762a1bSJed Brown 32c4762a1bSJed Brown dx = (xp - xm) / (Nx - 1); 33c4762a1bSJed Brown dy = (yp - ym) / (Ny - 1); 34c4762a1bSJed Brown 35c4762a1bSJed Brown /* Initialize the Petsc context */ 36327415f7SBarry Smith PetscFunctionBeginUser; 379566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 389566063dSJacob Faibussowitsch 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)); 399566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(da2D)); 409566063dSJacob Faibussowitsch PetscCall(DMSetUp(da2D)); 41c4762a1bSJed Brown 42c4762a1bSJed Brown /* Set the coordinates */ 439566063dSJacob Faibussowitsch PetscCall(DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0)); 44c4762a1bSJed Brown 45c4762a1bSJed Brown /* Declare gauss as a DMDA component */ 469566063dSJacob Faibussowitsch PetscCall(DMCreateGlobalVector(da2D, &gauss)); 479566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)gauss, "pressure")); 48c4762a1bSJed Brown 49c4762a1bSJed Brown /* Initialize vector gauss with a constant value (=1) */ 509566063dSJacob Faibussowitsch PetscCall(VecSet(gauss, value)); 51c4762a1bSJed Brown 52c4762a1bSJed Brown /* Get the coordinates of the corners for each process */ 539566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0)); 54c4762a1bSJed Brown 55c4762a1bSJed Brown /* Build the gaussian profile (exp(-x^2-y^2)) */ 569566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da2D, gauss, &gauss_ptr)); 57c4762a1bSJed Brown for (j = iys; j < iys + iym; j++) { 58*ad540459SPierre Jolivet for (i = ixs; i < ixs + ixm; i++) gauss_ptr[j][i] = PetscExpScalar(-(xm + i * dx) * (xm + i * dx) - (ym + j * dy) * (ym + j * dy)); 59c4762a1bSJed Brown } 609566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da2D, gauss, &gauss_ptr)); 61c4762a1bSJed Brown 62c4762a1bSJed Brown /* Create the HDF5 viewer */ 639566063dSJacob Faibussowitsch PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_WRITE, &H5viewer)); 649566063dSJacob Faibussowitsch PetscCall(PetscViewerSetFromOptions(H5viewer)); 65c4762a1bSJed Brown 66c4762a1bSJed Brown /* Write the H5 file */ 679566063dSJacob Faibussowitsch PetscCall(VecView(gauss, H5viewer)); 68c4762a1bSJed Brown 69c4762a1bSJed Brown /* Close the viewer */ 709566063dSJacob Faibussowitsch PetscCall(PetscViewerDestroy(&H5viewer)); 71c4762a1bSJed Brown 729566063dSJacob Faibussowitsch PetscCall(VecDuplicate(gauss, &input)); 739566063dSJacob Faibussowitsch PetscCall(PetscObjectGetName((PetscObject)gauss, &vecname)); 749566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)input, vecname)); 75c4762a1bSJed Brown 76c4762a1bSJed Brown /* Create the HDF5 viewer for reading */ 779566063dSJacob Faibussowitsch PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_READ, &H5viewer)); 789566063dSJacob Faibussowitsch PetscCall(PetscViewerSetFromOptions(H5viewer)); 799566063dSJacob Faibussowitsch PetscCall(VecLoad(input, H5viewer)); 809566063dSJacob Faibussowitsch PetscCall(PetscViewerDestroy(&H5viewer)); 81c4762a1bSJed Brown 829566063dSJacob Faibussowitsch PetscCall(VecAXPY(input, -1.0, gauss)); 839566063dSJacob Faibussowitsch PetscCall(VecNorm(input, NORM_2, &norm)); 8408401ef6SPierre Jolivet PetscCheck(norm <= 1.e-6, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Vec read in does not match vector written out"); 85c4762a1bSJed Brown 869566063dSJacob Faibussowitsch PetscCall(VecDestroy(&input)); 879566063dSJacob Faibussowitsch PetscCall(VecDestroy(&gauss)); 889566063dSJacob Faibussowitsch PetscCall(DMDestroy(&da2D)); 899566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 90b122ec5aSJacob Faibussowitsch return 0; 91c4762a1bSJed Brown } 92c4762a1bSJed Brown 93c4762a1bSJed Brown /*TEST 94c4762a1bSJed Brown 95c4762a1bSJed Brown build: 96dfd57a17SPierre Jolivet requires: hdf5 !defined(PETSC_USE_CXXCOMPLEX) 97c4762a1bSJed Brown 98c4762a1bSJed Brown test: 99c4762a1bSJed Brown nsize: 4 100c4762a1bSJed Brown 101c4762a1bSJed Brown test: 102c4762a1bSJed Brown nsize: 4 103c4762a1bSJed Brown suffix: 2 104c4762a1bSJed Brown args: -viewer_hdf5_base_dimension2 105c4762a1bSJed Brown output_file: output/ex10_1.out 106c4762a1bSJed Brown 107c4762a1bSJed Brown test: 108c4762a1bSJed Brown nsize: 4 109c4762a1bSJed Brown suffix: 3 110c4762a1bSJed Brown args: -viewer_hdf5_sp_output 111c4762a1bSJed Brown output_file: output/ex10_1.out 112c4762a1bSJed Brown 113c4762a1bSJed Brown TEST*/ 114