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));
86: PetscCall(VecSet(u, 0.0));
87: PetscCall(VecSet(b, 0.0));
89: /*
90: Assemble right-hand-side vector
91: */
92: for (i = start; i < end; i++) {
93: /* location of lower left corner of element */
94: x = h * (i % m);
95: y = h * (i / m);
96: /* node numbers for the four corners of element */
97: idx[0] = (m + 1) * (i / m) + (i % m);
98: idx[1] = idx[0] + 1;
99: idx[2] = idx[1] + m + 1;
100: idx[3] = idx[2] - 1;
101: PetscCall(FormElementRhs(x, y, h * h, r));
102: PetscCall(VecSetValues(b, 4, idx, r, ADD_VALUES));
103: }
104: PetscCall(VecAssemblyBegin(b));
105: PetscCall(VecAssemblyEnd(b));
107: /*
108: Modify matrix and right-hand side for Dirichlet boundary conditions
109: */
110: PetscCall(PetscMalloc1(4 * m, &rows));
111: for (i = 0; i < m + 1; i++) {
112: rows[i] = i; /* bottom */
113: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
114: }
115: count = m + 1; /* left side */
116: for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
117: count = 2 * m; /* left side */
118: for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
119: for (i = 0; i < 4 * m; i++) {
120: y = h * (rows[i] / (m + 1));
121: PetscCall(VecSetValues(u, 1, &rows[i], &y, INSERT_VALUES));
122: PetscCall(VecSetValues(b, 1, &rows[i], &y, INSERT_VALUES));
123: }
124: PetscCall(MatZeroRows(A, 4 * m, rows, 1.0, 0, 0));
125: PetscCall(PetscFree(rows));
127: PetscCall(VecAssemblyBegin(u));
128: PetscCall(VecAssemblyEnd(u));
129: PetscCall(VecAssemblyBegin(b));
130: PetscCall(VecAssemblyEnd(b));
132: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: Create the linear solver and set various options
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
137: PetscCall(KSPSetOperators(ksp, A, A));
138: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
139: PetscCall(KSPSetFromOptions(ksp));
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Solve the linear system
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: PetscCall(KSPSolve(ksp, b, u));
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Check solution and clean up
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: /* Check error */
152: PetscCall(VecGetOwnershipRange(ustar, &start, &end));
153: for (i = start; i < end; i++) {
154: y = h * (i / (m + 1));
155: PetscCall(VecSetValues(ustar, 1, &i, &y, INSERT_VALUES));
156: }
157: PetscCall(VecAssemblyBegin(ustar));
158: PetscCall(VecAssemblyEnd(ustar));
159: PetscCall(VecAXPY(u, -1.0, ustar));
160: PetscCall(VecNorm(u, NORM_2, &norm));
161: PetscCall(KSPGetIterationNumber(ksp, &its));
162: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its));
164: /*
165: Free work space. All PETSc objects should be destroyed when they
166: are no longer needed.
167: */
168: PetscCall(KSPDestroy(&ksp));
169: PetscCall(VecDestroy(&u));
170: PetscCall(VecDestroy(&ustar));
171: PetscCall(VecDestroy(&b));
172: PetscCall(MatDestroy(&A));
174: /*
175: Always call PetscFinalize() before exiting a program. This routine
176: - finalizes the PETSc libraries as well as MPI
177: - provides summary and diagnostic information if certain runtime
178: options are chosen (e.g., -log_view).
179: */
180: PetscCall(PetscFinalize());
181: return 0;
182: }
184: /* --------------------------------------------------------------------- */
185: /* element stiffness for Laplacian */
186: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
187: {
188: PetscFunctionBeginUser;
189: Ke[0] = H / 6.0;
190: Ke[1] = -.125 * H;
191: Ke[2] = H / 12.0;
192: Ke[3] = -.125 * H;
193: Ke[4] = -.125 * H;
194: Ke[5] = H / 6.0;
195: Ke[6] = -.125 * H;
196: Ke[7] = H / 12.0;
197: Ke[8] = H / 12.0;
198: Ke[9] = -.125 * H;
199: Ke[10] = H / 6.0;
200: Ke[11] = -.125 * H;
201: Ke[12] = -.125 * H;
202: Ke[13] = H / 12.0;
203: Ke[14] = -.125 * H;
204: Ke[15] = H / 6.0;
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
207: /* --------------------------------------------------------------------- */
208: PetscErrorCode FormElementRhs(PetscScalar x, PetscScalar y, PetscReal H, PetscScalar *r)
209: {
210: PetscFunctionBeginUser;
211: r[0] = 0.;
212: r[1] = 0.;
213: r[2] = 0.;
214: r[3] = 0.0;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*TEST
220: test:
221: suffix: 1
222: nsize: 2
223: args: -ksp_monitor_short
225: TEST*/