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;
41: PetscFunctionBeginUser;
42: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
43: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
44: PetscCall(PetscOptionsGetBool(NULL, NULL, "-second_order", &second_order, NULL));
45: h = 1.0 / (n - 1);
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Create nonlinear solver context
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create vector data structures; set function evaluation routine
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
58: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
59: PetscCall(VecSetFromOptions(x));
60: PetscCall(VecDuplicate(x, &r));
61: PetscCall(VecDuplicate(x, &F));
63: PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
65: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: Create matrix data structures; set Jacobian evaluation routine
67: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69: PetscCall(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: /* PetscCall(SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0)); */
77: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, 0));
78: /* PetscCall(SNESSetJacobian(snes,J,JPrec,FormJacobian,0)); */
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Customize nonlinear solver; set runtime options
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: PetscCall(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: PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
97: v2 = a * xp + (b - a) * PetscPowScalar(xp, k);
98: PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES));
99: xp += h;
100: }
102: /* perturb initial guess */
103: PetscCall(VecGetArray(x, &xx));
104: for (i = 0; i < n; i++) {
105: v2 = xx[i] * sperturb;
106: PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES));
107: }
108: PetscCall(VecRestoreArray(x, &xx));
110: PetscCall(SNESSolve(snes, NULL, x));
111: PetscCall(SNESGetIterationNumber(snes, &it));
112: PetscCall(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: PetscCall(VecDestroy(&x));
120: PetscCall(VecDestroy(&r));
121: PetscCall(VecDestroy(&F));
122: PetscCall(MatDestroy(&J));
123: PetscCall(SNESDestroy(&snes));
124: PetscCall(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: PetscFunctionBeginUser;
135: PetscCall(VecGetArrayRead(x, &xx));
136: PetscCall(VecGetArray(f, &ff));
137: PetscCall(VecGetArray((Vec)dummy, &FF));
138: PetscCall(VecGetSize(x, &n));
139: d = (PetscReal)(n - 1);
140: d2 = d * d;
142: if (second_order) ff[0] = d * (0.5 * d * (-xx[2] + 4. * xx[1] - 3. * xx[0]) - X0DOT);
143: else ff[0] = d * (d * (xx[1] - xx[0]) - X0DOT);
145: 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];
147: ff[n - 1] = d * d * (xx[n - 1] - X1);
148: PetscCall(VecRestoreArrayRead(x, &xx));
149: PetscCall(VecRestoreArray(f, &ff));
150: PetscCall(VecRestoreArray((Vec)dummy, &FF));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat prejac, void *dummy)
155: {
156: const PetscScalar *xx;
157: PetscScalar A[3], d, d2;
158: PetscInt i, n, j[3];
160: PetscFunctionBeginUser;
161: PetscCall(VecGetSize(x, &n));
162: PetscCall(VecGetArrayRead(x, &xx));
163: d = (PetscReal)(n - 1);
164: d2 = d * d;
166: i = 0;
167: if (second_order) {
168: j[0] = 0;
169: j[1] = 1;
170: j[2] = 2;
171: A[0] = -3. * d * d * 0.5;
172: A[1] = 4. * d * d * 0.5;
173: A[2] = -1. * d * d * 0.5;
174: PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES));
175: } else {
176: j[0] = 0;
177: j[1] = 1;
178: A[0] = -d * d;
179: A[1] = d * d;
180: PetscCall(MatSetValues(prejac, 1, &i, 2, j, A, INSERT_VALUES));
181: }
182: for (i = 1; i < n - 1; i++) {
183: j[0] = i - 1;
184: j[1] = i;
185: j[2] = i + 1;
186: A[0] = d2;
187: A[1] = -2. * d2 + 2. * xx[i];
188: A[2] = d2;
189: PetscCall(MatSetValues(prejac, 1, &i, 3, j, A, INSERT_VALUES));
190: }
192: i = n - 1;
193: A[0] = d * d;
194: PetscCall(MatSetValues(prejac, 1, &i, 1, &i, &A[0], INSERT_VALUES));
196: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
197: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
198: PetscCall(MatAssemblyBegin(prejac, MAT_FINAL_ASSEMBLY));
199: PetscCall(MatAssemblyEnd(prejac, MAT_FINAL_ASSEMBLY));
201: PetscCall(VecRestoreArrayRead(x, &xx));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: /*TEST
207: test:
208: args: -n 14 -snes_monitor_short -snes_converged_reason
209: requires: !single
211: test:
212: suffix: 2
213: args: -n 15 -snes_monitor_short -snes_converged_reason
214: requires: !single
216: test:
217: suffix: 3
218: args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
219: requires: !single
221: TEST*/