1 static char help[] = "Test DMStag 3d star stencil\n\n";
2 #include <petscdm.h>
3 #include <petscdmstag.h>
4
main(int argc,char ** argv)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,sum;
11 PetscInt startx,starty,startz,nx,ny,nz,i,j,k,d,is,js,ks,dof0,dof1,dof2,dof3,dofTotal,stencilWidth,ngx,ngy,ngz;
12 DMBoundaryType boundaryTypex,boundaryTypey,boundaryTypez;
13 PetscMPIInt rank;
14
15 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
16 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
17 dof0 = 1;
18 dof1 = 1;
19 dof2 = 1;
20 dof3 = 1;
21 stencilWidth = 2;
22 ierr = 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_STAR,stencilWidth,NULL,NULL,NULL,&dm);CHKERRQ(ierr);
23 ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
24 ierr = DMSetUp(dm);CHKERRQ(ierr);
25 ierr = DMStagGetDOF(dm,&dof0,&dof1,&dof2,&dof3);CHKERRQ(ierr);
26 dofTotal = dof0 + 3*dof1 + 3*dof2 + dof3;
27 ierr = DMStagGetStencilWidth(dm,&stencilWidth);CHKERRQ(ierr);
28
29 ierr = DMCreateLocalVector(dm,&vecLocal1);CHKERRQ(ierr);
30 ierr = VecDuplicate(vecLocal1,&vecLocal2);CHKERRQ(ierr);
31
32 ierr = DMCreateGlobalVector(dm,&vec);CHKERRQ(ierr);
33 ierr = VecSet(vec,1.0);CHKERRQ(ierr);
34 ierr = VecSet(vecLocal1,0.0);CHKERRQ(ierr);
35 ierr = DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal1);CHKERRQ(ierr);
36 ierr = DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal1);CHKERRQ(ierr);
37
38 ierr = DMStagGetCorners(dm,&startx,&starty,&startz,&nx,&ny,&nz,NULL,NULL,NULL);CHKERRQ(ierr);
39 ierr = DMStagVecGetArrayRead(dm,vecLocal1,&a1);CHKERRQ(ierr);
40 ierr = DMStagVecGetArray(dm,vecLocal2,&a2);CHKERRQ(ierr);
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 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,(double)PetscRealPart(a1[k][j][i][d]),1.0);CHKERRQ(ierr);
47 }
48 a2[k][j][i][d] = 0.0;
49 for (ks = -stencilWidth; ks <= stencilWidth; ++ks) {
50 a2[k][j][i][d] += a1[k+ks][j][i][d];
51 }
52 for (js = -stencilWidth; js <= stencilWidth; ++js) {
53 a2[k][j][i][d] += a1[k][j+js][i][d];
54 }
55 for (is = -stencilWidth; is <= stencilWidth; ++is) {
56 a2[k][j][i][d] += a1[k][j][i+is][d];
57 }
58 a2[k][j][i][d] -= 2.0*a1[k][j][i][d];
59 }
60 }
61 }
62 }
63 ierr = DMStagVecRestoreArrayRead(dm,vecLocal1,&a1);CHKERRQ(ierr);
64 ierr = DMStagVecRestoreArray(dm,vecLocal2,&a2);CHKERRQ(ierr);
65
66 DMLocalToGlobalBegin(dm,vecLocal2,INSERT_VALUES,vec);CHKERRQ(ierr);
67 DMLocalToGlobalEnd(dm,vecLocal2,INSERT_VALUES,vec);CHKERRQ(ierr);
68
69 /* For the all-periodic case, some additional checks */
70 ierr = DMStagGetBoundaryTypes(dm,&boundaryTypex,&boundaryTypey,&boundaryTypez);CHKERRQ(ierr);
71 if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
72
73 ierr = DMStagGetGhostCorners(dm,NULL,NULL,NULL,&ngx,&ngy,&ngz);CHKERRQ(ierr);
74 expected = (ngx*ngy*ngz - 8*stencilWidth*stencilWidth*stencilWidth - 4*stencilWidth*stencilWidth*(nx + ny + nz))*dofTotal;
75 ierr = VecSum(vecLocal1,&sum);CHKERRQ(ierr);
76 if (sum != expected) {
77 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected sum of local entries %g (expected %g)\n",rank,(double)PetscRealPart(sum),(double)PetscRealPart(expected));CHKERRQ(ierr);
78 }
79
80 ierr = VecGetArray(vec,&a);CHKERRQ(ierr);
81 expected = 1 + 6*stencilWidth;
82 for (i=0; i<nz*ny*nx*dofTotal; ++i) {
83 if (a[i] != expected) {
84 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,(double)PetscRealPart(a[i]),(double)PetscRealPart(expected));CHKERRQ(ierr);
85 }
86 }
87 ierr = VecRestoreArray(vec,&a);CHKERRQ(ierr);
88 }
89
90 ierr = VecDestroy(&vec);CHKERRQ(ierr);
91 ierr = VecDestroy(&vecLocal1);CHKERRQ(ierr);
92 ierr = VecDestroy(&vecLocal2);CHKERRQ(ierr);
93 ierr = DMDestroy(&dm);CHKERRQ(ierr);
94 ierr = PetscFinalize();
95 return ierr;
96 }
97
98 /*TEST
99
100 test:
101 suffix: 1
102 nsize: 8
103 args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_stencil_width 1
104
105 test:
106 suffix: 2
107 nsize: 12
108 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 2 -stag_grid_x 6
109
110 test:
111 suffix: 3
112 nsize: 8
113 args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_dof_3 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6 -stag_grid_z 7
114
115 test:
116 suffix: 4
117 nsize: 8
118 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_x ghosted
119
120 test:
121 suffix: 5
122 nsize: 8
123 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_y ghosted
124
125 test:
126 suffix: 6
127 nsize: 8
128 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_z ghosted -stag_dof_0 2 -stag_dof_1 3 -stag_dof_2 2 -stag_dof_3 2
129
130 test:
131 suffix: 7
132 nsize: 8
133 args: -stag_stencil_width 1 -stag_grid_x 3 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_x ghosted -stag_boundary_type_y ghosted
134
135 test:
136 suffix: 8
137 nsize: 8
138 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 5 -stag_grid_z 2 -stag_boundary_type_x ghosted -stag_boundary_type_z ghosted
139
140 test:
141 suffix: 9
142 nsize: 8
143 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 3 -stag_boundary_type_y ghosted -stag_boundary_type_z ghosted
144
145 test:
146 suffix: 10
147 nsize: 5
148 args: -stag_stencil_width 1 -stag_grid_y 2 -stag_grid_z 3 -stag_grid_x 17 -stag_boundary_type_y ghosted -stag_boundary_type_z ghosted -stag_ranks_x 5
149 TEST*/
150