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: /*
8: Include "petscksp.h" so that we can use KSP solvers. Note that this file
9: automatically includes:
10: petscsys.h - base PETSc routines petscvec.h - vectors
11: petscmat.h - matrices
12: petscis.h - index sets petscksp.h - Krylov subspace methods
13: petscviewer.h - viewers petscpc.h - preconditioners
14: */
15: #include <petscksp.h>
17: /* Declare user-defined routines */
18: extern PetscErrorCode FormElementStiffness(PetscReal, PetscScalar *);
19: extern PetscErrorCode FormElementRhs(PetscScalar, PetscScalar, PetscReal, PetscScalar *);
21: int main(int argc, char **args)
22: {
23: Vec u, b, ustar; /* approx solution, RHS, exact solution */
24: Mat A; /* linear system matrix */
25: KSP ksp; /* Krylov subspace method context */
26: PetscInt N; /* dimension of system (global) */
27: PetscInt M; /* number of elements (global) */
28: PetscMPIInt rank; /* processor rank */
29: PetscMPIInt size; /* size of communicator */
30: PetscScalar Ke[16]; /* element matrix */
31: PetscScalar r[4]; /* element vector */
32: PetscReal h; /* mesh width */
33: PetscReal norm; /* norm of solution error */
34: PetscScalar x, y;
35: PetscInt idx[4], count, *rows, i, m = 5, start, end, its;
38: PetscInitialize(&argc, &args, (char *)0, help);
39: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
40: N = (m + 1) * (m + 1);
41: M = m * m;
42: h = 1.0 / m;
43: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
44: MPI_Comm_size(PETSC_COMM_WORLD, &size);
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Compute the matrix and right-hand-side vector that define
48: the linear system, Au = b.
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: /*
52: Create stiffness matrix
53: */
54: MatCreate(PETSC_COMM_WORLD, &A);
55: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N);
56: MatSetFromOptions(A);
57: MatSeqAIJSetPreallocation(A, 9, NULL);
58: MatMPIAIJSetPreallocation(A, 9, NULL, 8, NULL);
59: start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
60: end = start + M / size + ((M % size) > rank);
62: /*
63: Assemble matrix
64: */
65: FormElementStiffness(h * h, Ke);
66: for (i = start; i < end; i++) {
67: /* node numbers for the four corners of element */
68: idx[0] = (m + 1) * (i / m) + (i % m);
69: idx[1] = idx[0] + 1;
70: idx[2] = idx[1] + m + 1;
71: idx[3] = idx[2] - 1;
72: MatSetValues(A, 4, idx, 4, idx, Ke, ADD_VALUES);
73: }
74: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
77: /*
78: Create right-hand-side and solution vectors
79: */
80: VecCreate(PETSC_COMM_WORLD, &u);
81: VecSetSizes(u, PETSC_DECIDE, N);
82: VecSetFromOptions(u);
83: PetscObjectSetName((PetscObject)u, "Approx. Solution");
84: VecDuplicate(u, &b);
85: PetscObjectSetName((PetscObject)b, "Right hand side");
86: VecDuplicate(b, &ustar);
87: VecSet(u, 0.0);
88: VecSet(b, 0.0);
90: /*
91: Assemble right-hand-side vector
92: */
93: for (i = start; i < end; i++) {
94: /* location of lower left corner of element */
95: x = h * (i % m);
96: y = h * (i / m);
97: /* node numbers for the four corners of element */
98: idx[0] = (m + 1) * (i / m) + (i % m);
99: idx[1] = idx[0] + 1;
100: idx[2] = idx[1] + m + 1;
101: idx[3] = idx[2] - 1;
102: FormElementRhs(x, y, h * h, r);
103: VecSetValues(b, 4, idx, r, ADD_VALUES);
104: }
105: VecAssemblyBegin(b);
106: VecAssemblyEnd(b);
108: /*
109: Modify matrix and right-hand-side for Dirichlet boundary conditions
110: */
111: PetscMalloc1(4 * m, &rows);
112: for (i = 0; i < m + 1; i++) {
113: rows[i] = i; /* bottom */
114: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
115: }
116: count = m + 1; /* left side */
117: for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
118: count = 2 * m; /* left side */
119: for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
120: for (i = 0; i < 4 * m; i++) {
121: y = h * (rows[i] / (m + 1));
122: VecSetValues(u, 1, &rows[i], &y, INSERT_VALUES);
123: VecSetValues(b, 1, &rows[i], &y, INSERT_VALUES);
124: }
125: MatZeroRows(A, 4 * m, rows, 1.0, 0, 0);
126: PetscFree(rows);
128: VecAssemblyBegin(u);
129: VecAssemblyEnd(u);
130: VecAssemblyBegin(b);
131: VecAssemblyEnd(b);
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Create the linear solver and set various options
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: KSPCreate(PETSC_COMM_WORLD, &ksp);
138: KSPSetOperators(ksp, A, A);
139: KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
140: KSPSetFromOptions(ksp);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Solve the linear system
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: KSPSolve(ksp, b, u);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Check solution and clean up
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: /* Check error */
153: VecGetOwnershipRange(ustar, &start, &end);
154: for (i = start; i < end; i++) {
155: y = h * (i / (m + 1));
156: VecSetValues(ustar, 1, &i, &y, INSERT_VALUES);
157: }
158: VecAssemblyBegin(ustar);
159: VecAssemblyEnd(ustar);
160: VecAXPY(u, -1.0, ustar);
161: VecNorm(u, NORM_2, &norm);
162: KSPGetIterationNumber(ksp, &its);
163: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its);
165: /*
166: Free work space. All PETSc objects should be destroyed when they
167: are no longer needed.
168: */
169: KSPDestroy(&ksp);
170: VecDestroy(&u);
171: VecDestroy(&ustar);
172: VecDestroy(&b);
173: MatDestroy(&A);
175: /*
176: Always call PetscFinalize() before exiting a program. This routine
177: - finalizes the PETSc libraries as well as MPI
178: - provides summary and diagnostic information if certain runtime
179: options are chosen (e.g., -log_view).
180: */
181: PetscFinalize();
182: return 0;
183: }
185: /* --------------------------------------------------------------------- */
186: /* element stiffness for Laplacian */
187: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
188: {
190: Ke[0] = H / 6.0;
191: Ke[1] = -.125 * H;
192: Ke[2] = H / 12.0;
193: Ke[3] = -.125 * H;
194: Ke[4] = -.125 * H;
195: Ke[5] = H / 6.0;
196: Ke[6] = -.125 * H;
197: Ke[7] = H / 12.0;
198: Ke[8] = H / 12.0;
199: Ke[9] = -.125 * H;
200: Ke[10] = H / 6.0;
201: Ke[11] = -.125 * H;
202: Ke[12] = -.125 * H;
203: Ke[13] = H / 12.0;
204: Ke[14] = -.125 * H;
205: Ke[15] = H / 6.0;
206: return 0;
207: }
208: /* --------------------------------------------------------------------- */
209: PetscErrorCode FormElementRhs(PetscScalar x, PetscScalar y, PetscReal H, PetscScalar *r)
210: {
212: r[0] = 0.;
213: r[1] = 0.;
214: r[2] = 0.;
215: r[3] = 0.0;
216: return 0;
217: }
219: /*TEST
221: test:
222: suffix: 1
223: nsize: 2
224: args: -ksp_monitor_short
226: TEST*/