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;

 23:   PetscFunctionBeginUser;
 24:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 25:   comm = PETSC_COMM_WORLD;

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

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

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

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

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

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

 71:   PetscCall(DMPlexCreateNaturalVector(dm, &nv));
 72:   PetscCall(PetscObjectSetName((PetscObject)nv, "NV"));
 73:   PetscCall(DMPlexGlobalToNaturalBegin(dm, v, nv));
 74:   PetscCall(DMPlexGlobalToNaturalEnd(dm, v, nv));

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

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

 85:   PetscCall(VecViewFromOptions(v, NULL, "-global_vec_view"));

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

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

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

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

159: TEST*/