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
7: petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets
10: 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, NULL, 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_CURRENT, PETSC_CURRENT, 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) 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, NULL, NULL));
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, PetscCtx 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 matrix used to construct the preconditioner
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_orthogonalization_cgs_refinement_type refine_always -snes_monitor -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 -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 -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 -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 -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 -prefix_pop
346: requires: !single
347: output_file: output/ex1_1.out
349: test:
350: suffix: X
351: args: -prefix_push mysolver_ -ksp_monitor -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor -snes_rtol 1.e-4 -prefix_pop
352: requires: !single x
354: TEST*/