Actual source code: ex7.c
1: static char help[] = "Test DMStag 3d periodic and ghosted boundary conditions\n\n";
2: #include <petscdm.h>
3: #include <petscdmstag.h>
5: int main(int argc, char **argv)
6: {
7: DM dm;
8: Vec vec, vecLocal1, vecLocal2;
9: PetscScalar *a, ****a1, ****a2, expected;
10: PetscInt startx, starty, startz, nx, ny, nz, i, j, k, d, is, js, ks, dof0, dof1, dof2, dof3, dofTotal, stencilWidth, Nx, Ny, Nz;
11: DMBoundaryType boundaryTypex, boundaryTypey, boundaryTypez;
12: PetscMPIInt rank;
15: PetscInitialize(&argc, &argv, (char *)0, help);
16: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
17: dof0 = 1;
18: dof1 = 1;
19: dof2 = 1;
20: dof3 = 1;
21: stencilWidth = 2;
22: DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, 4, 4, 4, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, dof3, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, NULL, &dm);
23: DMSetFromOptions(dm);
24: DMSetUp(dm);
25: DMStagGetDOF(dm, &dof0, &dof1, &dof2, &dof3);
26: dofTotal = dof0 + 3 * dof1 + 3 * dof2 + dof3;
27: DMStagGetStencilWidth(dm, &stencilWidth);
29: DMCreateLocalVector(dm, &vecLocal1);
30: VecDuplicate(vecLocal1, &vecLocal2);
32: DMCreateGlobalVector(dm, &vec);
33: VecSet(vec, 1.0);
34: VecSet(vecLocal1, 0.0);
35: DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1);
36: DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1);
38: DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL);
39: DMStagVecGetArrayRead(dm, vecLocal1, &a1);
40: DMStagVecGetArray(dm, vecLocal2, &a2);
41: for (k = startz; k < startz + nz; ++k) {
42: for (j = starty; j < starty + ny; ++j) {
43: for (i = startx; i < startx + nx; ++i) {
44: for (d = 0; d < dofTotal; ++d) {
45: if (a1[k][j][i][d] != 1.0) PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a1[k][j][i][d]), 1.0);
46: a2[k][j][i][d] = 0.0;
47: for (ks = -stencilWidth; ks <= stencilWidth; ++ks) {
48: for (js = -stencilWidth; js <= stencilWidth; ++js) {
49: for (is = -stencilWidth; is <= stencilWidth; ++is) a2[k][j][i][d] += a1[k + ks][j + js][i + is][d];
50: }
51: }
52: }
53: }
54: }
55: }
56: DMStagVecRestoreArrayRead(dm, vecLocal1, &a1);
57: DMStagVecRestoreArray(dm, vecLocal2, &a2);
59: DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec);
60: DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec);
62: /* For the all-periodic case, all values are the same . Otherwise, just check the local version */
63: DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, &boundaryTypez);
64: if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
65: VecGetArray(vec, &a);
66: expected = 1.0;
67: for (d = 0; d < 3; ++d) expected *= (2 * stencilWidth + 1);
68: for (i = 0; i < nz * ny * nx * dofTotal; ++i) {
69: if (a[i] != expected) PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected));
70: }
71: VecRestoreArray(vec, &a);
72: } else {
73: DMStagVecGetArrayRead(dm, vecLocal2, &a2);
74: DMStagGetGlobalSizes(dm, &Nx, &Ny, &Nz);
76: for (k = startz; k < startz + nz; ++k) {
77: for (j = starty; j < starty + ny; ++j) {
78: for (i = startx; i < startx + nx; ++i) {
79: PetscInt dd, extra[3];
80: PetscBool bnd[3];
81: bnd[0] = (PetscBool)((i == 0 || i == Nx - 1) && boundaryTypex != DM_BOUNDARY_PERIODIC);
82: bnd[1] = (PetscBool)((j == 0 || j == Ny - 1) && boundaryTypey != DM_BOUNDARY_PERIODIC);
83: bnd[2] = (PetscBool)((k == 0 || k == Nz - 1) && boundaryTypez != DM_BOUNDARY_PERIODIC);
84: extra[0] = i == Nx - 1 && boundaryTypex != DM_BOUNDARY_PERIODIC ? 1 : 0;
85: extra[1] = j == Ny - 1 && boundaryTypey != DM_BOUNDARY_PERIODIC ? 1 : 0;
86: extra[2] = k == Nz - 1 && boundaryTypez != DM_BOUNDARY_PERIODIC ? 1 : 0;
87: { /* vertices */
88: PetscScalar expected = 1.0;
89: for (dd = 0; dd < 3; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
90: for (d = 0; d < dof0; ++d) {
91: if (a2[k][j][i][d] != expected) {
92: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
93: }
94: }
95: }
96: { /* back down edges */
97: PetscScalar expected = ((bnd[0] ? 1 : 2) * stencilWidth + 1);
98: for (dd = 1; dd < 3; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
99: for (d = dof0; d < dof0 + dof1; ++d) {
100: if (a2[k][j][i][d] != expected) {
101: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
102: }
103: }
104: }
105: { /* back left edges */
106: PetscScalar expected = ((bnd[1] ? 1 : 2) * stencilWidth + 1);
107: for (dd = 0; dd < 3; dd += 2) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
108: for (d = dof0 + dof1; d < dof0 + 2 * dof1; ++d) {
109: if (a2[k][j][i][d] != expected) {
110: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
111: }
112: }
113: }
114: { /* back faces */
115: PetscScalar expected = (bnd[2] ? stencilWidth + 1 + extra[2] : 2 * stencilWidth + 1);
116: for (dd = 0; dd < 2; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
117: for (d = dof0 + 2 * dof1; d < dof0 + 2 * dof1 + dof2; ++d) {
118: if (a2[k][j][i][d] != expected) {
119: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
120: }
121: }
122: }
123: { /* down left edges */
124: PetscScalar expected = ((bnd[2] ? 1 : 2) * stencilWidth + 1);
125: for (dd = 0; dd < 2; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
126: for (d = dof0 + 2 * dof1 + dof2; d < dof0 + 3 * dof1 + dof2; ++d) {
127: if (a2[k][j][i][d] != expected) {
128: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
129: }
130: }
131: }
132: { /* down faces */
133: PetscScalar expected = (bnd[1] ? stencilWidth + 1 + extra[1] : 2 * stencilWidth + 1);
134: for (dd = 0; dd < 3; dd += 2) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
135: for (d = dof0 + 3 * dof1 + dof2; d < dof0 + 3 * dof1 + 2 * dof2; ++d) {
136: if (a2[k][j][i][d] != expected) {
137: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
138: }
139: }
140: }
141: { /* left faces */
142: PetscScalar expected = (bnd[0] ? stencilWidth + 1 + extra[0] : 2 * stencilWidth + 1);
143: for (dd = 1; dd < 3; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
144: for (d = dof0 + 3 * dof1 + 2 * dof2; d < dof0 + 3 * dof1 + 3 * dof2; ++d) {
145: if (a2[k][j][i][d] != expected) {
146: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
147: }
148: }
149: }
150: { /* elements */
151: PetscScalar expected = 1.0;
152: for (dd = 0; dd < 3; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
153: for (d = dofTotal - dof3; d < dofTotal; ++d) {
154: if (a2[k][j][i][d] != expected) {
155: PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected));
156: }
157: }
158: }
159: }
160: }
161: }
162: DMStagVecRestoreArrayRead(dm, vecLocal2, &a2);
163: }
165: VecDestroy(&vec);
166: VecDestroy(&vecLocal1);
167: VecDestroy(&vecLocal2);
168: DMDestroy(&dm);
169: PetscFinalize();
170: return 0;
171: }
173: /*TEST
175: test:
176: suffix: 1
177: nsize: 8
178: args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_stencil_width 1 -stag_dof_3 2 -stag_grid_z 3
180: test:
181: suffix: 2
182: nsize: 8
183: args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_2 2 -stag_grid_y 5
185: test:
186: suffix: 3
187: nsize: 12
188: args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 2 -stag_grid_x 6
190: test:
191: suffix: 4
192: nsize: 12
193: args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 0 -stag_dof_1 0 -stag_dof_2 0 -stag_grid_x 4 -stag_boundary_type_x ghosted -stag_boundary_type_y ghosted -stag_boundary_type_z ghosted -stag_stencil_width 1
195: test:
196: suffix: 5
197: nsize: 12
198: args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 0 -stag_dof_1 0 -stag_dof_2 0 -stag_grid_x 4 -stag_boundary_type_x ghosted -stag_boundary_type_z ghosted -stag_stencil_width 1
200: test:
201: suffix: 6
202: nsize: 8
203: args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_dof_3 2 -stag_boundary_type_y ghosted -stag_boundary_type_z ghosted -stag_stencil_width 1
204: TEST*/