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: PetscErrorCode ierr;
33: PetscInt it,n = 11,i;
34: PetscReal h,xp = 0.0;
35: PetscScalar v;
36: const PetscScalar a = X0DOT;
37: const PetscScalar b = X1;
38: const PetscScalar k = KPOW;
39: PetscScalar v2;
40: PetscScalar *xx;
42: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
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: {
96: 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.);
97: VecSetValues(F,1,&i,&v,INSERT_VALUES);
98: v2 = a*xp + (b-a)*PetscPowScalar(xp,k);
99: VecSetValues(x,1,&i,&v2,INSERT_VALUES);
100: xp += h;
101: }
103: /* perturb initial guess */
104: VecGetArray(x,&xx);
105: for (i=0; i<n; i++) {
106: v2 = xx[i]*sperturb;
107: VecSetValues(x,1,&i,&v2,INSERT_VALUES);
108: }
109: VecRestoreArray(x,&xx);
111: SNESSolve(snes,NULL,x);
112: SNESGetIterationNumber(snes,&it);
113: PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Free work space. All PETSc objects should be destroyed when they
117: are no longer needed.
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: VecDestroy(&x); VecDestroy(&r);
121: VecDestroy(&F); MatDestroy(&J);
122: SNESDestroy(&snes);
123: PetscFinalize();
124: return ierr;
125: }
127: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
128: {
129: const PetscScalar *xx;
130: PetscScalar *ff,*FF,d,d2;
131: PetscErrorCode ierr;
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); d2 = d*d;
140: if (second_order) ff[0] = d*(0.5*d*(-xx[2] + 4.*xx[1] - 3.*xx[0]) - X0DOT);
141: else ff[0] = d*(d*(xx[1] - xx[0]) - X0DOT);
143: 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: ff[n-1] = d*d*(xx[n-1] - X1);
146: VecRestoreArrayRead(x,&xx);
147: VecRestoreArray(f,&ff);
148: VecRestoreArray((Vec)dummy,&FF);
149: return 0;
150: }
152: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat prejac,void *dummy)
153: {
154: const PetscScalar *xx;
155: PetscScalar A[3],d,d2;
156: PetscInt i,n,j[3];
157: PetscErrorCode ierr;
159: VecGetSize(x,&n);
160: VecGetArrayRead(x,&xx);
161: d = (PetscReal)(n - 1); d2 = d*d;
163: i = 0;
164: if (second_order) {
165: j[0] = 0; j[1] = 1; j[2] = 2;
166: A[0] = -3.*d*d*0.5; A[1] = 4.*d*d*0.5; A[2] = -1.*d*d*0.5;
167: MatSetValues(prejac,1,&i,3,j,A,INSERT_VALUES);
168: } else {
169: j[0] = 0; j[1] = 1;
170: A[0] = -d*d; A[1] = d*d;
171: MatSetValues(prejac,1,&i,2,j,A,INSERT_VALUES);
172: }
173: for (i=1; i<n-1; i++) {
174: j[0] = i - 1; j[1] = i; j[2] = i + 1;
175: A[0] = d2; A[1] = -2.*d2 + 2.*xx[i]; A[2] = d2;
176: MatSetValues(prejac,1,&i,3,j,A,INSERT_VALUES);
177: }
179: i = n-1;
180: A[0] = d*d;
181: MatSetValues(prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
183: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
184: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
185: MatAssemblyBegin(prejac,MAT_FINAL_ASSEMBLY);
186: MatAssemblyEnd(prejac,MAT_FINAL_ASSEMBLY);
188: VecRestoreArrayRead(x,&xx);
189: return 0;
190: }
192: /*TEST
194: test:
195: args: -n 14 -snes_monitor_short -snes_converged_reason
196: requires: !single
198: test:
199: suffix: 2
200: args: -n 15 -snes_monitor_short -snes_converged_reason
201: requires: !single
203: test:
204: suffix: 3
205: args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
206: requires: !single
208: TEST*/