1 /*
2 Laplacian in 3D. Use for testing BAIJ matrix.
3 Modeled by the partial differential equation
4
5 - Laplacian u = 1,0 < x,y,z < 1,
6
7 with boundary conditions
8 u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
9 */
10
11 static char help[] = "Solves 3D Laplacian using wirebasket based multigrid.\n\n";
12
13 #include <petscdm.h>
14 #include <petscdmda.h>
15 #include <petscksp.h>
16
17 extern PetscErrorCode ComputeMatrix(DM, Mat);
18 extern PetscErrorCode ComputeRHS(DM, Vec);
19
main(int argc,char ** argv)20 int main(int argc, char **argv)
21 {
22 KSP ksp;
23 PC pc;
24 Vec x, b;
25 DM da;
26 Mat A, Atrans;
27 PetscInt dof = 1, M = 8;
28 PetscBool flg, trans = PETSC_FALSE, sbaij = PETSC_FALSE;
29
30 PetscFunctionBeginUser;
31 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
32 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
33 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
34 PetscCall(PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL));
35 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sbaij", &sbaij, NULL));
36
37 PetscCall(DMDACreate(PETSC_COMM_WORLD, &da));
38 PetscCall(DMSetDimension(da, 3));
39 PetscCall(DMDASetBoundaryType(da, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE));
40 PetscCall(DMDASetStencilType(da, DMDA_STENCIL_STAR));
41 PetscCall(DMDASetSizes(da, M, M, M));
42 PetscCall(DMDASetNumProcs(da, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
43 PetscCall(DMDASetDof(da, dof));
44 PetscCall(DMDASetStencilWidth(da, 1));
45 PetscCall(DMDASetOwnershipRanges(da, NULL, NULL, NULL));
46 PetscCall(DMSetFromOptions(da));
47 PetscCall(DMSetUp(da));
48
49 PetscCall(DMCreateGlobalVector(da, &x));
50 PetscCall(DMCreateGlobalVector(da, &b));
51 PetscCall(ComputeRHS(da, b));
52 PetscCall(DMSetMatType(da, MATBAIJ));
53 if (sbaij) PetscCall(DMSetMatType(da, MATSBAIJ));
54 PetscCall(DMCreateMatrix(da, &A));
55 PetscCall(ComputeMatrix(da, A));
56
57 /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
58 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
59 PetscCall(MatAXPY(A, 1.0, Atrans, DIFFERENT_NONZERO_PATTERN));
60 PetscCall(MatScale(A, 0.5));
61 PetscCall(MatDestroy(&Atrans));
62
63 /* Test sbaij matrix */
64 flg = PETSC_FALSE;
65 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_sbaij1", &flg, NULL));
66 if (flg) {
67 Mat sA;
68 PetscBool issymm;
69 PetscCall(MatIsTranspose(A, A, 0.0, &issymm));
70 if (issymm) PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
71 else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: A is non-symmetric\n"));
72 PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
73 PetscCall(MatDestroy(&A));
74 A = sA;
75 }
76
77 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
78 PetscCall(KSPSetFromOptions(ksp));
79 PetscCall(KSPSetOperators(ksp, A, A));
80 PetscCall(KSPGetPC(ksp, &pc));
81 PetscCall(PCSetDM(pc, da));
82
83 if (trans) {
84 PetscCall(KSPSolveTranspose(ksp, b, x));
85 } else {
86 PetscCall(KSPSolve(ksp, b, x));
87 }
88
89 /* check final residual */
90 flg = PETSC_FALSE;
91 PetscCall(PetscOptionsGetBool(NULL, NULL, "-check_final_residual", &flg, NULL));
92 if (flg) {
93 Vec b1;
94 PetscReal norm;
95 PetscCall(KSPGetSolution(ksp, &x));
96 PetscCall(VecDuplicate(b, &b1));
97 PetscCall(MatMult(A, x, b1));
98 PetscCall(VecAXPY(b1, -1.0, b));
99 PetscCall(VecNorm(b1, NORM_2, &norm));
100 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final residual %g\n", (double)norm));
101 PetscCall(VecDestroy(&b1));
102 }
103
104 PetscCall(KSPDestroy(&ksp));
105 PetscCall(VecDestroy(&x));
106 PetscCall(VecDestroy(&b));
107 PetscCall(MatDestroy(&A));
108 PetscCall(DMDestroy(&da));
109 PetscCall(PetscFinalize());
110 return 0;
111 }
112
ComputeRHS(DM da,Vec b)113 PetscErrorCode ComputeRHS(DM da, Vec b)
114 {
115 PetscInt mx, my, mz;
116 PetscScalar h;
117
118 PetscFunctionBeginUser;
119 PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
120 h = 1.0 / ((mx - 1) * (my - 1) * (mz - 1));
121 PetscCall(VecSet(b, h));
122 PetscFunctionReturn(PETSC_SUCCESS);
123 }
124
ComputeMatrix(DM da,Mat B)125 PetscErrorCode ComputeMatrix(DM da, Mat B)
126 {
127 PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs, dof, k1, k2, k3;
128 PetscScalar *v, *v_neighbor, Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy;
129 MatStencil row, col;
130
131 PetscFunctionBeginUser;
132 PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
133 /* For simplicity, this example only works on mx=my=mz */
134 PetscCheck(mx == my && mx == mz, PETSC_COMM_SELF, PETSC_ERR_SUP, "This example only works with mx %" PetscInt_FMT " = my %" PetscInt_FMT " = mz %" PetscInt_FMT, mx, my, mz);
135
136 Hx = 1.0 / (PetscReal)(mx - 1);
137 Hy = 1.0 / (PetscReal)(my - 1);
138 Hz = 1.0 / (PetscReal)(mz - 1);
139 HxHydHz = Hx * Hy / Hz;
140 HxHzdHy = Hx * Hz / Hy;
141 HyHzdHx = Hy * Hz / Hx;
142
143 PetscCall(PetscMalloc1(2 * dof * dof + 1, &v));
144 v_neighbor = v + dof * dof;
145 PetscCall(PetscArrayzero(v, 2 * dof * dof + 1));
146 k3 = 0;
147 for (k1 = 0; k1 < dof; k1++) {
148 for (k2 = 0; k2 < dof; k2++) {
149 if (k1 == k2) {
150 v[k3] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx);
151 v_neighbor[k3] = -HxHydHz;
152 } else {
153 v[k3] = k1 / (dof * dof);
154 v_neighbor[k3] = k2 / (dof * dof);
155 }
156 k3++;
157 }
158 }
159 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
160
161 for (k = zs; k < zs + zm; k++) {
162 for (j = ys; j < ys + ym; j++) {
163 for (i = xs; i < xs + xm; i++) {
164 row.i = i;
165 row.j = j;
166 row.k = k;
167 if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) { /* boundary points */
168 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &row, v, INSERT_VALUES));
169 } else { /* interior points */
170 /* center */
171 col.i = i;
172 col.j = j;
173 col.k = k;
174 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v, INSERT_VALUES));
175
176 /* x neighbors */
177 col.i = i - 1;
178 col.j = j;
179 col.k = k;
180 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
181 col.i = i + 1;
182 col.j = j;
183 col.k = k;
184 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
185
186 /* y neighbors */
187 col.i = i;
188 col.j = j - 1;
189 col.k = k;
190 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
191 col.i = i;
192 col.j = j + 1;
193 col.k = k;
194 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
195
196 /* z neighbors */
197 col.i = i;
198 col.j = j;
199 col.k = k - 1;
200 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
201 col.i = i;
202 col.j = j;
203 col.k = k + 1;
204 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
205 }
206 }
207 }
208 }
209 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
210 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
211 PetscCall(PetscFree(v));
212 PetscFunctionReturn(PETSC_SUCCESS);
213 }
214
215 /*TEST
216
217 test:
218 args: -ksp_monitor_short -sbaij -ksp_monitor_short -pc_type cholesky -ksp_view
219
220 TEST*/
221