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: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
12: {
13: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
34: {
35: PetscFunctionBeginUser;
36: r[0] = 0.;
37: r[1] = 0.;
38: r[2] = 0.;
39: r[3] = 0.0;
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: int main(int argc, char **args)
44: {
45: Mat C;
46: PetscMPIInt rank, size;
47: PetscInt i, m = 5, N, start, end, M, its;
48: PetscScalar val, Ke[16], r[4];
49: PetscReal x, y, h, norm;
50: PetscInt idx[4], count, *rows;
51: Vec u, ustar, b;
52: KSP ksp;
53: PetscBool viewkspest = PETSC_FALSE;
55: PetscFunctionBeginUser;
56: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
57: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
58: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_est_view", &viewkspest, NULL));
59: N = (m + 1) * (m + 1); /* dimension of matrix */
60: M = m * m; /* number of elements */
61: h = 1.0 / m; /* mesh width */
62: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
63: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
65: /* Create stiffness matrix */
66: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
67: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
68: PetscCall(MatSetFromOptions(C));
69: PetscCall(MatSetUp(C));
70: start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
71: end = start + M / size + ((M % size) > rank);
73: /* Assemble matrix */
74: PetscCall(FormElementStiffness(h * h, Ke)); /* element stiffness for Laplacian */
75: for (i = start; i < end; i++) {
76: /* node numbers for the four corners of element */
77: idx[0] = (m + 1) * (i / m) + (i % m);
78: idx[1] = idx[0] + 1;
79: idx[2] = idx[1] + m + 1;
80: idx[3] = idx[2] - 1;
81: PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES));
82: }
83: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
84: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
86: /* Create right-hand-side and solution vectors */
87: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
88: PetscCall(VecSetSizes(u, PETSC_DECIDE, N));
89: PetscCall(VecSetFromOptions(u));
90: PetscCall(PetscObjectSetName((PetscObject)u, "Approx. Solution"));
91: PetscCall(VecDuplicate(u, &b));
92: PetscCall(PetscObjectSetName((PetscObject)b, "Right hand side"));
93: PetscCall(VecDuplicate(b, &ustar));
94: PetscCall(VecSet(u, 0.0));
95: PetscCall(VecSet(b, 0.0));
97: /* Assemble right-hand-side vector */
98: for (i = start; i < end; i++) {
99: /* location of lower left corner of element */
100: x = h * (i % m);
101: y = h * (i / m);
102: /* node numbers for the four corners of element */
103: idx[0] = (m + 1) * (i / m) + (i % m);
104: idx[1] = idx[0] + 1;
105: idx[2] = idx[1] + m + 1;
106: idx[3] = idx[2] - 1;
107: PetscCall(FormElementRhs(x, y, h * h, r));
108: PetscCall(VecSetValues(b, 4, idx, r, ADD_VALUES));
109: }
110: PetscCall(VecAssemblyBegin(b));
111: PetscCall(VecAssemblyEnd(b));
113: /* Modify matrix and right-hand-side for Dirichlet boundary conditions */
114: PetscCall(PetscMalloc1(4 * m, &rows));
115: for (i = 0; i < m + 1; i++) {
116: rows[i] = i; /* bottom */
117: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
118: }
119: count = m + 1; /* left side */
120: for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
122: count = 2 * m; /* left side */
123: for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
124: for (i = 0; i < 4 * m; i++) {
125: val = h * (rows[i] / (m + 1));
126: PetscCall(VecSetValues(u, 1, &rows[i], &val, INSERT_VALUES));
127: PetscCall(VecSetValues(b, 1, &rows[i], &val, INSERT_VALUES));
128: }
129: PetscCall(MatZeroRows(C, 4 * m, rows, 1.0, 0, 0));
131: PetscCall(PetscFree(rows));
132: PetscCall(VecAssemblyBegin(u));
133: PetscCall(VecAssemblyEnd(u));
134: PetscCall(VecAssemblyBegin(b));
135: PetscCall(VecAssemblyEnd(b));
137: {
138: Mat A;
139: PetscCall(MatConvert(C, MATSAME, MAT_INITIAL_MATRIX, &A));
140: PetscCall(MatDestroy(&C));
141: PetscCall(MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C));
142: PetscCall(MatDestroy(&A));
143: }
145: /* Solve linear system */
146: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
147: PetscCall(KSPSetOperators(ksp, C, C));
148: PetscCall(KSPSetFromOptions(ksp));
149: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
150: PetscCall(KSPSolve(ksp, b, u));
152: if (viewkspest) {
153: KSP kspest;
155: PetscCall(KSPChebyshevEstEigGetKSP(ksp, &kspest));
156: if (kspest) PetscCall(KSPView(kspest, PETSC_VIEWER_STDOUT_WORLD));
157: }
159: /* Check error */
160: PetscCall(VecGetOwnershipRange(ustar, &start, &end));
161: for (i = start; i < end; i++) {
162: val = h * (i / (m + 1));
163: PetscCall(VecSetValues(ustar, 1, &i, &val, INSERT_VALUES));
164: }
165: PetscCall(VecAssemblyBegin(ustar));
166: PetscCall(VecAssemblyEnd(ustar));
167: PetscCall(VecAXPY(u, -1.0, ustar));
168: PetscCall(VecNorm(u, NORM_2, &norm));
169: PetscCall(KSPGetIterationNumber(ksp, &its));
170: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its));
172: /* Free work space */
173: PetscCall(KSPDestroy(&ksp));
174: PetscCall(VecDestroy(&ustar));
175: PetscCall(VecDestroy(&u));
176: PetscCall(VecDestroy(&b));
177: PetscCall(MatDestroy(&C));
178: PetscCall(PetscFinalize());
179: return 0;
180: }
182: /*TEST
184: test:
185: args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
187: test:
188: suffix: 2
189: nsize: 2
190: args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
192: test:
193: suffix: 2_kokkos
194: nsize: 2
195: args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always -mat_type aijkokkos -vec_type kokkos
196: output_file: output/ex3_2.out
197: requires: kokkos_kernels
199: test:
200: suffix: nocheby
201: args: -ksp_est_view
203: test:
204: suffix: chebynoest
205: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_eigenvalues 0.1,1.0
207: test:
208: suffix: chebyest
209: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_esteig
210: filter: sed -e "s/Iterations 19/Iterations 20/g"
212: TEST*/