Actual source code: ex33.c


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

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

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

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

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

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

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

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

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

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

 85:   if (flg2) {
 86:     PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, &da2));
 87:   } else if (flg3) {
 88:     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));
 89:   } else {
 90:     PetscCall(DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da2));
 91:   }
 92:   PetscCall(DMSetFromOptions(da2));
 93:   PetscCall(DMSetUp(da2));

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

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

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

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