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