Actual source code: ex15.c

  1: static char help[] = "An example of writing a global Vec from a DMPlex with HDF5 format.\n";

  3: #include <petscdmplex.h>
  4: #include <petscviewerhdf5.h>

  6: int main(int argc, char **argv)
  7: {
  8:   MPI_Comm       comm;
  9:   DM             dm;
 10:   Vec            v, nv, rv, coord;
 11:   PetscBool      test_read = PETSC_FALSE, verbose = PETSC_FALSE, flg;
 12:   PetscViewer    hdf5Viewer;
 13:   PetscInt       numFields   = 1;
 14:   PetscInt       numBC       = 0;
 15:   PetscInt       numComp[1]  = {2};
 16:   PetscInt       numDof[3]   = {2, 0, 0};
 17:   PetscInt       bcFields[1] = {0};
 18:   IS             bcPoints[1] = {NULL};
 19:   PetscSection   section;
 20:   PetscReal      norm;
 21:   PetscInt       dim;

 24:   PetscInitialize(&argc, &argv, (char *) 0, help);if (ierr) return ierr;
 25:   comm = PETSC_COMM_WORLD;

 27:   PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");
 28:   PetscOptionsBool("-test_read","Test reading from the HDF5 file","",PETSC_FALSE,&test_read,NULL);
 29:   PetscOptionsBool("-verbose","print the Vecs","",PETSC_FALSE,&verbose,NULL);
 30:   PetscOptionsEnd();

 32:   DMCreate(comm, &dm);
 33:   DMSetType(dm, DMPLEX);
 34:   DMSetFromOptions(dm);
 35:   DMViewFromOptions(dm, NULL, "-dm_view");
 36:   DMGetDimension(dm, &dim);
 37:   numDof[0] = dim;
 38:   DMSetNumFields(dm, numFields);
 39:   DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, &section);
 40:   DMSetLocalSection(dm, section);
 41:   PetscSectionDestroy(&section);
 42:   DMSetUseNatural(dm, PETSC_TRUE);
 43:   {
 44:     PetscPartitioner part;
 45:     DM               dmDist;

 47:     DMPlexGetPartitioner(dm,&part);
 48:     PetscPartitionerSetFromOptions(part);
 49:     DMPlexDistribute(dm, 0, NULL, &dmDist);
 50:     if (dmDist) {
 51:       DMDestroy(&dm);
 52:       dm   = dmDist;
 53:     }
 54:   }

 56:   DMCreateGlobalVector(dm, &v);
 57:   PetscObjectSetName((PetscObject) v, "V");
 58:   DMGetCoordinates(dm, &coord);
 59:   VecCopy(coord, v);

 61:   if (verbose) {
 62:     PetscInt size, bs;

 64:     VecGetSize(v, &size);
 65:     VecGetBlockSize(v, &bs);
 66:     PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d\tblock size=%d\n", size, bs);
 67:     VecView(v, PETSC_VIEWER_STDOUT_WORLD);
 68:   }

 70:   DMCreateGlobalVector(dm, &nv);
 71:   PetscObjectSetName((PetscObject) nv, "NV");
 72:   DMPlexGlobalToNaturalBegin(dm, v, nv);
 73:   DMPlexGlobalToNaturalEnd(dm, v, nv);

 75:   if (verbose) {
 76:     PetscInt size, bs;

 78:     VecGetSize(nv, &size);
 79:     VecGetBlockSize(nv, &bs);
 80:     PetscPrintf(PETSC_COMM_WORLD, "====  V in natural ordering. size==%d\tblock size=%d\n", size, bs);
 81:     VecView(nv, PETSC_VIEWER_STDOUT_WORLD);
 82:   }

 84:   VecViewFromOptions(v, NULL, "-global_vec_view");

 86:   if (test_read) {
 87:     DMCreateGlobalVector(dm, &rv);
 88:     PetscObjectSetName((PetscObject) rv, "V");
 89:     /* Test native read */
 90:     PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
 91:     PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE);
 92:     VecLoad(rv, hdf5Viewer);
 93:     PetscViewerPopFormat(hdf5Viewer);
 94:     PetscViewerDestroy(&hdf5Viewer);
 95:     if (verbose) {
 96:       PetscInt size, bs;

 98:       VecGetSize(rv, &size);
 99:       VecGetBlockSize(rv, &bs);
100:       PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
101:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
102:     }
103:     VecEqual(rv, v, &flg);
104:     if (flg) {
105:       PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal\n");
106:     } else {
107:       PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal\n\n");
108:       VecAXPY(rv, -1.0, v);
109:       VecNorm(rv, NORM_INFINITY, &norm);
110:       PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
111:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
112:     }
113:     /* Test raw read */
114:     PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
115:     VecLoad(rv, hdf5Viewer);
116:     PetscViewerDestroy(&hdf5Viewer);
117:     if (verbose) {
118:       PetscInt size, bs;

120:       VecGetSize(rv, &size);
121:       VecGetBlockSize(rv, &bs);
122:       PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
123:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
124:     }
125:     VecEqual(rv, nv, &flg);
126:     if (flg) {
127:       PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n");
128:     } else {
129:       PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n");
130:       VecAXPY(rv, -1.0, v);
131:       VecNorm(rv, NORM_INFINITY, &norm);
132:       PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
133:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
134:     }
135:     VecDestroy(&rv);
136:   }
137:   VecDestroy(&nv);
138:   VecDestroy(&v);
139:   DMDestroy(&dm);
140:   PetscFinalize();
141:   return ierr;
142: }

144: /*TEST
145:   build:
146:     requires: triangle hdf5
147:   test:
148:     suffix: 0
149:     requires: triangle hdf5
150:     nsize: 2
151:     args: -petscpartitioner_type simple -verbose -globaltonatural_sf_view
152:   test:
153:     suffix: 1
154:     requires: triangle hdf5
155:     nsize: 2
156:     args: -petscpartitioner_type simple -verbose -global_vec_view hdf5:V.h5:native -test_read

158: TEST*/