Actual source code: ex3.c
2: static char help[] = "Bilinear elements on the unit square for Laplacian. To test the parallel\n\
3: matrix assembly, the matrix is intentionally laid out across processors\n\
4: differently from the way it is assembled. Input arguments are:\n\
5: -m <size> : problem size\n\n";
7: /* Addendum: piggy-backing on this example to test KSPChebyshev methods */
9: #include <petscksp.h>
11: int FormElementStiffness(PetscReal H, PetscScalar *Ke)
12: {
14: Ke[0] = H / 6.0;
15: Ke[1] = -.125 * H;
16: Ke[2] = H / 12.0;
17: Ke[3] = -.125 * H;
18: Ke[4] = -.125 * H;
19: Ke[5] = H / 6.0;
20: Ke[6] = -.125 * H;
21: Ke[7] = H / 12.0;
22: Ke[8] = H / 12.0;
23: Ke[9] = -.125 * H;
24: Ke[10] = H / 6.0;
25: Ke[11] = -.125 * H;
26: Ke[12] = -.125 * H;
27: Ke[13] = H / 12.0;
28: Ke[14] = -.125 * H;
29: Ke[15] = H / 6.0;
30: return 0;
31: }
32: int FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
33: {
35: r[0] = 0.;
36: r[1] = 0.;
37: r[2] = 0.;
38: r[3] = 0.0;
39: return 0;
40: }
42: int main(int argc, char **args)
43: {
44: Mat C;
45: PetscMPIInt rank, size;
46: PetscInt i, m = 5, N, start, end, M, its;
47: PetscScalar val, Ke[16], r[4];
48: PetscReal x, y, h, norm;
49: PetscInt idx[4], count, *rows;
50: Vec u, ustar, b;
51: KSP ksp;
52: PetscBool viewkspest = PETSC_FALSE;
55: PetscInitialize(&argc, &args, (char *)0, help);
56: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
57: PetscOptionsGetBool(NULL, NULL, "-ksp_est_view", &viewkspest, NULL);
58: N = (m + 1) * (m + 1); /* dimension of matrix */
59: M = m * m; /* number of elements */
60: h = 1.0 / m; /* mesh width */
61: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
62: MPI_Comm_size(PETSC_COMM_WORLD, &size);
64: /* Create stiffness matrix */
65: MatCreate(PETSC_COMM_WORLD, &C);
66: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N);
67: MatSetFromOptions(C);
68: MatSetUp(C);
69: start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
70: end = start + M / size + ((M % size) > rank);
72: /* Assemble matrix */
73: FormElementStiffness(h * h, Ke); /* element stiffness for Laplacian */
74: for (i = start; i < end; i++) {
75: /* node numbers for the four corners of element */
76: idx[0] = (m + 1) * (i / m) + (i % m);
77: idx[1] = idx[0] + 1;
78: idx[2] = idx[1] + m + 1;
79: idx[3] = idx[2] - 1;
80: MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
81: }
82: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
85: /* Create right-hand-side and solution vectors */
86: VecCreate(PETSC_COMM_WORLD, &u);
87: VecSetSizes(u, PETSC_DECIDE, N);
88: VecSetFromOptions(u);
89: PetscObjectSetName((PetscObject)u, "Approx. Solution");
90: VecDuplicate(u, &b);
91: PetscObjectSetName((PetscObject)b, "Right hand side");
92: VecDuplicate(b, &ustar);
93: VecSet(u, 0.0);
94: VecSet(b, 0.0);
96: /* Assemble right-hand-side vector */
97: for (i = start; i < end; i++) {
98: /* location of lower left corner of element */
99: x = h * (i % m);
100: y = h * (i / m);
101: /* node numbers for the four corners of element */
102: idx[0] = (m + 1) * (i / m) + (i % m);
103: idx[1] = idx[0] + 1;
104: idx[2] = idx[1] + m + 1;
105: idx[3] = idx[2] - 1;
106: FormElementRhs(x, y, h * h, r);
107: VecSetValues(b, 4, idx, r, ADD_VALUES);
108: }
109: VecAssemblyBegin(b);
110: VecAssemblyEnd(b);
112: /* Modify matrix and right-hand-side for Dirichlet boundary conditions */
113: PetscMalloc1(4 * m, &rows);
114: for (i = 0; i < m + 1; i++) {
115: rows[i] = i; /* bottom */
116: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
117: }
118: count = m + 1; /* left side */
119: for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
121: count = 2 * m; /* left side */
122: for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
123: for (i = 0; i < 4 * m; i++) {
124: val = h * (rows[i] / (m + 1));
125: VecSetValues(u, 1, &rows[i], &val, INSERT_VALUES);
126: VecSetValues(b, 1, &rows[i], &val, INSERT_VALUES);
127: }
128: MatZeroRows(C, 4 * m, rows, 1.0, 0, 0);
130: PetscFree(rows);
131: VecAssemblyBegin(u);
132: VecAssemblyEnd(u);
133: VecAssemblyBegin(b);
134: VecAssemblyEnd(b);
136: {
137: Mat A;
138: MatConvert(C, MATSAME, MAT_INITIAL_MATRIX, &A);
139: MatDestroy(&C);
140: MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C);
141: MatDestroy(&A);
142: }
144: /* Solve linear system */
145: KSPCreate(PETSC_COMM_WORLD, &ksp);
146: KSPSetOperators(ksp, C, C);
147: KSPSetFromOptions(ksp);
148: KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
149: KSPSolve(ksp, b, u);
151: if (viewkspest) {
152: KSP kspest;
154: KSPChebyshevEstEigGetKSP(ksp, &kspest);
155: if (kspest) KSPView(kspest, PETSC_VIEWER_STDOUT_WORLD);
156: }
158: /* Check error */
159: VecGetOwnershipRange(ustar, &start, &end);
160: for (i = start; i < end; i++) {
161: val = h * (i / (m + 1));
162: VecSetValues(ustar, 1, &i, &val, INSERT_VALUES);
163: }
164: VecAssemblyBegin(ustar);
165: VecAssemblyEnd(ustar);
166: VecAXPY(u, -1.0, ustar);
167: VecNorm(u, NORM_2, &norm);
168: KSPGetIterationNumber(ksp, &its);
169: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its);
171: /* Free work space */
172: KSPDestroy(&ksp);
173: VecDestroy(&ustar);
174: VecDestroy(&u);
175: VecDestroy(&b);
176: MatDestroy(&C);
177: PetscFinalize();
178: return 0;
179: }
181: /*TEST
183: test:
184: args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
186: test:
187: suffix: 2
188: nsize: 2
189: args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
191: test:
192: suffix: 2_kokkos
193: nsize: 2
194: args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always -mat_type aijkokkos -vec_type kokkos
195: output_file: output/ex3_2.out
196: requires: kokkos_kernels
198: test:
199: suffix: nocheby
200: args: -ksp_est_view
202: test:
203: suffix: chebynoest
204: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_eigenvalues 0.1,1.0
206: test:
207: suffix: chebyest
208: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_esteig
209: filter: sed -e "s/Iterations 19/Iterations 20/g"
211: TEST*/