Actual source code: ex10.c

      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, NULL, 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*/