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;
36: PetscFunctionBeginUser;
37: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
38: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
39: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
40: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
41: h = 1.0 / (n - 1);
43: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Create nonlinear solver context
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
48: PetscCall(SNESGetKSP(snes, &ksp));
49: PetscCall(KSPGetPC(ksp, &pc));
50: PetscCall(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: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
58: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
59: PetscCall(VecSetFromOptions(x));
60: PetscCall(VecDuplicate(x, &r));
61: PetscCall(VecDuplicate(x, &F));
62: PetscCall(VecDuplicate(x, &U));
64: /*
65: Set function evaluation routine and vector
66: */
67: PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Create matrix data structure; set Jacobian evaluation routine
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
74: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
75: PetscCall(MatSetFromOptions(J));
76: PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));
77: PetscCall(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: PetscCall(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: PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
100: PetscCall(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: PetscCall(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: PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
114: 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));
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: PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
125: v = xp * xp * xp;
126: PetscCall(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: PetscCall(FormInitialGuess(x));
140: PetscCall(SNESSolve(snes, NULL, x));
141: PetscCall(SNESGetIterationNumber(snes, &its));
142: PetscCall(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: PetscCall(VecAXPY(x, none, U));
152: PetscCall(VecNorm(x, NORM_2, &norm));
153: PetscCall(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: PetscCall(VecDestroy(&x));
160: PetscCall(VecDestroy(&r));
161: PetscCall(VecDestroy(&U));
162: PetscCall(VecDestroy(&F));
163: PetscCall(MatDestroy(&J));
164: PetscCall(SNESDestroy(&snes));
165: PetscCall(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;
179: PetscFunctionBeginUser;
180: PetscCall(VecSet(x, pfive));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*
185: FormFunction - Evaluates nonlinear function, F(x).
187: Input Parameters:
188: . snes - the SNES context
189: . x - input vector
190: . ctx - optional user-defined context, as set by SNESSetFunction()
192: Output Parameter:
193: . f - function vector
195: Note:
196: The user-defined context can contain any application-specific data
197: needed for the function evaluation (such as various parameters, work
198: vectors, and grid information). In this program the context is just
199: a vector containing the right-hand-side of the discretized PDE.
200: */
202: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
203: {
204: Vec g = (Vec)ctx;
205: const PetscScalar *xx, *gg;
206: PetscScalar *ff, d;
207: PetscInt i, n;
209: PetscFunctionBeginUser;
210: /*
211: Get pointers to vector data.
212: - For default PETSc vectors, VecGetArray() returns a pointer to
213: the data array. Otherwise, the routine is implementation dependent.
214: - You MUST call VecRestoreArray() when you no longer need access to
215: the array.
216: */
217: PetscCall(VecGetArrayRead(x, &xx));
218: PetscCall(VecGetArray(f, &ff));
219: PetscCall(VecGetArrayRead(g, &gg));
221: /*
222: Compute function
223: */
224: PetscCall(VecGetSize(x, &n));
225: d = (PetscReal)(n - 1);
226: d = d * d;
227: ff[0] = xx[0];
228: 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];
229: ff[n - 1] = xx[n - 1] - 1.0;
231: /*
232: Restore vectors
233: */
234: PetscCall(VecRestoreArrayRead(x, &xx));
235: PetscCall(VecRestoreArray(f, &ff));
236: PetscCall(VecRestoreArrayRead(g, &gg));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*
241: FormJacobian - Evaluates Jacobian matrix.
243: Input Parameters:
244: . snes - the SNES context
245: . x - input vector
246: . dummy - optional user-defined context (not used here)
248: Output Parameters:
249: . jac - Jacobian matrix
250: . B - optionally different preconditioning matrix
252: */
254: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
255: {
256: const PetscScalar *xx;
257: PetscScalar A[3], d;
258: PetscInt i, n, j[3];
260: PetscFunctionBeginUser;
261: /*
262: Get pointer to vector data
263: */
264: PetscCall(VecGetArrayRead(x, &xx));
266: /*
267: Compute Jacobian entries and insert into matrix.
268: - Note that in this case we set all elements for a particular
269: row at once.
270: */
271: PetscCall(VecGetSize(x, &n));
272: d = (PetscReal)(n - 1);
273: d = d * d;
275: /*
276: Interior grid points
277: */
278: for (i = 1; i < n - 1; i++) {
279: j[0] = i - 1;
280: j[1] = i;
281: j[2] = i + 1;
282: A[0] = d;
283: A[1] = -2.0 * d + 2.0 * xx[i];
284: A[2] = d;
285: PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
286: }
288: /*
289: Boundary points
290: */
291: i = 0;
292: A[0] = 1.0;
294: PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
296: i = n - 1;
297: A[0] = 1.0;
299: PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
301: /*
302: Restore vector
303: */
304: PetscCall(VecRestoreArrayRead(x, &xx));
306: /*
307: Assemble matrix
308: */
309: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
310: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
311: if (jac != B) {
312: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
313: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
314: }
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*TEST
320: testset:
321: args: -snes_monitor_short -snes_view -ksp_monitor
322: output_file: output/ex5_1.out
323: filter: grep -v "type: seqaij"
325: test:
326: suffix: 1
328: test:
329: suffix: cuda
330: requires: cuda
331: args: -mat_type aijcusparse -vec_type cuda
333: test:
334: suffix: kok
335: requires: kokkos_kernels
336: args: -mat_type aijkokkos -vec_type kokkos
338: # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly
339: # the solution is wrong on purpose
340: test:
341: requires: !single !complex
342: suffix: transpose_only
343: 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
345: TEST*/