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 { 8 MPI_Comm comm; 9 DM dm, ddm; 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; 22 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 46 PetscCall(DMPlexGetPartitioner(dm,&part)); 47 PetscCall(PetscPartitionerSetFromOptions(part)); 48 PetscCall(DMPlexDistribute(dm, 0, NULL, &ddm)); 49 } 50 51 PetscCall(DMCreateGlobalVector(ddm, &v)); 52 PetscCall(PetscObjectSetName((PetscObject) v, "V")); 53 PetscCall(DMGetCoordinates(ddm, &coord)); 54 PetscCall(VecCopy(coord, v)); 55 56 if (verbose) { 57 PetscInt size, bs; 58 59 PetscCall(VecGetSize(v, &size)); 60 PetscCall(VecGetBlockSize(v, &bs)); 61 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d\tblock size=%d\n", size, bs)); 62 PetscCall(VecView(v, PETSC_VIEWER_STDOUT_WORLD)); 63 } 64 65 PetscCall(DMCreateGlobalVector(dm, &nv)); 66 PetscCall(PetscObjectSetName((PetscObject) nv, "NV")); 67 PetscCall(DMPlexGlobalToNaturalBegin(ddm, v, nv)); 68 PetscCall(DMPlexGlobalToNaturalEnd(ddm, v, nv)); 69 70 if (verbose) { 71 PetscInt size, bs; 72 73 PetscCall(VecGetSize(nv, &size)); 74 PetscCall(VecGetBlockSize(nv, &bs)); 75 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== V in natural ordering. size==%d\tblock size=%d\n", size, bs)); 76 PetscCall(VecView(nv, PETSC_VIEWER_STDOUT_WORLD)); 77 } 78 79 PetscCall(VecViewFromOptions(v, NULL, "-global_vec_view")); 80 81 if (0 && test_read) { 82 PetscCall(DMCreateGlobalVector(ddm, &rv)); 83 PetscCall(PetscObjectSetName((PetscObject) rv, "V")); 84 /* Test native read */ 85 PetscCall(PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer)); 86 PetscCall(PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE)); 87 PetscCall(VecLoad(rv, hdf5Viewer)); 88 PetscCall(PetscViewerPopFormat(hdf5Viewer)); 89 PetscCall(PetscViewerDestroy(&hdf5Viewer)); 90 if (verbose) { 91 PetscInt size, bs; 92 93 PetscCall(VecGetSize(rv, &size)); 94 PetscCall(VecGetBlockSize(rv, &bs)); 95 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs)); 96 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 97 } 98 PetscCall(VecEqual(rv, v, &flg)); 99 if (flg) { 100 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal\n")); 101 } else { 102 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal\n\n")); 103 PetscCall(VecAXPY(rv, -1.0, v)); 104 PetscCall(VecNorm(rv, NORM_INFINITY, &norm)); 105 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm)); 106 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 107 } 108 /* Test raw read */ 109 PetscCall(PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer)); 110 PetscCall(VecLoad(rv, hdf5Viewer)); 111 PetscCall(PetscViewerDestroy(&hdf5Viewer)); 112 if (verbose) { 113 PetscInt size, bs; 114 115 PetscCall(VecGetSize(rv, &size)); 116 PetscCall(VecGetBlockSize(rv, &bs)); 117 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs)); 118 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 119 } 120 PetscCall(VecEqual(rv, nv, &flg)); 121 if (flg) { 122 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n")); 123 } else { 124 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n")); 125 PetscCall(VecAXPY(rv, -1.0, v)); 126 PetscCall(VecNorm(rv, NORM_INFINITY, &norm)); 127 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm)); 128 PetscCall(VecView(rv, PETSC_VIEWER_STDOUT_WORLD)); 129 } 130 PetscCall(VecDestroy(&rv)); 131 } 132 PetscCall(VecDestroy(&nv)); 133 PetscCall(VecDestroy(&v)); 134 PetscCall(DMDestroy(&ddm)); 135 PetscCall(DMDestroy(&dm)); 136 PetscCall(PetscFinalize()); 137 return 0; 138 } 139 140 /*TEST 141 #build: 142 # requires: triangle hdf5 143 #test: 144 # suffix: 0 145 # requires: triangle hdf5 146 # nsize: 2 147 # args: -petscpartitioner_type simple -verbose -globaltonatural_sf_view 148 #test: 149 # suffix: 1 150 # requires: triangle hdf5 151 # nsize: 2 152 # args: -petscpartitioner_type simple -verbose -global_vec_view hdf5:V.h5:native -test_read 153 154 TEST*/ 155