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: }