Actual source code: ex59.c
petsc-3.3-p7 2013-05-11
2: #include <stdlib.h>
4: static const char help[] = "Tries to solve u`` + u^{2} = f for an easy case and an impossible case.\n\n";
6: /*
7: This example was contributed by Peter Graf to show how SNES fails when given a nonlinear problem with no solution.
9: Run with -n 14 to see fail to converge and -n 15 to see convergence
11: The option -second_order uses a different discretization of the Neumann boundary condition and always converges
12:
13: */
15: #include <petscsnes.h>
17: PetscBool second_order = PETSC_FALSE;
18: #define X0DOT -2.0
19: #define X1 5.0
20: #define KPOW 2.0
21: const PetscScalar sperturb = 1.1;
23: /*
24: User-defined routines
25: */
26: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
27: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
29: int main(int argc,char **argv)
30: {
31: SNES snes; /* SNES context */
32: Vec x,r,F; /* vectors */
33: Mat J; /* Jacobian */
34: PetscErrorCode ierr;
35: PetscInt it,n = 11,i;
36: PetscReal h,xp = 0.0;
37: PetscScalar v;
38: const PetscScalar a = X0DOT;
39: const PetscScalar b = X1;
40: const PetscScalar k = KPOW;
41: PetscScalar v2;
42: PetscScalar *xx;
43:
45: PetscInitialize(&argc,&argv,(char *)0,help);
46: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
47: PetscOptionsGetBool(PETSC_NULL,"-second_order",&second_order,PETSC_NULL);
48: h = 1.0/(n-1);
50: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: Create nonlinear solver context
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: SNESCreate(PETSC_COMM_WORLD,&snes);
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Create vector data structures; set function evaluation routine
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: VecCreate(PETSC_COMM_SELF,&x);
61: VecSetSizes(x,PETSC_DECIDE,n);
62: VecSetFromOptions(x);
63: VecDuplicate(x,&r);
64: VecDuplicate(x,&F);
66: SNESSetFunction(snes,r,FormFunction,(void*)F);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create matrix data structures; set Jacobian evaluation routine
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&J);
74: /*
75: Note that in this case we create separate matrices for the Jacobian
76: and preconditioner matrix. Both of these are computed in the
77: routine FormJacobian()
78: */
79: // SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0);
80: SNESSetJacobian(snes,J,J,FormJacobian,0);
81: /// SNESSetJacobian(snes,J,JPrec,FormJacobian,0);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Customize nonlinear solver; set runtime options
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: SNESSetFromOptions(snes);
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Initialize application:
91: Store right-hand-side of PDE and exact solution
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: /* set right hand side and initial guess to be exact solution of continuim problem */
95: #define SQR(x) ((x)*(x))
96: xp = 0.0;
97: for (i=0; i<n; i++)
98: {
99: 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.);
100: VecSetValues(F,1,&i,&v,INSERT_VALUES);
101: v2 = a*xp + (b-a)*PetscPowScalar(xp,k);
102: VecSetValues(x,1,&i,&v2,INSERT_VALUES);
103: xp += h;
104: }
106: /* perturb initial guess */
107: VecGetArray(x,&xx);
108: for (i=0; i<n; i++) {
109: v2 = xx[i]*sperturb;
110: VecSetValues(x,1,&i,&v2,INSERT_VALUES);
111: }
112: VecRestoreArray(x,&xx);
114:
115: SNESSolve(snes,PETSC_NULL,x);
116: SNESGetIterationNumber(snes,&it);
117: PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Free work space. All PETSc objects should be destroyed when they
121: are no longer needed.
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: VecDestroy(&x); VecDestroy(&r);
125: VecDestroy(&F); MatDestroy(&J);
126: SNESDestroy(&snes);
127: PetscFinalize();
129: return 0;
130: }
132: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
133: {
134: PetscScalar *xx,*ff,*FF,d,d2;
136: PetscInt i,n;
138: VecGetArray(x,&xx);
139: VecGetArray(f,&ff);
140: VecGetArray((Vec)dummy,&FF);
141: VecGetSize(x,&n);
142: d = (PetscReal)(n - 1); d2 = d*d;
144: if (second_order){
145: ff[0] = d*(0.5*d*(-xx[2] + 4.*xx[1] - 3.*xx[0]) - X0DOT);
146: } else {
147: ff[0] = d*(d*(xx[1] - xx[0]) - X0DOT);
148: }
149: for (i=1; i<n-1; i++) {
150: ff[i] = d2*(xx[i-1] - 2.*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
151: }
152: ff[n-1] = d*d*(xx[n-1] - X1);
153: VecRestoreArray(x,&xx);
154: VecRestoreArray(f,&ff);
155: VecRestoreArray((Vec)dummy,&FF);
156: return 0;
157: }
159: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *prejac,MatStructure *flag,void *dummy)
160: {
161: PetscScalar *xx,A[3],d,d2;
162: PetscInt i,n,j[3];
165: VecGetSize(x,&n);
166: VecGetArray(x,&xx);
167: d = (PetscReal)(n - 1); d2 = d*d;
169: i = 0;
170: if (second_order) {
171: j[0] = 0; j[1] = 1; j[2] = 2;
172: A[0] = -3.*d*d*0.5; A[1] = 4.*d*d*0.5; A[2] = -1.*d*d*0.5;
173: MatSetValues(*prejac,1,&i,3,j,A,INSERT_VALUES);
174: } else {
175: j[0] = 0; j[1] = 1;
176: A[0] = -d*d; A[1] = d*d;
177: MatSetValues(*prejac,1,&i,2,j,A,INSERT_VALUES);
178: }
179: for (i=1; i<n-1; i++) {
180: j[0] = i - 1; j[1] = i; j[2] = i + 1;
181: A[0] = d2; A[1] = -2.*d2 + 2.*xx[i]; A[2] = d2;
182: MatSetValues(*prejac,1,&i,3,j,A,INSERT_VALUES);
183: }
185: i = n-1;
186: A[0] = d*d;
187: MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
189: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
190: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
191: MatAssemblyBegin(*prejac,MAT_FINAL_ASSEMBLY);
192: MatAssemblyEnd(*prejac,MAT_FINAL_ASSEMBLY);
194: VecRestoreArray(x,&xx);
195: *flag = SAME_NONZERO_PATTERN;
196: return 0;
197: }