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 DM da2D; 21 PetscInt i, j, ixs, ixm, iys, iym; 22 PetscViewer H5viewer; 23 PetscScalar xm = -1.0, xp = 1.0; 24 PetscScalar ym = -1.0, yp = 1.0; 25 PetscScalar value = 1.0, dx, dy; 26 PetscInt Nx = 40, Ny = 40; 27 Vec gauss, input; 28 PetscScalar **gauss_ptr; 29 PetscReal norm; 30 const char *vecname; 31 32 dx = (xp - xm) / (Nx - 1); 33 dy = (yp - ym) / (Ny - 1); 34 35 /* Initialize the Petsc context */ 36 PetscFunctionBeginUser; 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++) { gauss_ptr[j][i] = PetscExpScalar(-(xm + i * dx) * (xm + i * dx) - (ym + j * dy) * (ym + j * dy)); } 59 } 60 PetscCall(DMDAVecRestoreArray(da2D, gauss, &gauss_ptr)); 61 62 /* Create the HDF5 viewer */ 63 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_WRITE, &H5viewer)); 64 PetscCall(PetscViewerSetFromOptions(H5viewer)); 65 66 /* Write the H5 file */ 67 PetscCall(VecView(gauss, H5viewer)); 68 69 /* Close the viewer */ 70 PetscCall(PetscViewerDestroy(&H5viewer)); 71 72 PetscCall(VecDuplicate(gauss, &input)); 73 PetscCall(PetscObjectGetName((PetscObject)gauss, &vecname)); 74 PetscCall(PetscObjectSetName((PetscObject)input, vecname)); 75 76 /* Create the HDF5 viewer for reading */ 77 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_READ, &H5viewer)); 78 PetscCall(PetscViewerSetFromOptions(H5viewer)); 79 PetscCall(VecLoad(input, H5viewer)); 80 PetscCall(PetscViewerDestroy(&H5viewer)); 81 82 PetscCall(VecAXPY(input, -1.0, gauss)); 83 PetscCall(VecNorm(input, NORM_2, &norm)); 84 PetscCheck(norm <= 1.e-6, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Vec read in does not match vector written out"); 85 86 PetscCall(VecDestroy(&input)); 87 PetscCall(VecDestroy(&gauss)); 88 PetscCall(DMDestroy(&da2D)); 89 PetscCall(PetscFinalize()); 90 return 0; 91 } 92 93 /*TEST 94 95 build: 96 requires: hdf5 !defined(PETSC_USE_CXXCOMPLEX) 97 98 test: 99 nsize: 4 100 101 test: 102 nsize: 4 103 suffix: 2 104 args: -viewer_hdf5_base_dimension2 105 output_file: output/ex10_1.out 106 107 test: 108 nsize: 4 109 suffix: 3 110 args: -viewer_hdf5_sp_output 111 output_file: output/ex10_1.out 112 113 TEST*/ 114