Actual source code: ex4.c
1: static char help[] = "Bilinear elements on the unit square for the Laplacian. Input arguments are:\n\
2: -m <size> : problem size\n\n";
4: #include <petscksp.h>
6: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
7: {
8: Ke[0] = H / 6.0;
9: Ke[1] = -.125 * H;
10: Ke[2] = H / 12.0;
11: Ke[3] = -.125 * H;
12: Ke[4] = -.125 * H;
13: Ke[5] = H / 6.0;
14: Ke[6] = -.125 * H;
15: Ke[7] = H / 12.0;
16: Ke[8] = H / 12.0;
17: Ke[9] = -.125 * H;
18: Ke[10] = H / 6.0;
19: Ke[11] = -.125 * H;
20: Ke[12] = -.125 * H;
21: Ke[13] = H / 12.0;
22: Ke[14] = -.125 * H;
23: Ke[15] = H / 6.0;
24: return PETSC_SUCCESS;
25: }
27: PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
28: {
29: r[0] = 0.;
30: r[1] = 0.;
31: r[2] = 0.;
32: r[3] = 0.0;
33: return PETSC_SUCCESS;
34: }
36: /* Note: this code is for testing purposes only. The assembly process is not scalable */
37: int main(int argc, char **args)
38: {
39: Mat C;
40: PetscInt i, m = 2, N, M, its, idx[4], count, *rows;
41: PetscScalar val, Ke[16], r[4];
42: PetscReal x, y, h, norm;
43: Vec u, ustar, b;
44: KSP ksp;
45: PetscMPIInt rank;
46: PetscBool usezerorows = PETSC_TRUE;
48: PetscFunctionBeginUser;
49: PetscCall(PetscInitialize(&argc, &args, NULL, help));
50: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
51: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
52: PetscCall(PetscOptionsGetBool(NULL, NULL, "-usezerorows", &usezerorows, NULL));
53: N = (m + 1) * (m + 1); /* dimension of matrix */
54: M = m * m; /* number of elements */
55: h = 1.0 / m; /* mesh width */
57: /* create stiffness matrix */
58: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
59: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
60: PetscCall(MatSetFromOptions(C));
61: PetscCall(MatMPIAIJSetPreallocation(C, 9, NULL, 9, NULL));
62: PetscCall(MatSeqAIJSetPreallocation(C, 9, NULL));
63: #if defined(PETSC_HAVE_HYPRE)
64: PetscCall(MatHYPRESetPreallocation(C, 9, NULL, 9, NULL));
65: #endif
67: /* forms the element stiffness for the Laplacian */
68: PetscCall(FormElementStiffness(h * h, Ke));
70: /* assemble the matrix: only process 0 adds the values, not scalable */
71: if (rank == 0) {
72: for (i = 0; i < M; i++) {
73: /* node numbers for the four corners of element */
74: idx[0] = (m + 1) * (i / m) + (i % m);
75: idx[1] = idx[0] + 1;
76: idx[2] = idx[1] + m + 1;
77: idx[3] = idx[2] - 1;
78: if (i == M - 1 && !usezerorows) { /* If MatZeroRows not supported -> make it non-singular */
79: for (PetscInt ii = 0; ii < 4; ii++) {
80: Ke[2 * 4 + ii] = 0.0;
81: Ke[ii * 4 + 2] = 0.0;
82: }
83: Ke[10] = 1.0;
84: }
85: PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES));
86: }
87: }
88: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
89: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
91: /* create right-hand side and solution */
92: PetscCall(MatCreateVecs(C, &u, &b));
93: PetscCall(VecDuplicate(u, &ustar));
94: PetscCall(VecSet(u, 0.0));
95: PetscCall(VecSet(b, 0.0));
97: /* assemble the right-hand side: only MPI process with rank 0 adds the values, this is not scalable */
98: if (rank == 0) {
99: for (i = 0; i < M; i++) {
100: /* location of lower left corner of element */
101: x = h * (i % m);
102: y = h * (i / m);
103: /* node numbers for the four corners of element */
104: idx[0] = (m + 1) * (i / m) + (i % m);
105: idx[1] = idx[0] + 1;
106: idx[2] = idx[1] + m + 1;
107: idx[3] = idx[2] - 1;
108: PetscCall(FormElementRhs(x, y, h * h, r));
109: PetscCall(VecSetValues(b, 4, idx, r, ADD_VALUES));
110: }
111: }
112: PetscCall(VecAssemblyBegin(b));
113: PetscCall(VecAssemblyEnd(b));
115: /* modify matrix and rhs for Dirichlet boundary conditions */
116: if (rank == 0) {
117: PetscCall(PetscMalloc1(4 * m + 1, &rows));
118: for (i = 0; i < m + 1; i++) {
119: rows[i] = i; /* bottom */
120: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
121: }
122: count = m + 1; /* left side */
123: for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
125: count = 2 * m; /* right side */
126: for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
128: for (i = 0; i < 4 * m; i++) {
129: val = h * (rows[i] / (m + 1));
130: PetscCall(VecSetValues(u, 1, &rows[i], &val, INSERT_VALUES));
131: PetscCall(VecSetValues(b, 1, &rows[i], &val, INSERT_VALUES));
132: }
133: if (usezerorows) PetscCall(MatZeroRows(C, 4 * m, rows, 1.0, 0, 0));
134: PetscCall(PetscFree(rows));
135: } else {
136: if (usezerorows) PetscCall(MatZeroRows(C, 0, NULL, 1.0, 0, 0));
137: }
138: PetscCall(VecAssemblyBegin(u));
139: PetscCall(VecAssemblyEnd(u));
140: PetscCall(VecAssemblyBegin(b));
141: PetscCall(VecAssemblyEnd(b));
143: if (!usezerorows) {
144: PetscCall(VecSet(ustar, 1.0));
145: PetscCall(MatMult(C, ustar, b));
146: }
148: /* solve linear system */
149: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
150: PetscCall(KSPSetOperators(ksp, C, C));
151: PetscCall(KSPSetFromOptions(ksp));
152: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
153: PetscCall(KSPSolve(ksp, b, u));
155: /* check error */
156: if (usezerorows) {
157: if (rank == 0) {
158: for (i = 0; i < N; i++) {
159: val = h * (i / (m + 1));
160: PetscCall(VecSetValues(ustar, 1, &i, &val, INSERT_VALUES));
161: }
162: }
163: PetscCall(VecAssemblyBegin(ustar));
164: PetscCall(VecAssemblyEnd(ustar));
165: }
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: { // Test getting Jacobi diag
173: PC pc;
174: PetscBool is_pcjacobi;
176: PetscCall(KSPGetPC(ksp, &pc));
177: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_pcjacobi));
178: if (is_pcjacobi) {
179: Vec diag;
181: PetscCall(MatCreateVecs(C, &diag, NULL));
182: PetscCall(PCJacobiGetDiagonal(pc, diag, NULL));
183: PetscCall(VecNorm(diag, NORM_2, &norm));
184: PetscCheck(norm > 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Jacobi preconditioner should have norm greater than 0");
185: PetscCall(VecDestroy(&diag));
186: }
187: }
189: PetscCall(KSPDestroy(&ksp));
190: PetscCall(VecDestroy(&ustar));
191: PetscCall(VecDestroy(&u));
192: PetscCall(VecDestroy(&b));
193: PetscCall(MatDestroy(&C));
194: PetscCall(PetscFinalize());
195: return 0;
196: }
198: /*TEST
200: test:
201: args: -ksp_monitor_short -m 5 -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
203: test:
204: suffix: 3
205: args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
207: test:
208: suffix: 5
209: args: -pc_type eisenstat -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
211: test:
212: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
213: suffix: hypre_device_none
214: output_file: output/ex4_hypre_none.out
215: nsize: {{1 2}}
216: args: -usezerorows 0 -mat_type hypre -pc_type none -m 5
218: test:
219: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
220: suffix: hypre_device_amg
221: nsize: {{1 2}}
222: args: -usezerorows 0 -mat_type hypre -pc_type hypre -m 25 -ksp_type cg -ksp_norm_type natural
224: test:
225: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
226: suffix: dev
227: nsize: 2
228: args: -m 25 -mat_type hypre -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_coarsen_type PMIS
230: TEST*/