Actual source code: ex4.c
2: static char help[] = "Bilinear elements on the unit square for the Laplacian. Input arguments are:\n\
3: -m <size> : problem size\n\n";
5: #include <petscksp.h>
7: int FormElementStiffness(PetscReal H, PetscScalar *Ke)
8: {
9: Ke[0] = H / 6.0;
10: Ke[1] = -.125 * H;
11: Ke[2] = H / 12.0;
12: Ke[3] = -.125 * H;
13: Ke[4] = -.125 * H;
14: Ke[5] = H / 6.0;
15: Ke[6] = -.125 * H;
16: Ke[7] = H / 12.0;
17: Ke[8] = H / 12.0;
18: Ke[9] = -.125 * H;
19: Ke[10] = H / 6.0;
20: Ke[11] = -.125 * H;
21: Ke[12] = -.125 * H;
22: Ke[13] = H / 12.0;
23: Ke[14] = -.125 * H;
24: Ke[15] = H / 6.0;
25: return 0;
26: }
28: int FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
29: {
30: r[0] = 0.;
31: r[1] = 0.;
32: r[2] = 0.;
33: r[3] = 0.0;
34: return 0;
35: }
37: /* Note: this code is for testing purposes only. The assembly process is not scalable */
38: int main(int argc, char **args)
39: {
40: Mat C;
41: PetscInt i, m = 2, N, M, its, idx[4], count, *rows;
42: PetscScalar val, Ke[16], r[4];
43: PetscReal x, y, h, norm;
44: Vec u, ustar, b;
45: KSP ksp;
46: PetscMPIInt rank;
47: PetscBool usezerorows = PETSC_TRUE;
50: PetscInitialize(&argc, &args, (char *)0, help);
51: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
52: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
53: PetscOptionsGetBool(NULL, NULL, "-usezerorows", &usezerorows, NULL);
54: N = (m + 1) * (m + 1); /* dimension of matrix */
55: M = m * m; /* number of elements */
56: h = 1.0 / m; /* mesh width */
58: /* create stiffness matrix */
59: MatCreate(PETSC_COMM_WORLD, &C);
60: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N);
61: MatSetFromOptions(C);
62: MatMPIAIJSetPreallocation(C, 9, NULL, 9, NULL);
63: MatSeqAIJSetPreallocation(C, 9, NULL);
64: #if defined(PETSC_HAVE_HYPRE)
65: MatHYPRESetPreallocation(C, 9, NULL, 9, NULL);
66: #endif
68: /* forms the element stiffness for the Laplacian */
69: FormElementStiffness(h * h, Ke);
71: /* assemble the matrix: only process 0 adds the values, not scalable */
72: if (rank == 0) {
73: for (i = 0; i < M; i++) {
74: /* node numbers for the four corners of element */
75: idx[0] = (m + 1) * (i / m) + (i % m);
76: idx[1] = idx[0] + 1;
77: idx[2] = idx[1] + m + 1;
78: idx[3] = idx[2] - 1;
79: if (i == M - 1 && !usezerorows) { /* If MatZeroRows not supported -> make it non-singular */
80: for (PetscInt ii = 0; ii < 4; ii++) {
81: Ke[2 * 4 + ii] = 0.0;
82: Ke[ii * 4 + 2] = 0.0;
83: }
84: Ke[10] = 1.0;
85: }
86: MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
87: }
88: }
89: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
92: /* create right hand side and solution */
93: MatCreateVecs(C, &u, &b);
94: VecDuplicate(u, &ustar);
95: VecSet(u, 0.0);
96: VecSet(b, 0.0);
98: /* assemble the right hand side: only process 0 adds the values, not scalable */
99: if (rank == 0) {
100: for (i = 0; i < M; i++) {
101: /* location of lower left corner of element */
102: x = h * (i % m);
103: y = h * (i / m);
104: /* node numbers for the four corners of element */
105: idx[0] = (m + 1) * (i / m) + (i % m);
106: idx[1] = idx[0] + 1;
107: idx[2] = idx[1] + m + 1;
108: idx[3] = idx[2] - 1;
109: FormElementRhs(x, y, h * h, r);
110: VecSetValues(b, 4, idx, r, ADD_VALUES);
111: }
112: }
113: VecAssemblyBegin(b);
114: VecAssemblyEnd(b);
116: /* modify matrix and rhs for Dirichlet boundary conditions */
117: if (rank == 0) {
118: PetscMalloc1(4 * m + 1, &rows);
119: for (i = 0; i < m + 1; i++) {
120: rows[i] = i; /* bottom */
121: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
122: }
123: count = m + 1; /* left side */
124: for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
126: count = 2 * m; /* right side */
127: for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
129: for (i = 0; i < 4 * m; i++) {
130: val = h * (rows[i] / (m + 1));
131: VecSetValues(u, 1, &rows[i], &val, INSERT_VALUES);
132: VecSetValues(b, 1, &rows[i], &val, INSERT_VALUES);
133: }
134: if (usezerorows) MatZeroRows(C, 4 * m, rows, 1.0, 0, 0);
135: PetscFree(rows);
136: } else {
137: if (usezerorows) MatZeroRows(C, 0, NULL, 1.0, 0, 0);
138: }
139: VecAssemblyBegin(u);
140: VecAssemblyEnd(u);
141: VecAssemblyBegin(b);
142: VecAssemblyEnd(b);
144: if (!usezerorows) {
145: VecSet(ustar, 1.0);
146: MatMult(C, ustar, b);
147: }
149: /* solve linear system */
150: KSPCreate(PETSC_COMM_WORLD, &ksp);
151: KSPSetOperators(ksp, C, C);
152: KSPSetFromOptions(ksp);
153: KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
154: KSPSolve(ksp, b, u);
156: /* check error */
157: if (usezerorows) {
158: if (rank == 0) {
159: for (i = 0; i < N; i++) {
160: val = h * (i / (m + 1));
161: VecSetValues(ustar, 1, &i, &val, INSERT_VALUES);
162: }
163: }
164: VecAssemblyBegin(ustar);
165: VecAssemblyEnd(ustar);
166: }
168: VecAXPY(u, -1.0, ustar);
169: VecNorm(u, NORM_2, &norm);
170: KSPGetIterationNumber(ksp, &its);
171: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its);
173: KSPDestroy(&ksp);
174: VecDestroy(&ustar);
175: VecDestroy(&u);
176: VecDestroy(&b);
177: MatDestroy(&C);
178: PetscFinalize();
179: return 0;
180: }
182: /*TEST
184: test:
185: args: -ksp_monitor_short -m 5 -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
187: test:
188: suffix: 3
189: args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
191: test:
192: suffix: 5
193: args: -pc_type eisenstat -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
195: test:
196: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
197: suffix: hypre_device_none
198: output_file: output/ex4_hypre_none.out
199: nsize: {{1 2}}
200: args: -usezerorows 0 -mat_type hypre -pc_type none -m 5
202: test:
203: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
204: suffix: hypre_device_amg
205: nsize: {{1 2}}
206: args: -usezerorows 0 -mat_type hypre -pc_type hypre -m 25 -ksp_type cg -ksp_norm_type natural
208: TEST*/