Actual source code: ex10.c
1: static char help[] = "Test DMStag 2d 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, nx, ny, i, j, d, is, js, dof0, dof1, dof2, dofTotal, stencilWidth, Nx, Ny;
11: DMBoundaryType boundaryTypex, boundaryTypey;
12: PetscMPIInt rank;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
16: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
17: dof0 = 1;
18: dof1 = 1;
19: dof2 = 1;
20: stencilWidth = 2;
21: PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, 4, 4, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, &dm));
22: PetscCall(DMSetFromOptions(dm));
23: PetscCall(DMSetUp(dm));
24: PetscCall(DMStagGetDOF(dm, &dof0, &dof1, &dof2, NULL));
25: dofTotal = dof0 + 2 * dof1 + dof2;
26: PetscCall(DMStagGetStencilWidth(dm, &stencilWidth));
28: PetscCall(DMCreateLocalVector(dm, &vecLocal1));
29: PetscCall(VecDuplicate(vecLocal1, &vecLocal2));
31: PetscCall(DMCreateGlobalVector(dm, &vec));
32: PetscCall(VecSet(vec, 1.0));
33: PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1));
34: PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1));
36: PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
37: PetscCall(DMStagVecGetArrayRead(dm, vecLocal1, &a1));
38: PetscCall(DMStagVecGetArray(dm, vecLocal2, &a2));
39: for (j = starty; j < starty + ny; ++j) {
40: for (i = startx; i < startx + nx; ++i) {
41: for (d = 0; d < dofTotal; ++d) {
42: if (a1[j][i][d] != 1.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a1[j][i][d]), 1.0));
43: a2[j][i][d] = 0.0;
44: for (js = -stencilWidth; js <= stencilWidth; ++js) {
45: for (is = -stencilWidth; is <= stencilWidth; ++is) a2[j][i][d] += a1[j + js][i + is][d];
46: }
47: }
48: }
49: }
50: PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal1, &a1));
51: PetscCall(DMStagVecRestoreArray(dm, vecLocal2, &a2));
53: PetscCall(DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec));
54: PetscCall(DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec));
56: /* For the all-periodic case, all values are the same . Otherwise, just check the local version */
57: PetscCall(DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, NULL));
58: if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {
59: PetscCall(VecGetArray(vec, &a));
60: expected = 1.0;
61: for (d = 0; d < 2; ++d) expected *= (2 * stencilWidth + 1);
62: for (i = 0; i < ny * nx * dofTotal; ++i) {
63: if (a[i] != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected)));
64: }
65: PetscCall(VecRestoreArray(vec, &a));
66: } else {
67: PetscCall(DMStagVecGetArrayRead(dm, vecLocal2, &a2));
68: PetscCall(DMStagGetGlobalSizes(dm, &Nx, &Ny, NULL));
69: PetscCheck(stencilWidth <= 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Non-periodic check implemented assuming stencilWidth = 1");
70: for (j = starty; j < starty + ny; ++j) {
71: for (i = startx; i < startx + nx; ++i) {
72: PetscInt dd, extra[2];
73: PetscBool bnd[2];
74: bnd[0] = (PetscBool)((i == 0 || i == Nx - 1) && boundaryTypex != DM_BOUNDARY_PERIODIC);
75: bnd[1] = (PetscBool)((j == 0 || j == Ny - 1) && boundaryTypey != DM_BOUNDARY_PERIODIC);
76: extra[0] = i == Nx - 1 && boundaryTypex != DM_BOUNDARY_PERIODIC ? 1 : 0;
77: extra[1] = j == Ny - 1 && boundaryTypey != DM_BOUNDARY_PERIODIC ? 1 : 0;
78: { /* vertices */
79: PetscScalar expected = 1.0;
80: for (dd = 0; dd < 2; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
81: for (d = 0; d < dof0; ++d) {
82: if (a2[j][i][d] != expected) {
83: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
84: }
85: }
86: }
87: { /* down edges */
88: PetscScalar expected = (bnd[1] ? stencilWidth + 1 + extra[1] : 2 * stencilWidth + 1);
89: expected *= ((bnd[0] ? 1 : 2) * stencilWidth + 1);
90: for (d = dof0; d < dof0 + dof1; ++d) {
91: if (a2[j][i][d] != expected) {
92: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
93: }
94: }
95: }
96: { /* left edges */
97: PetscScalar expected = (bnd[0] ? stencilWidth + 1 + extra[0] : 2 * stencilWidth + 1);
98: expected *= ((bnd[1] ? 1 : 2) * stencilWidth + 1);
99: for (d = dof0 + dof1; d < dof0 + 2 * dof1; ++d) {
100: if (a2[j][i][d] != expected) {
101: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
102: }
103: }
104: }
105: { /* elements */
106: PetscScalar expected = 1.0;
107: for (dd = 0; dd < 2; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
108: for (d = dofTotal - dof2; d < dofTotal; ++d) {
109: if (a2[j][i][d] != expected) {
110: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
111: }
112: }
113: }
114: }
115: }
116: PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal2, &a2));
117: }
119: PetscCall(VecDestroy(&vec));
120: PetscCall(VecDestroy(&vecLocal1));
121: PetscCall(VecDestroy(&vecLocal2));
122: PetscCall(DMDestroy(&dm));
123: PetscCall(PetscFinalize());
124: return 0;
125: }
127: /*TEST
129: test:
130: suffix: 1
131: nsize: 4
132: args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_stencil_width 1 -stag_dof_2 2
133: output_file: output/empty.out
135: test:
136: suffix: 2
137: nsize: 4
138: args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_dof_1 2 -stag_grid_y 5
139: output_file: output/empty.out
141: test:
142: suffix: 3
143: nsize: 6
144: args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 2 -stag_grid_x 6
145: output_file: output/empty.out
147: test:
148: suffix: 4
149: nsize: 6
150: args: -stag_ranks_x 3 -stag_ranks_y 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_stencil_width 1
151: output_file: output/empty.out
153: test:
154: suffix: 5
155: nsize: 6
156: args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 0 -stag_dof_1 0 -stag_dof_2 0 -stag_grid_x 4 -stag_boundary_type_x ghosted -stag_stencil_width 1
157: output_file: output/empty.out
159: test:
160: suffix: 6
161: nsize: 9
162: args: -stag_dof_0 2 -stag_dof_1 2 -stag_dof_2 1 -stag_dof_2 1 -stag_boundary_type_y ghosted -stag_grid_x 9 -stag_grid_y 13 -stag_ranks_x 3 -stag_ranks_y 3 -stag_stencil_width 1
163: output_file: output/empty.out
165: test:
166: suffix: 7
167: nsize: 1
168: args: -stag_dof_0 2 -stag_dof_1 2 -stag_dof_2 1 -stag_dof_2 1 stag_grid_x 9 -stag_grid_y 13 -stag_stencil_width 1
169: output_file: output/empty.out
171: TEST*/