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) {
116: PetscCall(VecSet(x, pfive));
117: } else {
118: PetscCall(VecGetArray(x, &xx));
119: xx[0] = 2.0;
120: xx[1] = 3.0;
121: PetscCall(VecRestoreArray(x, &xx));
122: }
123: /*
124: Note: The user should initialize the vector, x, with the initial guess
125: for the nonlinear solver prior to calling SNESSolve(). In particular,
126: to employ an initial guess of zero, the user should explicitly set
127: this vector to zero by calling VecSet().
128: */
130: PetscCall(SNESSolve(snes, NULL, x));
131: if (flg) {
132: Vec f;
133: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
134: PetscCall(SNESGetFunction(snes, &f, 0, 0));
135: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
136: }
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Free work space. All PETSc objects should be destroyed when they
140: are no longer needed.
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: PetscCall(VecDestroy(&x));
144: PetscCall(VecDestroy(&r));
145: PetscCall(MatDestroy(&J));
146: PetscCall(SNESDestroy(&snes));
147: PetscCall(PetscFinalize());
148: return 0;
149: }
150: /* ------------------------------------------------------------------- */
151: /*
152: FormFunction1 - Evaluates nonlinear function, F(x).
154: Input Parameters:
155: . snes - the SNES context
156: . x - input vector
157: . ctx - optional user-defined context
159: Output Parameter:
160: . f - function vector
161: */
162: PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx)
163: {
164: const PetscScalar *xx;
165: PetscScalar *ff;
167: PetscFunctionBeginUser;
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: PetscCall(VecGetArrayRead(x, &xx));
176: PetscCall(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: PetscCall(VecRestoreArrayRead(x, &xx));
184: PetscCall(VecRestoreArray(f, &ff));
185: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBeginUser;
208: /*
209: Get pointer to vector data
210: */
211: PetscCall(VecGetArrayRead(x, &xx));
213: /*
214: Compute Jacobian entries and insert into matrix.
215: - Since this is such a small problem, we set all entries for
216: the matrix at once.
217: */
218: A[0] = 2.0 * xx[0] + xx[1];
219: A[1] = xx[0];
220: A[2] = xx[1];
221: A[3] = xx[0] + 2.0 * xx[1];
222: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
224: /*
225: Restore vector
226: */
227: PetscCall(VecRestoreArrayRead(x, &xx));
229: /*
230: Assemble matrix
231: */
232: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
233: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
234: if (jac != B) {
235: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
236: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
237: }
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /* ------------------------------------------------------------------- */
242: PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
243: {
244: const PetscScalar *xx;
245: PetscScalar *ff;
247: PetscFunctionBeginUser;
248: /*
249: Get pointers to vector data.
250: - For default PETSc vectors, VecGetArray() returns a pointer to
251: the data array. Otherwise, the routine is implementation dependent.
252: - You MUST call VecRestoreArray() when you no longer need access to
253: the array.
254: */
255: PetscCall(VecGetArrayRead(x, &xx));
256: PetscCall(VecGetArray(f, &ff));
258: /*
259: Compute function
260: */
261: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
262: ff[1] = xx[1];
264: /*
265: Restore vectors
266: */
267: PetscCall(VecRestoreArrayRead(x, &xx));
268: PetscCall(VecRestoreArray(f, &ff));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
271: /* ------------------------------------------------------------------- */
272: PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
273: {
274: const PetscScalar *xx;
275: PetscScalar A[4];
276: PetscInt idx[2] = {0, 1};
278: PetscFunctionBeginUser;
279: /*
280: Get pointer to vector data
281: */
282: PetscCall(VecGetArrayRead(x, &xx));
284: /*
285: Compute Jacobian entries and insert into matrix.
286: - Since this is such a small problem, we set all entries for
287: the matrix at once.
288: */
289: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
290: A[1] = 0.0;
291: A[2] = 0.0;
292: A[3] = 1.0;
293: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
295: /*
296: Restore vector
297: */
298: PetscCall(VecRestoreArrayRead(x, &xx));
300: /*
301: Assemble matrix
302: */
303: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
304: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
305: if (jac != B) {
306: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
307: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
308: }
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*TEST
314: test:
315: args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop
316: requires: !single
318: # test harness puts {{ }} options always at the end, need to specify the prefix explicitly
319: test:
320: suffix: 2
321: requires: !single
322: args: -prefix_push mysolver_ -snes_monitor_short -prefix_pop -mysolver_snes_ksp_ew {{0 1}}
323: output_file: output/ex1_1.out
325: test:
326: suffix: 3
327: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short -prefix_pop
328: requires: !single
329: output_file: output/ex1_1.out
331: test:
332: suffix: 4
333: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short -prefix_pop
334: requires: !single
335: output_file: output/ex1_1.out
337: test:
338: suffix: 5
339: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short -prefix_pop
340: requires: !single
341: output_file: output/ex1_1.out
343: test:
344: suffix: 6
345: args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short -prefix_pop
346: requires: !single
347: output_file: output/ex1_1.out
349: test:
350: suffix: X
351: args: -prefix_push mysolver_ -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 -prefix_pop
352: requires: !single x
354: TEST*/