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