1 static const char help[] = "Tries to solve u`` + u^{2} = f for an easy case and an impossible case.\n\n";
2
3 /*
4 This example was contributed by Peter Graf to show how SNES fails when given a nonlinear problem with no solution.
5
6 Run with -n 14 to see fail to converge and -n 15 to see convergence
7
8 The option -second_order uses a different discretization of the Neumann boundary condition and always converges
9
10 */
11
12 #include <petscsnes.h>
13
14 PetscBool second_order = PETSC_FALSE;
15 #define X0DOT -2.0
16 #define X1 5.0
17 #define KPOW 2.0
18 const PetscScalar sperturb = 1.1;
19
20 /*
21 User-defined routines
22 */
23 PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
24 PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
25
main(int argc,char ** argv)26 int main(int argc, char **argv)
27 {
28 SNES snes; /* SNES context */
29 Vec x, r, F; /* vectors */
30 Mat J; /* Jacobian */
31 PetscInt it, n = 11, i;
32 PetscReal h, xp = 0.0;
33 PetscScalar v;
34 const PetscScalar a = X0DOT;
35 const PetscScalar b = X1;
36 const PetscScalar k = KPOW;
37 PetscScalar v2;
38 PetscScalar *xx;
39
40 PetscFunctionBeginUser;
41 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
42 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
43 PetscCall(PetscOptionsGetBool(NULL, NULL, "-second_order", &second_order, NULL));
44 h = 1.0 / (n - 1);
45
46 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47 Create nonlinear solver context
48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49
50 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
51
52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53 Create vector data structures; set function evaluation routine
54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55
56 PetscCall(VecCreate(PETSC_COMM_SELF, &x));
57 PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
58 PetscCall(VecSetFromOptions(x));
59 PetscCall(VecDuplicate(x, &r));
60 PetscCall(VecDuplicate(x, &F));
61
62 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
63
64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65 Create matrix data structures; set Jacobian evaluation routine
66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67
68 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, &J));
69
70 /*
71 Note that in this case we create separate matrices for the Jacobian
72 and matrix used to compute the preconditioner. Both of these are computed in the
73 routine FormJacobian()
74 */
75 /* PetscCall(SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0)); */
76 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, 0));
77 /* PetscCall(SNESSetJacobian(snes,J,JPrec,FormJacobian,0)); */
78
79 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80 Customize nonlinear solver; set runtime options
81 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82
83 PetscCall(SNESSetFromOptions(snes));
84
85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 Initialize application:
87 Store right-hand side of PDE and exact solution
88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89
90 /* set right-hand side and initial guess to be exact solution of continuum problem */
91 #define SQR(x) ((x) * (x))
92 xp = 0.0;
93 for (i = 0; i < n; i++) {
94 v = k * (k - 1.) * (b - a) * PetscPowScalar(xp, k - 2.) + SQR(a * xp) + SQR(b - a) * PetscPowScalar(xp, 2. * k) + 2. * a * (b - a) * PetscPowScalar(xp, k + 1.);
95 PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
96 v2 = a * xp + (b - a) * PetscPowScalar(xp, k);
97 PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES));
98 xp += h;
99 }
100
101 /* perturb initial guess */
102 PetscCall(VecGetArray(x, &xx));
103 for (i = 0; i < n; i++) {
104 v2 = xx[i] * sperturb;
105 PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES));
106 }
107 PetscCall(VecRestoreArray(x, &xx));
108
109 PetscCall(SNESSolve(snes, NULL, x));
110 PetscCall(SNESGetIterationNumber(snes, &it));
111 PetscCall(PetscPrintf(PETSC_COMM_SELF, "SNES iterations = %" PetscInt_FMT "\n\n", it));
112
113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114 Free work space. All PETSc objects should be destroyed when they
115 are no longer needed.
116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117
118 PetscCall(VecDestroy(&x));
119 PetscCall(VecDestroy(&r));
120 PetscCall(VecDestroy(&F));
121 PetscCall(MatDestroy(&J));
122 PetscCall(SNESDestroy(&snes));
123 PetscCall(PetscFinalize());
124 return 0;
125 }
126
FormFunction(SNES snes,Vec x,Vec f,void * dummy)127 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *dummy)
128 {
129 const PetscScalar *xx;
130 PetscScalar *ff, *FF, d, d2;
131 PetscInt i, n;
132
133 PetscFunctionBeginUser;
134 PetscCall(VecGetArrayRead(x, &xx));
135 PetscCall(VecGetArray(f, &ff));
136 PetscCall(VecGetArray((Vec)dummy, &FF));
137 PetscCall(VecGetSize(x, &n));
138 d = (PetscReal)(n - 1);
139 d2 = d * d;
140
141 if (second_order) ff[0] = d * (0.5 * d * (-xx[2] + 4. * xx[1] - 3. * xx[0]) - X0DOT);
142 else ff[0] = d * (d * (xx[1] - xx[0]) - X0DOT);
143
144 for (i = 1; i < n - 1; i++) ff[i] = d2 * (xx[i - 1] - 2. * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i];
145
146 ff[n - 1] = d * d * (xx[n - 1] - X1);
147 PetscCall(VecRestoreArrayRead(x, &xx));
148 PetscCall(VecRestoreArray(f, &ff));
149 PetscCall(VecRestoreArray((Vec)dummy, &FF));
150 PetscFunctionReturn(PETSC_SUCCESS);
151 }
152
FormJacobian(SNES snes,Vec x,Mat jac,Mat prejac,void * dummy)153 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat prejac, void *dummy)
154 {
155 const PetscScalar *xx;
156 PetscScalar A[3], d, d2;
157 PetscInt i, n, j[3];
158
159 PetscFunctionBeginUser;
160 PetscCall(VecGetSize(x, &n));
161 PetscCall(VecGetArrayRead(x, &xx));
162 d = (PetscReal)(n - 1);
163 d2 = d * d;
164
165 i = 0;
166 if (second_order) {
167 j[0] = 0;
168 j[1] = 1;
169 j[2] = 2;
170 A[0] = -3. * d * d * 0.5;
171 A[1] = 4. * d * d * 0.5;
172 A[2] = -1. * d * d * 0.5;
173 PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES));
174 } else {
175 j[0] = 0;
176 j[1] = 1;
177 A[0] = -d * d;
178 A[1] = d * d;
179 PetscCall(MatSetValues(prejac, 1, &i, 2, j, A, INSERT_VALUES));
180 }
181 for (i = 1; i < n - 1; i++) {
182 j[0] = i - 1;
183 j[1] = i;
184 j[2] = i + 1;
185 A[0] = d2;
186 A[1] = -2. * d2 + 2. * xx[i];
187 A[2] = d2;
188 PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES));
189 }
190
191 i = n - 1;
192 A[0] = d * d;
193 PetscCall(MatSetValues(prejac, 1, &i, 1, &i, &A[0], INSERT_VALUES));
194
195 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
196 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
197 PetscCall(MatAssemblyBegin(prejac, MAT_FINAL_ASSEMBLY));
198 PetscCall(MatAssemblyEnd(prejac, MAT_FINAL_ASSEMBLY));
199
200 PetscCall(VecRestoreArrayRead(x, &xx));
201 PetscFunctionReturn(PETSC_SUCCESS);
202 }
203
204 /*TEST
205
206 test:
207 args: -n 14 -snes_monitor_short -snes_converged_reason
208 requires: !single
209
210 test:
211 suffix: 2
212 args: -n 15 -snes_monitor_short -snes_converged_reason
213 requires: !single
214
215 test:
216 suffix: 3
217 args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
218 requires: !single
219
220 TEST*/
221