Actual source code: ex1.c
1: static char help[] = "Newton's method for a two-variable system, sequential.\n\n";
3: /*
4: Include "petscsnes.h" so that we can use SNES solvers. Note that this
5: file automatically includes:
6: petscsys.h - base PETSc routines petscvec.h - vectors
7: petscmat.h - matrices
8: petscis.h - index sets petscksp.h - Krylov subspace methods
9: petscviewer.h - viewers petscpc.h - preconditioners
10: petscksp.h - linear solvers
11: */
12: /*F
13: This examples solves either
14: \begin{equation}
15: F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{x^2_0 + x_0 x_1 - 3}{x_0 x_1 + x^2_1 - 6}
16: \end{equation}
17: or if the {\tt -hard} options is given
18: \begin{equation}
19: F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
20: \end{equation}
21: F*/
22: #include <petscsnes.h>
24: /*
25: User-defined routines
26: */
27: extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
28: extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
29: extern PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *);
30: extern PetscErrorCode FormFunction2(SNES, Vec, Vec, void *);
32: int main(int argc, char **argv)
33: {
34: SNES snes; /* nonlinear solver context */
35: KSP ksp; /* linear solver context */
36: PC pc; /* preconditioner context */
37: Vec x, r; /* solution, residual vectors */
38: Mat J; /* Jacobian matrix */
39: PetscMPIInt size;
40: PetscScalar pfive = .5, *xx;
41: PetscBool flg;
43: PetscFunctionBeginUser;
44: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
45: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
46: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs");
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Create nonlinear solver context
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
52: PetscCall(SNESSetType(snes, SNESNEWTONLS));
53: PetscCall(SNESSetOptionsPrefix(snes, "mysolver_"));
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Create matrix and vector data structures; set corresponding routines
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: /*
59: Create vectors for solution and nonlinear function
60: */
61: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
62: PetscCall(VecSetSizes(x, PETSC_DECIDE, 2));
63: PetscCall(VecSetFromOptions(x));
64: PetscCall(VecDuplicate(x, &r));
66: /*
67: Create Jacobian matrix data structure
68: */
69: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
70: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
71: PetscCall(MatSetFromOptions(J));
72: PetscCall(MatSetUp(J));
74: PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg));
75: if (!flg) {
76: /*
77: Set function evaluation routine and vector.
78: */
79: PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));
81: /*
82: Set Jacobian matrix data structure and Jacobian evaluation routine
83: */
84: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
85: } else {
86: PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL));
87: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL));
88: }
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Customize nonlinear solver; set runtime options
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: /*
94: Set linear solver defaults for this problem. By extracting the
95: KSP and PC contexts from the SNES context, we can then
96: directly call any KSP and PC routines to set various options.
97: */
98: PetscCall(SNESGetKSP(snes, &ksp));
99: PetscCall(KSPGetPC(ksp, &pc));
100: PetscCall(PCSetType(pc, PCNONE));
101: PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_CURRENT, PETSC_CURRENT, 20));
103: /*
104: Set SNES/KSP/KSP/PC runtime options, e.g.,
105: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
106: These options will override those specified above as long as
107: SNESSetFromOptions() is called _after_ any other customization
108: routines.
109: */
110: PetscCall(SNESSetFromOptions(snes));
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Evaluate initial guess; then solve nonlinear system
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: if (!flg) PetscCall(VecSet(x, pfive));
116: else {
117: PetscCall(VecGetArray(x, &xx));
118: xx[0] = 2.0;
119: xx[1] = 3.0;
120: PetscCall(VecRestoreArray(x, &xx));
121: }
122: /*
123: Note: The user should initialize the vector, x, with the initial guess
124: for the nonlinear solver prior to calling SNESSolve(). In particular,
125: to employ an initial guess of zero, the user should explicitly set
126: this vector to zero by calling VecSet().
127: */
129: PetscCall(SNESSolve(snes, NULL, x));
130: if (flg) {
131: Vec f;
132: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
133: PetscCall(SNESGetFunction(snes, &f, 0, 0));
134: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
135: }
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Free work space. All PETSc objects should be destroyed when they
139: are no longer needed.
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: PetscCall(VecDestroy(&x));
143: PetscCall(VecDestroy(&r));
144: PetscCall(MatDestroy(&J));
145: PetscCall(SNESDestroy(&snes));
146: PetscCall(PetscFinalize());
147: return 0;
148: }
149: /* ------------------------------------------------------------------- */
150: /*
151: FormFunction1 - Evaluates nonlinear function, F(x).
153: Input Parameters:
154: . snes - the SNES context
155: . x - input vector
156: . ctx - optional user-defined context
158: Output Parameter:
159: . f - function vector
160: */
161: PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, PetscCtx ctx)
162: {
163: const PetscScalar *xx;
164: PetscScalar *ff;
166: PetscFunctionBeginUser;
167: /*
168: Get pointers to vector data.
169: - For default PETSc vectors, VecGetArray() returns a pointer to
170: the data array. Otherwise, the routine is implementation dependent.
171: - You MUST call VecRestoreArray() when you no longer need access to
172: the array.
173: */
174: PetscCall(VecGetArrayRead(x, &xx));
175: PetscCall(VecGetArray(f, &ff));
177: /* Compute function */
178: ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0;
179: ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0;
181: /* Restore vectors */
182: PetscCall(VecRestoreArrayRead(x, &xx));
183: PetscCall(VecRestoreArray(f, &ff));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
186: /* ------------------------------------------------------------------- */
187: /*
188: FormJacobian1 - Evaluates Jacobian matrix.
190: Input Parameters:
191: . snes - the SNES context
192: . x - input vector
193: . dummy - optional user-defined context (not used here)
195: Output Parameters:
196: . jac - Jacobian matrix
197: . B - optionally different matrix used to construct the preconditioner
199: */
200: PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
201: {
202: const PetscScalar *xx;
203: PetscScalar A[4];
204: PetscInt idx[2] = {0, 1};
206: PetscFunctionBeginUser;
207: /*
208: Get pointer to vector data
209: */
210: PetscCall(VecGetArrayRead(x, &xx));
212: /*
213: Compute Jacobian entries and insert into matrix.
214: - Since this is such a small problem, we set all entries for
215: the matrix at once.
216: */
217: A[0] = 2.0 * xx[0] + xx[1];
218: A[1] = xx[0];
219: A[2] = xx[1];
220: A[3] = xx[0] + 2.0 * xx[1];
221: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
223: /*
224: Restore vector
225: */
226: PetscCall(VecRestoreArrayRead(x, &xx));
228: /*
229: Assemble matrix
230: */
231: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
232: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
233: if (jac != B) {
234: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
235: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
236: }
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /* ------------------------------------------------------------------- */
241: PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
242: {
243: const PetscScalar *xx;
244: PetscScalar *ff;
246: PetscFunctionBeginUser;
247: /*
248: Get pointers to vector data.
249: - For default PETSc vectors, VecGetArray() returns a pointer to
250: the data array. Otherwise, the routine is implementation dependent.
251: - You MUST call VecRestoreArray() when you no longer need access to
252: the array.
253: */
254: PetscCall(VecGetArrayRead(x, &xx));
255: PetscCall(VecGetArray(f, &ff));
257: /*
258: Compute function
259: */
260: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
261: ff[1] = xx[1];
263: /*
264: Restore vectors
265: */
266: PetscCall(VecRestoreArrayRead(x, &xx));
267: PetscCall(VecRestoreArray(f, &ff));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
270: /* ------------------------------------------------------------------- */
271: PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
272: {
273: const PetscScalar *xx;
274: PetscScalar A[4];
275: PetscInt idx[2] = {0, 1};
277: PetscFunctionBeginUser;
278: /*
279: Get pointer to vector data
280: */
281: PetscCall(VecGetArrayRead(x, &xx));
283: /*
284: Compute Jacobian entries and insert into matrix.
285: - Since this is such a small problem, we set all entries for
286: the matrix at once.
287: */
288: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
289: A[1] = 0.0;
290: A[2] = 0.0;
291: A[3] = 1.0;
292: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
294: /*
295: Restore vector
296: */
297: PetscCall(VecRestoreArrayRead(x, &xx));
299: /*
300: Assemble matrix
301: */
302: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
303: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
304: if (jac != B) {
305: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
306: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
307: }
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*TEST
313: test:
314: args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop
315: requires: !single
317: # test harness puts {{ }} options always at the end, need to specify the prefix explicitly
318: test:
319: suffix: 2
320: requires: !single
321: args: -prefix_push mysolver_ -snes_monitor_short -prefix_pop -mysolver_snes_ksp_ew {{0 1}}
322: output_file: output/ex1_1.out
324: test:
325: suffix: 3
326: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short -prefix_pop
327: requires: !single
328: output_file: output/ex1_1.out
330: test:
331: suffix: 4
332: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short -prefix_pop
333: requires: !single
334: output_file: output/ex1_1.out
336: test:
337: suffix: 5
338: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short -prefix_pop
339: requires: !single
340: output_file: output/ex1_1.out
342: test:
343: suffix: 6
344: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short -prefix_pop
345: requires: !single
346: output_file: output/ex1_1.out
348: test:
349: suffix: X
350: args: -prefix_push mysolver_ -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 -prefix_pop
351: requires: !single x
353: TEST*/