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

  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);

  9:    The file storage of the vector is independent of the number of processes used.
 10:  */

 12: #include <petscdm.h>
 13: #include <petscdmda.h>
 14: #include <petscsys.h>
 15: #include <petscviewerhdf5.h>

 17: static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n";

 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;

 33:   dx = (xp - xm) / (Nx - 1);
 34:   dy = (yp - ym) / (Ny - 1);

 36:   /* Initialize the Petsc context */
 37:   PetscFunctionBeginUser;
 38:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 39:   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));
 40:   PetscCall(DMSetFromOptions(da2D));
 41:   PetscCall(DMSetUp(da2D));

 43:   /* Set the coordinates */
 44:   PetscCall(DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));

 46:   /* Declare gauss as a DMDA component */
 47:   PetscCall(DMCreateGlobalVector(da2D, &gauss));
 48:   PetscCall(PetscObjectSetName((PetscObject)gauss, "pressure"));

 50:   /* Initialize vector gauss with a constant value (=1) */
 51:   PetscCall(VecSet(gauss, value));

 53:   /* Get the coordinates of the corners for each process */
 54:   PetscCall(DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0));

 56:   /* Build the gaussian profile (exp(-x^2-y^2)) */
 57:   PetscCall(DMDAVecGetArray(da2D, gauss, &gauss_ptr));
 58:   for (j = iys; j < iys + iym; j++) {
 59:     for (i = ixs; i < ixs + ixm; i++) gauss_ptr[j][i] = PetscExpScalar(-(xm + i * dx) * (xm + i * dx) - (ym + j * dy) * (ym + j * dy));
 60:   }
 61:   PetscCall(DMDAVecRestoreArray(da2D, gauss, &gauss_ptr));

 63:   /* Create the HDF5 viewer */
 64:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_WRITE, &H5viewer));
 65:   PetscCall(PetscViewerSetFromOptions(H5viewer));

 67:   /* Write the H5 file */
 68:   PetscCall(VecView(gauss, H5viewer));

 70:   /* Close the viewer */
 71:   PetscCall(PetscViewerDestroy(&H5viewer));

 73:   PetscCall(VecDuplicate(gauss, &input));
 74:   PetscCall(PetscObjectGetName((PetscObject)gauss, &vecname));
 75:   PetscCall(PetscObjectSetName((PetscObject)input, vecname));

 77:   /* Create the HDF5 viewer for reading */
 78:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_READ, &H5viewer));
 79:   PetscCall(PetscViewerSetFromOptions(H5viewer));
 80:   PetscCall(VecLoad(input, H5viewer));
 81:   PetscCall(PetscViewerDestroy(&H5viewer));

 83:   PetscCall(VecAXPY(input, -1.0, gauss));
 84:   PetscCall(VecNorm(input, NORM_2, &norm));
 85:   PetscCheck(norm <= 1.e-6, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Vec read in does not match vector written out");

 87:   PetscCall(VecDestroy(&input));
 88:   PetscCall(VecDestroy(&gauss));
 89:   PetscCall(DMDestroy(&da2D));
 90:   PetscCall(PetscFinalize());
 91:   return 0;
 92: }

 94: /*TEST

 96:       build:
 97:          requires: hdf5 !defined(PETSC_USE_CXXCOMPLEX)

 99:       test:
100:          nsize: 4

102:       test:
103:          nsize: 4
104:          suffix: 2
105:          args: -viewer_hdf5_base_dimension2
106:          output_file: output/ex10_1.out

108:       test:
109:          nsize: 4
110:          suffix: 3
111:          args: -viewer_hdf5_sp_output
112:          output_file: output/ex10_1.out

114: TEST*/