Actual source code: ex59.c
1: static const char help[] = "Tries to solve u`` + u^{2} = f for an easy case and an impossible case.\n\n";
3: /*
4: This example was contributed by Peter Graf to show how SNES fails when given a nonlinear problem with no solution.
6: Run with -n 14 to see fail to converge and -n 15 to see convergence
8: The option -second_order uses a different discretization of the Neumann boundary condition and always converges
10: */
12: #include <petscsnes.h>
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;
20: /*
21: User-defined routines
22: */
23: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
24: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
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;
40: PetscFunctionBeginUser;
41: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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);
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Create nonlinear solver context
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Create vector data structures; set function evaluation routine
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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));
62: PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
64: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: Create matrix data structures; set Jacobian evaluation routine
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, &J));
70: /*
71: Note that in this case we create separate matrices for the Jacobian
72: and preconditioner matrix. 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)); */
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Customize nonlinear solver; set runtime options
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83: PetscCall(SNESSetFromOptions(snes));
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Initialize application:
87: Store right-hand side of PDE and exact solution
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: /* set right-hand side and initial guess to be exact solution of continuim 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: }
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));
109: PetscCall(SNESSolve(snes, NULL, x));
110: PetscCall(SNESGetIterationNumber(snes, &it));
111: PetscCall(PetscPrintf(PETSC_COMM_SELF, "SNES iterations = %" PetscInt_FMT "\n\n", it));
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Free work space. All PETSc objects should be destroyed when they
115: are no longer needed.
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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: }
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;
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;
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: PetscCall(VecRestoreArrayRead(x, &xx));
148: PetscCall(VecRestoreArray(f, &ff));
149: PetscCall(VecRestoreArray((Vec)dummy, &FF));
150: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBeginUser;
160: PetscCall(VecGetSize(x, &n));
161: PetscCall(VecGetArrayRead(x, &xx));
162: d = (PetscReal)(n - 1);
163: d2 = d * d;
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: }
191: i = n - 1;
192: A[0] = d * d;
193: PetscCall(MatSetValues(prejac, 1, &i, 1, &i, &A[0], INSERT_VALUES));
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));
200: PetscCall(VecRestoreArrayRead(x, &xx));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*TEST
206: test:
207: args: -n 14 -snes_monitor_short -snes_converged_reason
208: requires: !single
210: test:
211: suffix: 2
212: args: -n 15 -snes_monitor_short -snes_converged_reason
213: requires: !single
215: test:
216: suffix: 3
217: args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
218: requires: !single
220: TEST*/