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