1 static char help[] = "Newton's method for a two-variable system, sequential.\n\n";
2
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>
23
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 *);
31
main(int argc,char ** argv)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;
42
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");
47
48 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49 Create nonlinear solver context
50 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
52 PetscCall(SNESSetType(snes, SNESNEWTONLS));
53 PetscCall(SNESSetOptionsPrefix(snes, "mysolver_"));
54
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));
65
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));
73
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));
80
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 }
89
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));
102
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));
111
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 */
129
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 }
137
138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139 Free work space. All PETSc objects should be destroyed when they
140 are no longer needed.
141 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142
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).
153
154 Input Parameters:
155 . snes - the SNES context
156 . x - input vector
157 . ctx - optional user-defined context
158
159 Output Parameter:
160 . f - function vector
161 */
FormFunction1(SNES snes,Vec x,Vec f,PetscCtx ctx)162 PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, PetscCtx ctx)
163 {
164 const PetscScalar *xx;
165 PetscScalar *ff;
166
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));
177
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;
181
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.
190
191 Input Parameters:
192 . snes - the SNES context
193 . x - input vector
194 . dummy - optional user-defined context (not used here)
195
196 Output Parameters:
197 . jac - Jacobian matrix
198 . B - optionally different matrix used to construct the preconditioner
199
200 */
FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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};
206
207 PetscFunctionBeginUser;
208 /*
209 Get pointer to vector data
210 */
211 PetscCall(VecGetArrayRead(x, &xx));
212
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));
223
224 /*
225 Restore vector
226 */
227 PetscCall(VecRestoreArrayRead(x, &xx));
228
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 }
240
241 /* ------------------------------------------------------------------- */
FormFunction2(SNES snes,Vec x,Vec f,void * dummy)242 PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
243 {
244 const PetscScalar *xx;
245 PetscScalar *ff;
246
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));
257
258 /*
259 Compute function
260 */
261 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
262 ff[1] = xx[1];
263
264 /*
265 Restore vectors
266 */
267 PetscCall(VecRestoreArrayRead(x, &xx));
268 PetscCall(VecRestoreArray(f, &ff));
269 PetscFunctionReturn(PETSC_SUCCESS);
270 }
271 /* ------------------------------------------------------------------- */
FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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};
277
278 PetscFunctionBeginUser;
279 /*
280 Get pointer to vector data
281 */
282 PetscCall(VecGetArrayRead(x, &xx));
283
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));
294
295 /*
296 Restore vector
297 */
298 PetscCall(VecRestoreArrayRead(x, &xx));
299
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 }
311
312 /*TEST
313
314 test:
315 args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop
316 requires: !single
317
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
324
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
330
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
336
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
342
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
348
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
353
354 TEST*/
355