Actual source code: ex23.c

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

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

  6: int main(int argc, char **argv)
  7: {
  8:   PetscMPIInt     size;
  9:   PetscInt        N = 6, m = PETSC_DECIDE, n = PETSC_DECIDE, p = PETSC_DECIDE, M = 8, dof = 1, stencil_width = 1, P = 5, pt = 0, st = 0;
 10:   PetscBool       flg2, flg3, native = PETSC_FALSE;
 11:   DMBoundaryType  bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE, bz = DM_BOUNDARY_NONE;
 12:   DMDAStencilType stencil_type = DMDA_STENCIL_STAR;
 13:   DM              da;
 14:   Vec             global1, global2, global3, global4;
 15:   PetscScalar     mone = -1.0;
 16:   PetscReal       norm;
 17:   PetscViewer     viewer;
 18:   PetscRandom     rdm;

 20:   PetscFunctionBeginUser;
 21:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 22:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 23:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
 24:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-P", &P, NULL));
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
 26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &stencil_width, NULL));
 27:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-periodic", &pt, NULL));
 28:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-native", &native, NULL));
 29:   if (pt == 1) bx = DM_BOUNDARY_PERIODIC;
 30:   if (pt == 2) by = DM_BOUNDARY_PERIODIC;
 31:   if (pt == 3) {
 32:     bx = DM_BOUNDARY_PERIODIC;
 33:     by = DM_BOUNDARY_PERIODIC;
 34:   }
 35:   if (pt == 4) bz = DM_BOUNDARY_PERIODIC;

 37:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_type", &st, NULL));
 38:   stencil_type = (DMDAStencilType)st;

 40:   PetscCall(PetscOptionsHasName(NULL, NULL, "-one", &flg2));
 41:   PetscCall(PetscOptionsHasName(NULL, NULL, "-two", &flg2));
 42:   PetscCall(PetscOptionsHasName(NULL, NULL, "-three", &flg3));
 43:   if (flg2) {
 44:     PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, m, n, dof, stencil_width, 0, 0, &da));
 45:   } else if (flg3) {
 46:     PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stencil_type, M, N, P, m, n, p, dof, stencil_width, 0, 0, 0, &da));
 47:   } else {
 48:     PetscCall(DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da));
 49:   }
 50:   PetscCall(DMSetFromOptions(da));
 51:   PetscCall(DMSetUp(da));

 53:   PetscCall(DMCreateGlobalVector(da, &global1));
 54:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
 55:   PetscCall(PetscRandomSetFromOptions(rdm));
 56:   PetscCall(DMCreateGlobalVector(da, &global2));
 57:   PetscCall(DMCreateGlobalVector(da, &global3));
 58:   PetscCall(DMCreateGlobalVector(da, &global4));

 60:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_WRITE, &viewer));
 61:   if (native) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
 62:   PetscCall(VecSetRandom(global1, rdm));
 63:   PetscCall(VecView(global1, viewer));
 64:   PetscCall(VecSetRandom(global3, rdm));
 65:   PetscCall(VecView(global3, viewer));
 66:   if (native) PetscCall(PetscViewerPopFormat(viewer));
 67:   PetscCall(PetscViewerDestroy(&viewer));

 69:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer));
 70:   if (native) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
 71:   PetscCall(VecLoad(global2, viewer));
 72:   PetscCall(VecLoad(global4, viewer));
 73:   if (native) PetscCall(PetscViewerPopFormat(viewer));
 74:   PetscCall(PetscViewerDestroy(&viewer));

 76:   if (native) {
 77:     Vec       filenative;
 78:     PetscBool same;
 79:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer));
 80:     PetscCall(DMDACreateNaturalVector(da, &filenative));
 81:     /* DMDA "natural" Vec does not commandeer VecLoad.  The following load will only work when run on the same process
 82:      * layout, where as the standard VecView/VecLoad (using DMDA and not PETSC_VIEWER_NATIVE) can be read on a different
 83:      * number of processors. */
 84:     PetscCall(VecLoad(filenative, viewer));
 85:     PetscCall(VecEqual(global2, filenative, &same));
 86:     if (!same) {
 87:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: global vector does not match contents of file\n"));
 88:       PetscCall(VecView(global2, 0));
 89:       PetscCall(VecView(filenative, 0));
 90:     }
 91:     PetscCall(PetscViewerDestroy(&viewer));
 92:     PetscCall(VecDestroy(&filenative));
 93:   }

 95:   PetscCall(VecAXPY(global2, mone, global1));
 96:   PetscCall(VecNorm(global2, NORM_MAX, &norm));
 97:   if (norm != 0.0) {
 98:     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 99:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm));
100:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Number of processors %d\n", size));
101:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof));
102:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)bx));
103:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  dimension %d\n", 1 + (int)flg2 + (int)flg3));
104:   }
105:   PetscCall(VecAXPY(global4, mone, global3));
106:   PetscCall(VecNorm(global4, NORM_MAX, &norm));
107:   if (norm != 0.0) {
108:     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
109:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm));
110:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Number of processors %d\n", size));
111:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof));
112:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)bx));
113:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  dimension %d\n", 1 + (int)flg2 + (int)flg3));
114:   }

116:   PetscCall(PetscRandomDestroy(&rdm));
117:   PetscCall(DMDestroy(&da));
118:   PetscCall(VecDestroy(&global1));
119:   PetscCall(VecDestroy(&global2));
120:   PetscCall(VecDestroy(&global3));
121:   PetscCall(VecDestroy(&global4));
122:   PetscCall(PetscFinalize());
123:   return 0;
124: }

126: /*TEST

128:    test:
129:       nsize: {{1 3}}
130:       args: -one -dof {{1 2 3}} -stencil_type {{0 1}}

132:    test:
133:       suffix: 3
134:       nsize: {{2 4}}
135:       args: -two -dof {{1 3}} -stencil_type {{0 1}}

137:    test:
138:       suffix: 4
139:       nsize: {{1 4}}
140:       args: -three -dof {{2 3}} -stencil_type {{0 1}}

142:    test:
143:       suffix: 2
144:       nsize: 2
145:       args: -two -native

147: TEST*/