Actual source code: ex15.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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       dim         = 2;
     14:   PetscInt       numFields   = 1;
     15:   PetscInt       numBC       = 0;
     16:   PetscInt       numComp[1]  = {2};
     17:   PetscInt       numDof[3]   = {2, 0, 0};
     18:   PetscInt       bcFields[1] = {0};
     19:   IS             bcPoints[1] = {NULL};
     20:   PetscSection   section;
     21:   PetscReal      norm;

     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:   PetscOptionsRangeInt("-dim","the dimension of the problem","",2,&dim,NULL,1,3);
     31:   PetscOptionsEnd();

     33:   DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, &dm);
     34:   DMGetDimension(dm, &dim);
     35:   numDof[0] = dim;
     36:   DMSetNumFields(dm, numFields);
     37:   DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, &section);
     38:   DMSetLocalSection(dm, section);
     39:   PetscSectionDestroy(&section);
     40:   DMSetUseNatural(dm, PETSC_TRUE);
     41:   {
     42:     PetscPartitioner part;
     43:     DM               dmDist;

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

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

     59:   if (verbose) {
     60:     PetscInt size, bs;

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

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

     73:   if (verbose) {
     74:     PetscInt size, bs;

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

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

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

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

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

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

    156: TEST*/