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);
 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:   DMPlexDistributeSetDefault(dm, PETSC_FALSE);
 35:   DMSetFromOptions(dm);
 36:   DMViewFromOptions(dm, NULL, "-dm_view");
 37:   DMGetDimension(dm, &dim);
 38:   numDof[0] = dim;
 39:   DMSetNumFields(dm, numFields);
 40:   DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, &section);
 41:   DMSetLocalSection(dm, section);
 42:   PetscSectionDestroy(&section);
 43:   DMSetUseNatural(dm, PETSC_TRUE);
 44:   {
 45:     PetscPartitioner part;
 46:     DM               dmDist;

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

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

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

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

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

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

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

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

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

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

121:       VecGetSize(rv, &size);
122:       VecGetBlockSize(rv, &bs);
123:       PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
124:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
125:     }
126:     VecEqual(rv, nv, &flg);
127:     if (flg) {
128:       PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n");
129:     } else {
130:       PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n");
131:       VecAXPY(rv, -1.0, v);
132:       VecNorm(rv, NORM_INFINITY, &norm);
133:       PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double)norm);
134:       VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
135:     }
136:     VecDestroy(&rv);
137:   }
138:   VecDestroy(&nv);
139:   VecDestroy(&v);
140:   DMDestroy(&dm);
141:   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:     TODO: broken
156:     requires: triangle hdf5
157:     nsize: 2
158:     args: -petscpartitioner_type simple -verbose -global_vec_view hdf5:V.h5:native -test_read

160: TEST*/