1 static char help[] = "An example of writing a global Vec from a DMPlex with HDF5 format.\n"; 2 3 #include <petscdmplex.h> 4 #include <petscviewerhdf5.h> 5 6 int main(int argc, char **argv) { 7 MPI_Comm comm; 8 DM dm; 9 Vec v, nv, rv, coord; 10 PetscBool test_read = PETSC_FALSE, verbose = PETSC_FALSE, flg; 11 PetscViewer hdf5Viewer; 12 PetscInt numFields = 1; 13 PetscInt numBC = 0; 14 PetscInt numComp[1] = {2}; 15 PetscInt numDof[3] = {2, 0, 0}; 16 PetscInt bcFields[1] = {0}; 17 IS bcPoints[1] = {NULL}; 18 PetscSection section; 19 PetscReal norm; 20 PetscInt dim; 21 22 PetscFunctionBeginUser; 23 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 24 comm = PETSC_COMM_WORLD; 25 26 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Test Options", "none"); 27 PetscCall(PetscOptionsBool("-test_read", "Test reading from the HDF5 file", "", PETSC_FALSE, &test_read, NULL)); 28 PetscCall(PetscOptionsBool("-verbose", "print the Vecs", "", PETSC_FALSE, &verbose, NULL)); 29 PetscOptionsEnd(); 30 31 PetscCall(DMCreate(comm, &dm)); 32 PetscCall(DMSetType(dm, DMPLEX)); 33 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE)); 34 PetscCall(DMSetFromOptions(dm)); 35 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 36 PetscCall(DMGetDimension(dm, &dim)); 37 numDof[0] = dim; 38 PetscCall(DMSetNumFields(dm, numFields)); 39 PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, §ion)); 40 PetscCall(DMSetLocalSection(dm, section)); 41 PetscCall(PetscSectionDestroy(§ion)); 42 PetscCall(DMSetUseNatural(dm, PETSC_TRUE)); 43 { 44 PetscPartitioner part; 45 DM dmDist; 46 47 PetscCall(DMPlexGetPartitioner(dm, &part)); 48 PetscCall(PetscPartitionerSetFromOptions(part)); 49 PetscCall(DMPlexDistribute(dm, 0, NULL, &dmDist)); 50 if (dmDist) { 51 PetscCall(DMDestroy(&dm)); 52 dm = dmDist; 53 } 54 } 55 56 PetscCall(DMCreateGlobalVector(dm, &v)); 57 PetscCall(PetscObjectSetName((PetscObject)v, "V")); 58 PetscCall(DMGetCoordinates(dm, &coord)); 59 PetscCall(VecCopy(coord, v)); 60 61 if (verbose) { 62 PetscInt size, bs; 63 64 PetscCall(VecGetSize(v, &size)); 65 PetscCall(VecGetBlockSize(v, &bs)); 66 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d\tblock size=%d\n", size, bs)); 67 PetscCall(VecView(v, PETSC_VIEWER_STDOUT_WORLD)); 68 } 69 70 PetscCall(DMPlexCreateNaturalVector(dm, &nv)); 71 PetscCall(PetscObjectSetName((PetscObject)nv, "NV")); 72 PetscCall(DMPlexGlobalToNaturalBegin(dm, v, nv)); 73 PetscCall(DMPlexGlobalToNaturalEnd(dm, v, nv)); 74 75 if (verbose) { 76 PetscInt size, bs; 77 78 PetscCall(VecGetSize(nv, &size)); 79 PetscCall(VecGetBlockSize(nv, &bs)); 80 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== V in natural ordering. size==%d\tblock size=%d\n", size, bs)); 81 PetscCall(VecView(nv, PETSC_VIEWER_STDOUT_WORLD)); 82 } 83 84 PetscCall(VecViewFromOptions(v, NULL, "-global_vec_view")); 85 86 if (0 && test_read) { 87 PetscCall(DMCreateGlobalVector(dm, &rv)); 88 PetscCall(PetscObjectSetName((PetscObject)rv, "V")); 89 /* Test native read */ 90 PetscCall(PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer)); 91 PetscCall(PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE)); 92 PetscCall(VecLoad(rv, hdf5Viewer)); 93 PetscCall(PetscViewerPopFormat(hdf5Viewer)); 94 PetscCall(PetscViewerDestroy(&hdf5Viewer)); 95 if (verbose) { 96 PetscInt size, bs; 97 98 PetscCall(VecGetSize(rv, &size)); 99 PetscCall(VecGetBlockSize(rv, &bs)); 100 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs)); 101 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 102 } 103 PetscCall(VecEqual(rv, v, &flg)); 104 if (flg) { 105 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal\n")); 106 } else { 107 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal\n\n")); 108 PetscCall(VecAXPY(rv, -1.0, v)); 109 PetscCall(VecNorm(rv, NORM_INFINITY, &norm)); 110 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double)norm)); 111 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 112 } 113 /* Test raw read */ 114 PetscCall(PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer)); 115 PetscCall(VecLoad(rv, hdf5Viewer)); 116 PetscCall(PetscViewerDestroy(&hdf5Viewer)); 117 if (verbose) { 118 PetscInt size, bs; 119 120 PetscCall(VecGetSize(rv, &size)); 121 PetscCall(VecGetBlockSize(rv, &bs)); 122 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs)); 123 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 124 } 125 PetscCall(VecEqual(rv, nv, &flg)); 126 if (flg) { 127 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n")); 128 } else { 129 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n")); 130 PetscCall(VecAXPY(rv, -1.0, v)); 131 PetscCall(VecNorm(rv, NORM_INFINITY, &norm)); 132 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double)norm)); 133 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 134 } 135 PetscCall(VecDestroy(&rv)); 136 } 137 PetscCall(VecDestroy(&nv)); 138 PetscCall(VecDestroy(&v)); 139 PetscCall(DMDestroy(&dm)); 140 PetscCall(PetscFinalize()); 141 return 0; 142 } 143 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 TODO: broken 155 requires: triangle hdf5 156 nsize: 2 157 args: -petscpartitioner_type simple -verbose -global_vec_view hdf5:V.h5:native -test_read 158 159 TEST*/ 160