Actual source code: ex5.c
2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
3: This example tests PCVPBJacobiSetBlocks().\n\n";
5: /*
6: Include "petscsnes.h" so that we can use SNES solvers. Note that this
7: file automatically includes:
8: petscsys.h - base PETSc routines petscvec.h - vectors
9: petscmat.h - matrices
10: petscis.h - index sets petscksp.h - Krylov subspace methods
11: petscviewer.h - viewers petscpc.h - preconditioners
12: petscksp.h - linear solvers
13: */
15: #include <petscsnes.h>
17: /*
18: User-defined routines
19: */
20: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
21: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
22: extern PetscErrorCode FormInitialGuess(Vec);
24: int main(int argc, char **argv)
25: {
26: SNES snes; /* SNES context */
27: Vec x, r, F, U; /* vectors */
28: Mat J; /* Jacobian matrix */
29: PetscInt its, n = 5, i, maxit, maxf, lens[3] = {1, 2, 2};
30: PetscMPIInt size;
31: PetscScalar h, xp, v, none = -1.0;
32: PetscReal abstol, rtol, stol, norm;
33: KSP ksp;
34: PC pc;
37: PetscInitialize(&argc, &argv, (char *)0, help);
38: MPI_Comm_size(PETSC_COMM_WORLD, &size);
40: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
41: h = 1.0 / (n - 1);
43: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Create nonlinear solver context
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47: SNESCreate(PETSC_COMM_WORLD, &snes);
48: SNESGetKSP(snes, &ksp);
49: KSPGetPC(ksp, &pc);
50: PCSetType(pc, PCVPBJACOBI);
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create vector data structures; set function evaluation routine
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: /*
55: Note that we form 1 vector from scratch and then duplicate as needed.
56: */
57: VecCreate(PETSC_COMM_WORLD, &x);
58: VecSetSizes(x, PETSC_DECIDE, n);
59: VecSetFromOptions(x);
60: VecDuplicate(x, &r);
61: VecDuplicate(x, &F);
62: VecDuplicate(x, &U);
64: /*
65: Set function evaluation routine and vector
66: */
67: SNESSetFunction(snes, r, FormFunction, (void *)F);
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Create matrix data structure; set Jacobian evaluation routine
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: MatCreate(PETSC_COMM_WORLD, &J);
74: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n);
75: MatSetFromOptions(J);
76: MatSeqAIJSetPreallocation(J, 3, NULL);
77: MatSetVariableBlockSizes(J, 3, lens);
79: /*
80: Set Jacobian matrix data structure and default Jacobian evaluation
81: routine. User can override with:
82: -snes_fd : default finite differencing approximation of Jacobian
83: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
84: (unless user explicitly sets preconditioner)
85: -snes_mf_operator : form preconditioning matrix as set by the user,
86: but use matrix-free approx for Jacobian-vector
87: products within Newton-Krylov method
88: */
90: SNESSetJacobian(snes, J, J, FormJacobian, NULL);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Customize nonlinear solver; set runtime options
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: /*
97: Set names for some vectors to facilitate monitoring (optional)
98: */
99: PetscObjectSetName((PetscObject)x, "Approximate Solution");
100: PetscObjectSetName((PetscObject)U, "Exact Solution");
102: /*
103: Set SNES/KSP/KSP/PC runtime options, e.g.,
104: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
105: */
106: SNESSetFromOptions(snes);
108: /*
109: Print parameters used for convergence testing (optional) ... just
110: to demonstrate this routine; this information is also printed with
111: the option -snes_view
112: */
113: SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf);
114: PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf);
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Initialize application:
118: Store right-hand-side of PDE and exact solution
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: xp = 0.0;
122: for (i = 0; i < n; i++) {
123: v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
124: VecSetValues(F, 1, &i, &v, INSERT_VALUES);
125: v = xp * xp * xp;
126: VecSetValues(U, 1, &i, &v, INSERT_VALUES);
127: xp += h;
128: }
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Evaluate initial guess; then solve nonlinear system
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: /*
134: Note: The user should initialize the vector, x, with the initial guess
135: for the nonlinear solver prior to calling SNESSolve(). In particular,
136: to employ an initial guess of zero, the user should explicitly set
137: this vector to zero by calling VecSet().
138: */
139: FormInitialGuess(x);
140: SNESSolve(snes, NULL, x);
141: SNESGetIterationNumber(snes, &its);
142: PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Check solution and clean up
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: /*
149: Check the error
150: */
151: VecAXPY(x, none, U);
152: VecNorm(x, NORM_2, &norm);
153: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
155: /*
156: Free work space. All PETSc objects should be destroyed when they
157: are no longer needed.
158: */
159: VecDestroy(&x);
160: VecDestroy(&r);
161: VecDestroy(&U);
162: VecDestroy(&F);
163: MatDestroy(&J);
164: SNESDestroy(&snes);
165: PetscFinalize();
166: return 0;
167: }
169: /*
170: FormInitialGuess - Computes initial guess.
172: Input/Output Parameter:
173: . x - the solution vector
174: */
175: PetscErrorCode FormInitialGuess(Vec x)
176: {
177: PetscScalar pfive = .50;
178: VecSet(x, pfive);
179: return 0;
180: }
182: /*
183: FormFunction - Evaluates nonlinear function, F(x).
185: Input Parameters:
186: . snes - the SNES context
187: . x - input vector
188: . ctx - optional user-defined context, as set by SNESSetFunction()
190: Output Parameter:
191: . f - function vector
193: Note:
194: The user-defined context can contain any application-specific data
195: needed for the function evaluation (such as various parameters, work
196: vectors, and grid information). In this program the context is just
197: a vector containing the right-hand-side of the discretized PDE.
198: */
200: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
201: {
202: Vec g = (Vec)ctx;
203: const PetscScalar *xx, *gg;
204: PetscScalar *ff, d;
205: PetscInt i, n;
207: /*
208: Get pointers to vector data.
209: - For default PETSc vectors, VecGetArray() returns a pointer to
210: the data array. Otherwise, the routine is implementation dependent.
211: - You MUST call VecRestoreArray() when you no longer need access to
212: the array.
213: */
214: VecGetArrayRead(x, &xx);
215: VecGetArray(f, &ff);
216: VecGetArrayRead(g, &gg);
218: /*
219: Compute function
220: */
221: VecGetSize(x, &n);
222: d = (PetscReal)(n - 1);
223: d = d * d;
224: ff[0] = xx[0];
225: for (i = 1; i < n - 1; i++) ff[i] = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - gg[i];
226: ff[n - 1] = xx[n - 1] - 1.0;
228: /*
229: Restore vectors
230: */
231: VecRestoreArrayRead(x, &xx);
232: VecRestoreArray(f, &ff);
233: VecRestoreArrayRead(g, &gg);
234: return 0;
235: }
237: /*
238: FormJacobian - Evaluates Jacobian matrix.
240: Input Parameters:
241: . snes - the SNES context
242: . x - input vector
243: . dummy - optional user-defined context (not used here)
245: Output Parameters:
246: . jac - Jacobian matrix
247: . B - optionally different preconditioning matrix
249: */
251: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
252: {
253: const PetscScalar *xx;
254: PetscScalar A[3], d;
255: PetscInt i, n, j[3];
257: /*
258: Get pointer to vector data
259: */
260: VecGetArrayRead(x, &xx);
262: /*
263: Compute Jacobian entries and insert into matrix.
264: - Note that in this case we set all elements for a particular
265: row at once.
266: */
267: VecGetSize(x, &n);
268: d = (PetscReal)(n - 1);
269: d = d * d;
271: /*
272: Interior grid points
273: */
274: for (i = 1; i < n - 1; i++) {
275: j[0] = i - 1;
276: j[1] = i;
277: j[2] = i + 1;
278: A[0] = A[2] = d;
279: A[1] = -2.0 * d + 2.0 * xx[i];
280: MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES);
281: }
283: /*
284: Boundary points
285: */
286: i = 0;
287: A[0] = 1.0;
289: MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES);
291: i = n - 1;
292: A[0] = 1.0;
294: MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES);
296: /*
297: Restore vector
298: */
299: VecRestoreArrayRead(x, &xx);
301: /*
302: Assemble matrix
303: */
304: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
305: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
306: if (jac != B) {
307: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
309: }
310: return 0;
311: }
313: /*TEST
315: testset:
316: args: -snes_monitor_short -snes_view -ksp_monitor
317: output_file: output/ex5_1.out
318: filter: grep -v "type: seqaij"
320: test:
321: suffix: 1
323: test:
324: suffix: cuda
325: requires: cuda
326: args: -mat_type aijcusparse -vec_type cuda
328: test:
329: suffix: kok
330: requires: kokkos_kernels
331: args: -mat_type aijkokkos -vec_type kokkos
333: # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly
334: # the solution is wrong on purpose
335: test:
336: requires: !single !complex
337: suffix: transpose_only
338: args: -snes_monitor_short -snes_view -ksp_monitor -snes_type ksptransposeonly -pc_type ilu -snes_test_jacobian -snes_test_jacobian_view -ksp_view_rhs -ksp_view_solution -ksp_view_mat_explicit -ksp_view_preconditioned_operator_explicit
340: TEST*/