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