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;
45: PetscInitialize(&argc, &argv, (char *)0, help);
46: MPI_Comm_size(PETSC_COMM_WORLD, &size);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create nonlinear solver context
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: SNESCreate(PETSC_COMM_WORLD, &snes);
53: SNESSetType(snes, SNESNEWTONLS);
54: 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: VecCreate(PETSC_COMM_WORLD, &x);
63: VecSetSizes(x, PETSC_DECIDE, 2);
64: VecSetFromOptions(x);
65: VecDuplicate(x, &r);
67: /*
68: Create Jacobian matrix data structure
69: */
70: MatCreate(PETSC_COMM_WORLD, &J);
71: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2);
72: MatSetFromOptions(J);
73: MatSetUp(J);
75: PetscOptionsHasName(NULL, NULL, "-hard", &flg);
76: if (!flg) {
77: /*
78: Set function evaluation routine and vector.
79: */
80: SNESSetFunction(snes, r, FormFunction1, NULL);
82: /*
83: Set Jacobian matrix data structure and Jacobian evaluation routine
84: */
85: SNESSetJacobian(snes, J, J, FormJacobian1, NULL);
86: } else {
87: SNESSetFunction(snes, r, FormFunction2, NULL);
88: 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: SNESGetKSP(snes, &ksp);
100: KSPGetPC(ksp, &pc);
101: PCSetType(pc, PCNONE);
102: 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: SNESSetFromOptions(snes);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Evaluate initial guess; then solve nonlinear system
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: if (!flg) {
117: VecSet(x, pfive);
118: } else {
119: VecGetArray(x, &xx);
120: xx[0] = 2.0;
121: xx[1] = 3.0;
122: 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: SNESSolve(snes, NULL, x);
132: if (flg) {
133: Vec f;
134: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
135: SNESGetFunction(snes, &f, 0, 0);
136: 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: VecDestroy(&x);
145: VecDestroy(&r);
146: MatDestroy(&J);
147: SNESDestroy(&snes);
148: 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: /*
169: Get pointers to vector data.
170: - For default PETSc vectors, VecGetArray() returns a pointer to
171: the data array. Otherwise, the routine is implementation dependent.
172: - You MUST call VecRestoreArray() when you no longer need access to
173: the array.
174: */
175: VecGetArrayRead(x, &xx);
176: VecGetArray(f, &ff);
178: /* Compute function */
179: ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0;
180: ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0;
182: /* Restore vectors */
183: VecRestoreArrayRead(x, &xx);
184: VecRestoreArray(f, &ff);
185: return 0;
186: }
187: /* ------------------------------------------------------------------- */
188: /*
189: FormJacobian1 - Evaluates Jacobian matrix.
191: Input Parameters:
192: . snes - the SNES context
193: . x - input vector
194: . dummy - optional user-defined context (not used here)
196: Output Parameters:
197: . jac - Jacobian matrix
198: . B - optionally different preconditioning matrix
199: . flag - flag indicating matrix structure
200: */
201: PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
202: {
203: const PetscScalar *xx;
204: PetscScalar A[4];
205: PetscInt idx[2] = {0, 1};
207: /*
208: Get pointer to vector data
209: */
210: 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: MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES);
223: /*
224: Restore vector
225: */
226: VecRestoreArrayRead(x, &xx);
228: /*
229: Assemble matrix
230: */
231: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
232: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
233: if (jac != B) {
234: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
235: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
236: }
237: return 0;
238: }
240: /* ------------------------------------------------------------------- */
241: PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
242: {
243: const PetscScalar *xx;
244: PetscScalar *ff;
246: /*
247: Get pointers to vector data.
248: - For default PETSc vectors, VecGetArray() returns a pointer to
249: the data array. Otherwise, the routine is implementation dependent.
250: - You MUST call VecRestoreArray() when you no longer need access to
251: the array.
252: */
253: VecGetArrayRead(x, &xx);
254: VecGetArray(f, &ff);
256: /*
257: Compute function
258: */
259: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
260: ff[1] = xx[1];
262: /*
263: Restore vectors
264: */
265: VecRestoreArrayRead(x, &xx);
266: VecRestoreArray(f, &ff);
267: return 0;
268: }
269: /* ------------------------------------------------------------------- */
270: PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
271: {
272: const PetscScalar *xx;
273: PetscScalar A[4];
274: PetscInt idx[2] = {0, 1};
276: /*
277: Get pointer to vector data
278: */
279: VecGetArrayRead(x, &xx);
281: /*
282: Compute Jacobian entries and insert into matrix.
283: - Since this is such a small problem, we set all entries for
284: the matrix at once.
285: */
286: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
287: A[1] = 0.0;
288: A[2] = 0.0;
289: A[3] = 1.0;
290: MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES);
292: /*
293: Restore vector
294: */
295: VecRestoreArrayRead(x, &xx);
297: /*
298: Assemble matrix
299: */
300: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
301: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
302: if (jac != B) {
303: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
304: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
305: }
306: return 0;
307: }
309: /*TEST
311: test:
312: args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop
313: requires: !single
315: # test harness puts {{ }} options always at the end, need to specify the prefix explicitly
316: test:
317: suffix: 2
318: requires: !single
319: args: -prefix_push mysolver_ -snes_monitor_short -prefix_pop -mysolver_snes_ksp_ew {{0 1}}
320: output_file: output/ex1_1.out
322: test:
323: suffix: 3
324: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short -prefix_pop
325: requires: !single
326: output_file: output/ex1_1.out
328: test:
329: suffix: 4
330: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short -prefix_pop
331: requires: !single
332: output_file: output/ex1_1.out
334: test:
335: suffix: 5
336: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short -prefix_pop
337: requires: !single
338: output_file: output/ex1_1.out
340: test:
341: suffix: 6
342: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short -prefix_pop
343: requires: !single
344: output_file: output/ex1_1.out
346: test:
347: suffix: X
348: args: -prefix_push mysolver_ -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 -prefix_pop
349: requires: !single x
351: TEST*/