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