Actual source code: ex7.c
petsc-3.11.4 2019-09-28
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: PetscErrorCode ierr;
8: DM dm;
9: Vec vec,vecLocal1,vecLocal2;
10: PetscScalar *a,****a1,****a2,expected;
11: PetscInt startx,starty,startz,nx,ny,nz,i,j,k,d,is,js,ks,dof0,dof1,dof2,dof3,dofTotal,stencilWidth,Nx,Ny,Nz;
12: DMBoundaryType boundaryTypex,boundaryTypey,boundaryTypez;
13: PetscMPIInt rank;
15: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
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: DMStagVecGetArrayDOFRead(dm,vecLocal1,&a1);
40: DMStagVecGetArrayDOF(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) {
46: PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,a1[k][j][i][d],1.0);
47: }
48: a2[k][j][i][d] = 0.0;
49: for (ks = -stencilWidth; ks <= stencilWidth; ++ks) {
50: for (js = -stencilWidth; js <= stencilWidth; ++js) {
51: for (is = -stencilWidth; is <= stencilWidth; ++is) {
52: a2[k][j][i][d] += a1[k+ks][j+js][i+is][d];
53: }
54: }
55: }
56: }
57: }
58: }
59: }
60: DMStagVecRestoreArrayDOFRead(dm,vecLocal1,&a1);
61: DMStagVecRestoreArrayDOF(dm,vecLocal2,&a2);
63: DMLocalToGlobalBegin(dm,vecLocal2,INSERT_VALUES,vec);
64: DMLocalToGlobalEnd(dm,vecLocal2,INSERT_VALUES,vec);
66: /* For the all-periodic case, all values are the same . Otherwise, just check the local version */
67: DMStagGetBoundaryTypes(dm,&boundaryTypex,&boundaryTypey,&boundaryTypez);
68: if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
69: VecGetArray(vec,&a);
70: expected = 1.0; for(d=0;d<3;++d) expected *= (2*stencilWidth+1);
71: for (i=0; i<nz*ny*nx*dofTotal; ++i) {
72: if (a[i] != expected) {
73: PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,a[i],expected);
74: }
75: }
76: VecRestoreArray(vec,&a);
77: } else {
78: DMStagVecGetArrayDOFRead(dm,vecLocal2,&a2);
79: DMStagGetGlobalSizes(dm,&Nx,&Ny,&Nz);
80: if (stencilWidth > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Check implemented assuming stencilWidth = 1");
81: for (k=startz; k<startz + nz; ++k) {
82: for (j=starty; j<starty + ny; ++j) {
83: for (i=startx; i<startx + nx; ++i) {
84: PetscInt dd,extra[3];
85: PetscBool bnd[3];
86: bnd[0] = (PetscBool)((i == 0 || i == Nx-1) && boundaryTypex != DM_BOUNDARY_PERIODIC);
87: bnd[1] = (PetscBool)((j == 0 || j == Ny-1) && boundaryTypey != DM_BOUNDARY_PERIODIC);
88: bnd[2] = (PetscBool)((k == 0 || k == Nz-1) && boundaryTypez != DM_BOUNDARY_PERIODIC);
89: extra[0] = i == Nx-1 && boundaryTypex != DM_BOUNDARY_PERIODIC ? 1 : 0;
90: extra[1] = j == Ny-1 && boundaryTypey != DM_BOUNDARY_PERIODIC ? 1 : 0;
91: extra[2] = k == Nz-1 && boundaryTypez != DM_BOUNDARY_PERIODIC ? 1 : 0;
92: { /* vertices */
93: PetscScalar expected = 1.0;
94: for (dd=0;dd<3;++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2*stencilWidth + 1);
95: for (d=0; d<dof0; ++d) {
96: if (a2[k][j][i][d] != expected) {
97: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
98: }
99: }
100: }
101: { /* back down edges */
102: PetscScalar expected = ((bnd[0] ? 1 : 2) * stencilWidth + 1);
103: for (dd=1;dd<3;++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2*stencilWidth + 1);
104: for (d=dof0; d<dof0+dof1; ++d) {
105: if (a2[k][j][i][d] != expected) {
106: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
107: }
108: }
109: }
110: { /* back left edges */
111: PetscScalar expected = ((bnd[1] ? 1 : 2) * stencilWidth + 1);
112: for (dd=0;dd<3;dd+=2) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2*stencilWidth + 1);
113: for (d=dof0+dof1; d<dof0+2*dof1; ++d) {
114: if (a2[k][j][i][d] != expected) {
115: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
116: }
117: }
118: }
119: { /* back faces */
120: PetscScalar expected = (bnd[2] ? stencilWidth + 1 + extra[2] : 2*stencilWidth + 1);
121: for (dd=0;dd<2;++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
122: for (d=dof0+2*dof1; d<dof0+2*dof1+dof2; ++d) {
123: if (a2[k][j][i][d] != expected) {
124: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
125: }
126: }
127: }
128: { /* down left edges */
129: PetscScalar expected = ((bnd[2] ? 1 : 2) * stencilWidth + 1);
130: for (dd=0;dd<2;++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2*stencilWidth + 1);
131: for (d=dof0+2*dof1+dof2; d<dof0+3*dof1+dof2; ++d) {
132: if (a2[k][j][i][d] != expected) {
133: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
134: }
135: }
136: }
137: { /* down faces */
138: PetscScalar expected = (bnd[1] ? stencilWidth + 1 + extra[1] : 2*stencilWidth + 1);
139: for (dd=0;dd<3;dd+=2) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
140: for (d=dof0+3*dof1+dof2; d<dof0+3*dof1+2*dof2; ++d) {
141: if (a2[k][j][i][d] != expected) {
142: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
143: }
144: }
145: }
146: { /* left faces */
147: PetscScalar expected = (bnd[0] ? stencilWidth + 1 + extra[0] : 2*stencilWidth + 1);
148: for (dd=1;dd<3;++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
149: for (d=dof0+3*dof1+2*dof2; d<dof0+3*dof1+3*dof2; ++d) {
150: if (a2[k][j][i][d] != expected) {
151: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
152: }
153: }
154: }
155: { /* elements */
156: PetscScalar expected = 1.0;
157: for (dd=0;dd<3;++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
158: for (d=dofTotal-dof3; d<dofTotal; ++d) {
159: if (a2[k][j][i][d] != expected) {
160: PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,k,d,a2[k][j][i][d],expected);
161: }
162: }
163: }
164: }
165: }
166: }
167: DMStagVecRestoreArrayDOFRead(dm,vecLocal2,&a2);
168: }
170: VecDestroy(&vec);
171: VecDestroy(&vecLocal1);
172: VecDestroy(&vecLocal2);
173: DMDestroy(&dm);
174: PetscFinalize();
175: return ierr;
176: }
178: /*TEST
180: test:
181: suffix: 1
182: nsize: 8
183: 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
185: test:
186: suffix: 2
187: nsize: 8
188: args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_2 2 -stag_grid_y 5
190: test:
191: suffix: 3
192: nsize: 12
193: args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 2 -stag_grid_x 6
195: test:
196: suffix: 4
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_y ghosted -stag_boundary_type_z ghosted -stag_stencil_width 1
200: test:
201: suffix: 5
202: nsize: 12
203: 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
205: test:
206: suffix: 6
207: nsize: 8
208: 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
209: TEST*/