Actual source code: ex23.c
2: static char help[] = "Tests VecView()/VecLoad() for DMDA vectors (this tests DMDAGlobalToNatural()).\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
7: int main(int argc, char **argv)
8: {
9: PetscMPIInt size;
10: 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;
11: PetscBool flg2, flg3, native = PETSC_FALSE;
12: DMBoundaryType bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE, bz = DM_BOUNDARY_NONE;
13: DMDAStencilType stencil_type = DMDA_STENCIL_STAR;
14: DM da;
15: Vec global1, global2, global3, global4;
16: PetscScalar mone = -1.0;
17: PetscReal norm;
18: PetscViewer viewer;
19: PetscRandom rdm;
22: PetscInitialize(&argc, &argv, (char *)0, help);
23: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
24: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
25: PetscOptionsGetInt(NULL, NULL, "-P", &P, NULL);
26: PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL);
27: PetscOptionsGetInt(NULL, NULL, "-stencil_width", &stencil_width, NULL);
28: PetscOptionsGetInt(NULL, NULL, "-periodic", &pt, NULL);
29: PetscOptionsGetBool(NULL, NULL, "-native", &native, NULL);
30: if (pt == 1) bx = DM_BOUNDARY_PERIODIC;
31: if (pt == 2) by = DM_BOUNDARY_PERIODIC;
32: if (pt == 3) {
33: bx = DM_BOUNDARY_PERIODIC;
34: by = DM_BOUNDARY_PERIODIC;
35: }
36: if (pt == 4) bz = DM_BOUNDARY_PERIODIC;
38: PetscOptionsGetInt(NULL, NULL, "-stencil_type", &st, NULL);
39: stencil_type = (DMDAStencilType)st;
41: PetscOptionsHasName(NULL, NULL, "-one", &flg2);
42: PetscOptionsHasName(NULL, NULL, "-two", &flg2);
43: PetscOptionsHasName(NULL, NULL, "-three", &flg3);
44: if (flg2) {
45: DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, m, n, dof, stencil_width, 0, 0, &da);
46: } else if (flg3) {
47: DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stencil_type, M, N, P, m, n, p, dof, stencil_width, 0, 0, 0, &da);
48: } else {
49: DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da);
50: }
51: DMSetFromOptions(da);
52: DMSetUp(da);
54: DMCreateGlobalVector(da, &global1);
55: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
56: PetscRandomSetFromOptions(rdm);
57: DMCreateGlobalVector(da, &global2);
58: DMCreateGlobalVector(da, &global3);
59: DMCreateGlobalVector(da, &global4);
61: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_WRITE, &viewer);
62: if (native) PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
63: VecSetRandom(global1, rdm);
64: VecView(global1, viewer);
65: VecSetRandom(global3, rdm);
66: VecView(global3, viewer);
67: if (native) PetscViewerPopFormat(viewer);
68: PetscViewerDestroy(&viewer);
70: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer);
71: if (native) PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
72: VecLoad(global2, viewer);
73: VecLoad(global4, viewer);
74: if (native) PetscViewerPopFormat(viewer);
75: PetscViewerDestroy(&viewer);
77: if (native) {
78: Vec filenative;
79: PetscBool same;
80: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer);
81: DMDACreateNaturalVector(da, &filenative);
82: /* DMDA "natural" Vec does not commandeer VecLoad. The following load will only work when run on the same process
83: * layout, where as the standard VecView/VecLoad (using DMDA and not PETSC_VIEWER_NATIVE) can be read on a different
84: * number of processors. */
85: VecLoad(filenative, viewer);
86: VecEqual(global2, filenative, &same);
87: if (!same) {
88: PetscPrintf(PETSC_COMM_WORLD, "ex23: global vector does not match contents of file\n");
89: VecView(global2, 0);
90: VecView(filenative, 0);
91: }
92: PetscViewerDestroy(&viewer);
93: VecDestroy(&filenative);
94: }
96: VecAXPY(global2, mone, global1);
97: VecNorm(global2, NORM_MAX, &norm);
98: if (norm != 0.0) {
99: MPI_Comm_size(PETSC_COMM_WORLD, &size);
100: PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm);
101: PetscPrintf(PETSC_COMM_WORLD, " Number of processors %d\n", size);
102: PetscPrintf(PETSC_COMM_WORLD, " M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof);
103: PetscPrintf(PETSC_COMM_WORLD, " stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)bx);
104: PetscPrintf(PETSC_COMM_WORLD, " dimension %d\n", 1 + (int)flg2 + (int)flg3);
105: }
106: VecAXPY(global4, mone, global3);
107: VecNorm(global4, NORM_MAX, &norm);
108: if (norm != 0.0) {
109: MPI_Comm_size(PETSC_COMM_WORLD, &size);
110: PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm);
111: PetscPrintf(PETSC_COMM_WORLD, " Number of processors %d\n", size);
112: PetscPrintf(PETSC_COMM_WORLD, " M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof);
113: PetscPrintf(PETSC_COMM_WORLD, " stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)bx);
114: PetscPrintf(PETSC_COMM_WORLD, " dimension %d\n", 1 + (int)flg2 + (int)flg3);
115: }
117: PetscRandomDestroy(&rdm);
118: DMDestroy(&da);
119: VecDestroy(&global1);
120: VecDestroy(&global2);
121: VecDestroy(&global3);
122: VecDestroy(&global4);
123: PetscFinalize();
124: return 0;
125: }
127: /*TEST
129: test:
130: nsize: {{1 3}}
131: args: -one -dof {{1 2 3}} -stencil_type {{0 1}}
133: test:
134: suffix: 3
135: nsize: {{2 4}}
136: args: -two -dof {{1 3}} -stencil_type {{0 1}}
138: test:
139: suffix: 4
140: nsize: {{1 4}}
141: args: -three -dof {{2 3}} -stencil_type {{0 1}}
143: test:
144: suffix: 2
145: nsize: 2
146: args: -two -native
148: TEST*/