xref: /petsc/src/dm/tests/noflux_check.c (revision 18f812f49c83380c8efeb1af74ed9e41c2862ea8) 
1*18f812f4SBarry Smith static char help[] = "Check to see of DM_BOUNDARY_MIRROR works in 3D for DMDA with star stencil\n";
2*18f812f4SBarry Smith 
3*18f812f4SBarry Smith #include "petscdmda.h"
4*18f812f4SBarry Smith 
5*18f812f4SBarry Smith /* Contributed by Gourav Kumbhojkar */
6*18f812f4SBarry Smith 
7*18f812f4SBarry Smith PetscErrorCode globalKMat_3d(Mat K, DMDALocalInfo info)
8*18f812f4SBarry Smith {
9*18f812f4SBarry Smith   MatStencil  row, col[7];
10*18f812f4SBarry Smith   PetscScalar vals[7];
11*18f812f4SBarry Smith   PetscInt    ncols;
12*18f812f4SBarry Smith 
13*18f812f4SBarry Smith   PetscFunctionBeginUser;
14*18f812f4SBarry Smith   for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
15*18f812f4SBarry Smith     for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
16*18f812f4SBarry Smith       for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
17*18f812f4SBarry Smith         ncols = 0;
18*18f812f4SBarry Smith         row.i = i;
19*18f812f4SBarry Smith         row.j = j;
20*18f812f4SBarry Smith         row.k = k;
21*18f812f4SBarry Smith 
22*18f812f4SBarry Smith         col[0].i      = i;
23*18f812f4SBarry Smith         col[0].j      = j;
24*18f812f4SBarry Smith         col[0].k      = k;
25*18f812f4SBarry Smith         vals[ncols++] = -6.; //ncols=1
26*18f812f4SBarry Smith 
27*18f812f4SBarry Smith         col[ncols].i  = i - 1;
28*18f812f4SBarry Smith         col[ncols].j  = j;
29*18f812f4SBarry Smith         col[ncols].k  = k;
30*18f812f4SBarry Smith         vals[ncols++] = 1.; //ncols=2
31*18f812f4SBarry Smith 
32*18f812f4SBarry Smith         col[ncols].i  = i + 1;
33*18f812f4SBarry Smith         col[ncols].j  = j;
34*18f812f4SBarry Smith         col[ncols].k  = k;
35*18f812f4SBarry Smith         vals[ncols++] = 1.; //ncols=3
36*18f812f4SBarry Smith 
37*18f812f4SBarry Smith         col[ncols].i  = i;
38*18f812f4SBarry Smith         col[ncols].j  = j - 1;
39*18f812f4SBarry Smith         col[ncols].k  = k;
40*18f812f4SBarry Smith         vals[ncols++] = 1.; //ncols=4
41*18f812f4SBarry Smith 
42*18f812f4SBarry Smith         col[ncols].i  = i;
43*18f812f4SBarry Smith         col[ncols].j  = j + 1;
44*18f812f4SBarry Smith         col[ncols].k  = k;
45*18f812f4SBarry Smith         vals[ncols++] = 1.; //ncols=5
46*18f812f4SBarry Smith 
47*18f812f4SBarry Smith         col[ncols].i  = i;
48*18f812f4SBarry Smith         col[ncols].j  = j;
49*18f812f4SBarry Smith         col[ncols].k  = k + 1;
50*18f812f4SBarry Smith         vals[ncols++] = 1.; //ncols=6
51*18f812f4SBarry Smith 
52*18f812f4SBarry Smith         col[ncols].i  = i;
53*18f812f4SBarry Smith         col[ncols].j  = j;
54*18f812f4SBarry Smith         col[ncols].k  = k - 1;
55*18f812f4SBarry Smith         vals[ncols++] = 1.; //ncols=7
56*18f812f4SBarry Smith 
57*18f812f4SBarry Smith         PetscCall(MatSetValuesStencil(K, 1, &row, ncols, col, vals, ADD_VALUES));
58*18f812f4SBarry Smith       }
59*18f812f4SBarry Smith     }
60*18f812f4SBarry Smith   }
61*18f812f4SBarry Smith   PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
62*18f812f4SBarry Smith   PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));
63*18f812f4SBarry Smith   PetscFunctionReturn(PETSC_SUCCESS);
64*18f812f4SBarry Smith }
65*18f812f4SBarry Smith 
66*18f812f4SBarry Smith PetscErrorCode globalKMat_2d(Mat K, DMDALocalInfo info)
67*18f812f4SBarry Smith {
68*18f812f4SBarry Smith   MatStencil  row, col[5];
69*18f812f4SBarry Smith   PetscScalar vals[5];
70*18f812f4SBarry Smith   PetscInt    ncols;
71*18f812f4SBarry Smith 
72*18f812f4SBarry Smith   PetscFunctionBeginUser;
73*18f812f4SBarry Smith   for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
74*18f812f4SBarry Smith     for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
75*18f812f4SBarry Smith       ncols = 0;
76*18f812f4SBarry Smith       row.i = i;
77*18f812f4SBarry Smith       row.j = j;
78*18f812f4SBarry Smith 
79*18f812f4SBarry Smith       col[0].i      = i;
80*18f812f4SBarry Smith       col[0].j      = j;
81*18f812f4SBarry Smith       vals[ncols++] = -4.; //ncols=1
82*18f812f4SBarry Smith 
83*18f812f4SBarry Smith       col[ncols].i  = i - 1;
84*18f812f4SBarry Smith       col[ncols].j  = j;
85*18f812f4SBarry Smith       vals[ncols++] = 1.; //ncols=2
86*18f812f4SBarry Smith 
87*18f812f4SBarry Smith       col[ncols].i  = i;
88*18f812f4SBarry Smith       col[ncols].j  = j - 1;
89*18f812f4SBarry Smith       vals[ncols++] = 1.; //ncols=3
90*18f812f4SBarry Smith 
91*18f812f4SBarry Smith       col[ncols].i  = i + 1;
92*18f812f4SBarry Smith       col[ncols].j  = j;
93*18f812f4SBarry Smith       vals[ncols++] = 1.; //ncols=4
94*18f812f4SBarry Smith 
95*18f812f4SBarry Smith       col[ncols].i  = i;
96*18f812f4SBarry Smith       col[ncols].j  = j + 1;
97*18f812f4SBarry Smith       vals[ncols++] = 1.; //ncols=5
98*18f812f4SBarry Smith 
99*18f812f4SBarry Smith       PetscCall(MatSetValuesStencil(K, 1, &row, ncols, col, vals, ADD_VALUES));
100*18f812f4SBarry Smith     }
101*18f812f4SBarry Smith   }
102*18f812f4SBarry Smith   PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
103*18f812f4SBarry Smith   PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));
104*18f812f4SBarry Smith   PetscFunctionReturn(PETSC_SUCCESS);
105*18f812f4SBarry Smith }
106*18f812f4SBarry Smith 
107*18f812f4SBarry Smith int main(int argc, char **argv)
108*18f812f4SBarry Smith {
109*18f812f4SBarry Smith   DM                     da3d, da2d;
110*18f812f4SBarry Smith   DMDALocalInfo          info3d, info2d;
111*18f812f4SBarry Smith   Mat                    K3d, K2d;
112*18f812f4SBarry Smith   PetscInt               ne, num_pts;
113*18f812f4SBarry Smith   ISLocalToGlobalMapping ltgm3d, ltgm2d;
114*18f812f4SBarry Smith   Vec                    row2d, row3d;
115*18f812f4SBarry Smith   PetscReal              norm2d, norm3d;
116*18f812f4SBarry Smith 
117*18f812f4SBarry Smith   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
118*18f812f4SBarry Smith   ne      = 8;
119*18f812f4SBarry Smith   num_pts = ne + 1;
120*18f812f4SBarry Smith 
121*18f812f4SBarry Smith   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DMDA_STENCIL_STAR, num_pts, num_pts + 1, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da2d));
122*18f812f4SBarry Smith   PetscCall(DMSetUp(da2d));
123*18f812f4SBarry Smith   PetscCall(DMDAGetLocalInfo(da2d, &info2d));
124*18f812f4SBarry Smith   PetscCall(DMCreateMatrix(da2d, &K2d));
125*18f812f4SBarry Smith   PetscCall(DMGetLocalToGlobalMapping(da2d, &ltgm2d));
126*18f812f4SBarry Smith   PetscCall(ISLocalToGlobalMappingView(ltgm2d, PETSC_VIEWER_STDOUT_WORLD));
127*18f812f4SBarry Smith   //PetscFinalize();
128*18f812f4SBarry Smith   PetscCall(globalKMat_2d(K2d, info2d));
129*18f812f4SBarry Smith   PetscCall(MatView(K2d, PETSC_VIEWER_STDOUT_WORLD));
130*18f812f4SBarry Smith   PetscCall(MatCreateVecs(K2d, &row2d, NULL));
131*18f812f4SBarry Smith 
132*18f812f4SBarry Smith   PetscCall(MatGetRowSum(K2d, row2d));
133*18f812f4SBarry Smith   PetscCall(VecNorm(row2d, NORM_2, &norm2d));
134*18f812f4SBarry Smith 
135*18f812f4SBarry Smith   PetscCheck(norm2d == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "2D atrix row sum should be zero");
136*18f812f4SBarry Smith   PetscCall(VecDestroy(&row2d));
137*18f812f4SBarry Smith   PetscCall(MatDestroy(&K2d));
138*18f812f4SBarry Smith   PetscCall(DMDestroy(&da2d));
139*18f812f4SBarry Smith 
140*18f812f4SBarry Smith   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DMDA_STENCIL_STAR, num_pts, num_pts + 1, num_pts + 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &da3d));
141*18f812f4SBarry Smith   PetscCall(DMSetUp(da3d));
142*18f812f4SBarry Smith   PetscCall(DMCreateMatrix(da3d, &K3d));
143*18f812f4SBarry Smith   PetscCall(DMDAGetLocalInfo(da3d, &info3d));
144*18f812f4SBarry Smith   PetscCall(DMGetLocalToGlobalMapping(da3d, &ltgm3d));
145*18f812f4SBarry Smith   PetscCall(ISLocalToGlobalMappingView(ltgm3d, PETSC_VIEWER_STDOUT_WORLD));
146*18f812f4SBarry Smith   PetscCall(globalKMat_3d(K3d, info3d));
147*18f812f4SBarry Smith   PetscCall(MatView(K3d, PETSC_VIEWER_STDOUT_WORLD));
148*18f812f4SBarry Smith   PetscCall(MatCreateVecs(K3d, &row3d, NULL));
149*18f812f4SBarry Smith   PetscCall(MatGetRowSum(K3d, row3d));
150*18f812f4SBarry Smith   PetscCall(VecNorm(row3d, NORM_2, &norm3d));
151*18f812f4SBarry Smith   PetscCheck(norm3d == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "3D atrix row sum should be zero");
152*18f812f4SBarry Smith   PetscCall(VecDestroy(&row3d));
153*18f812f4SBarry Smith 
154*18f812f4SBarry Smith   PetscCall(DMDestroy(&da3d));
155*18f812f4SBarry Smith   PetscCall(MatDestroy(&K3d));
156*18f812f4SBarry Smith   return PetscFinalize();
157*18f812f4SBarry Smith }
158*18f812f4SBarry Smith 
159*18f812f4SBarry Smith /*TEST
160*18f812f4SBarry Smith 
161*18f812f4SBarry Smith    test:
162*18f812f4SBarry Smith 
163*18f812f4SBarry Smith    test:
164*18f812f4SBarry Smith       suffix: 2
165*18f812f4SBarry Smith       nsize: 2
166*18f812f4SBarry Smith 
167*18f812f4SBarry Smith    test:
168*18f812f4SBarry Smith       suffix: 4
169*18f812f4SBarry Smith       nsize: 4
170*18f812f4SBarry Smith 
171*18f812f4SBarry Smith    test:
172*18f812f4SBarry Smith       suffix: 8
173*18f812f4SBarry Smith       nsize: 8
174*18f812f4SBarry Smith 
175*18f812f4SBarry Smith TEST*/
176