Actual source code: ex33.c

  1: static char help[] = "Tests VecView()/VecLoad() for DMDA vectors (this tests DMDAGlobalToNatural()).\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmda.h>
  5: #include <petscviewerhdf5.h>

  7: int main(int argc, char **argv)
  8: {
  9:   PetscMPIInt     rank, size;
 10:   PetscInt        N = 6, M = 8, P = 5, dof = 1;
 11:   PetscInt        stencil_width = 1, pt = 0, st = 0;
 12:   PetscBool       flg2, flg3, isbinary, mpiio;
 13:   DMBoundaryType  bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE, bz = DM_BOUNDARY_NONE;
 14:   DMDAStencilType stencil_type = DMDA_STENCIL_STAR;
 15:   DM              da, da2;
 16:   Vec             global1, global2;
 17:   PetscScalar     mone = -1.0;
 18:   PetscReal       norm;
 19:   PetscViewer     viewer;
 20:   PetscRandom     rdm;
 21: #if defined(PETSC_HAVE_HDF5)
 22:   PetscBool ishdf5;
 23: #endif

 25:   PetscFunctionBeginUser;
 26:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 27:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 28:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 30:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 31:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
 32:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-P", &P, NULL));
 33:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
 34:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &stencil_width, NULL));
 35:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-periodic", &pt, NULL));
 36:   if (pt == 1) bx = DM_BOUNDARY_PERIODIC;
 37:   if (pt == 2) by = DM_BOUNDARY_PERIODIC;
 38:   if (pt == 4) {
 39:     bx = DM_BOUNDARY_PERIODIC;
 40:     by = DM_BOUNDARY_PERIODIC;
 41:   }

 43:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_type", &st, NULL));
 44:   stencil_type = (DMDAStencilType)st;

 46:   PetscCall(PetscOptionsHasName(NULL, NULL, "-oned", &flg2));
 47:   PetscCall(PetscOptionsHasName(NULL, NULL, "-twod", &flg2));
 48:   PetscCall(PetscOptionsHasName(NULL, NULL, "-threed", &flg3));

 50:   PetscCall(PetscOptionsHasName(NULL, NULL, "-binary", &isbinary));
 51: #if defined(PETSC_HAVE_HDF5)
 52:   PetscCall(PetscOptionsHasName(NULL, NULL, "-hdf5", &ishdf5));
 53: #endif
 54:   PetscCall(PetscOptionsHasName(NULL, NULL, "-mpiio", &mpiio));
 55:   if (flg2) {
 56:     PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, &da));
 57:   } else if (flg3) {
 58:     PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stencil_type, M, N, P, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, 0, &da));
 59:   } else {
 60:     PetscCall(DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da));
 61:   }
 62:   PetscCall(DMSetFromOptions(da));
 63:   PetscCall(DMSetUp(da));

 65:   PetscCall(DMCreateGlobalVector(da, &global1));
 66:   PetscCall(PetscObjectSetName((PetscObject)global1, "Test_Vec"));
 67:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
 68:   PetscCall(PetscRandomSetFromOptions(rdm));
 69:   PetscCall(VecSetRandom(global1, rdm));
 70:   if (isbinary) {
 71:     if (mpiio) PetscCall(PetscOptionsSetValue(NULL, "-viewer_binary_mpiio", ""));
 72:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_WRITE, &viewer));
 73: #if defined(PETSC_HAVE_HDF5)
 74:   } else if (ishdf5) {
 75:     PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "temp", FILE_MODE_WRITE, &viewer));
 76: #endif
 77:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid Viewer : Run with -binary or -hdf5 option");
 78:   PetscCall(VecView(global1, viewer));
 79:   PetscCall(PetscViewerDestroy(&viewer));

 81:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Global vector written to temp file is \n"));
 82:   PetscCall(VecView(global1, PETSC_VIEWER_STDOUT_WORLD));

 84:   if (flg2) {
 85:     PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, &da2));
 86:   } else if (flg3) {
 87:     PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stencil_type, M, N, P, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, 0, &da2));
 88:   } else {
 89:     PetscCall(DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da2));
 90:   }
 91:   PetscCall(DMSetFromOptions(da2));
 92:   PetscCall(DMSetUp(da2));

 94:   if (isbinary) {
 95:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer));
 96: #if defined(PETSC_HAVE_HDF5)
 97:   } else if (ishdf5) {
 98:     PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer));
 99: #endif
100:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid Viewer : Run with -binary or -hdf5 option");

102:   PetscCall(DMCreateGlobalVector(da2, &global2));
103:   PetscCall(PetscObjectSetName((PetscObject)global2, "Test_Vec"));
104:   PetscCall(VecLoad(global2, viewer));
105:   PetscCall(PetscViewerDestroy(&viewer));

107:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Global vector read from temp file is \n"));
108:   PetscCall(VecView(global2, PETSC_VIEWER_STDOUT_WORLD));
109:   PetscCall(VecAXPY(global2, mone, global1));
110:   PetscCall(VecNorm(global2, NORM_MAX, &norm));
111:   if (norm != 0.0) {
112:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm));
113:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Number of processors %d\n", size));
114:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof));
115:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)pt));
116:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  dimension %d\n", 1 + (int)flg2 + (int)flg3));
117:   }

119:   PetscCall(PetscRandomDestroy(&rdm));
120:   PetscCall(DMDestroy(&da));
121:   PetscCall(DMDestroy(&da2));
122:   PetscCall(VecDestroy(&global1));
123:   PetscCall(VecDestroy(&global2));
124:   PetscCall(PetscFinalize());
125:   return 0;
126: }