Actual source code: ex59.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  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);
 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);


112:   SNESSolve(snes,NULL,x);
113:   SNESGetIterationNumber(snes,&it);
114:   PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it);

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Free work space.  All PETSc objects should be destroyed when they
118:      are no longer needed.
119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   VecDestroy(&x);     VecDestroy(&r);
122:   VecDestroy(&F);     MatDestroy(&J);
123:   SNESDestroy(&snes);
124:   PetscFinalize();
125:   return ierr;
126: }

128: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
129: {
130:   const PetscScalar *xx;
131:   PetscScalar       *ff,*FF,d,d2;
132:   PetscErrorCode    ierr;
133:   PetscInt          i,n;

135:   VecGetArrayRead(x,&xx);
136:   VecGetArray(f,&ff);
137:   VecGetArray((Vec)dummy,&FF);
138:   VecGetSize(x,&n);
139:   d    = (PetscReal)(n - 1); 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:   VecRestoreArrayRead(x,&xx);
148:   VecRestoreArray(f,&ff);
149:   VecRestoreArray((Vec)dummy,&FF);
150:   return 0;
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];
158:   PetscErrorCode    ierr;

160:   VecGetSize(x,&n);
161:   VecGetArrayRead(x,&xx);
162:   d    = (PetscReal)(n - 1); d2 = d*d;

164:   i = 0;
165:   if (second_order) {
166:     j[0] = 0; j[1] = 1; j[2] = 2;
167:     A[0] = -3.*d*d*0.5; A[1] = 4.*d*d*0.5;  A[2] = -1.*d*d*0.5;
168:     MatSetValues(prejac,1,&i,3,j,A,INSERT_VALUES);
169:   } else {
170:     j[0] = 0; j[1] = 1;
171:     A[0] = -d*d; A[1] = d*d;
172:     MatSetValues(prejac,1,&i,2,j,A,INSERT_VALUES);
173:   }
174:   for (i=1; i<n-1; i++) {
175:     j[0] = i - 1; j[1] = i;                   j[2] = i + 1;
176:     A[0] = d2;    A[1] = -2.*d2 + 2.*xx[i];  A[2] = d2;
177:     MatSetValues(prejac,1,&i,3,j,A,INSERT_VALUES);
178:   }

180:   i    = n-1;
181:   A[0] = d*d;
182:   MatSetValues(prejac,1,&i,1,&i,&A[0],INSERT_VALUES);

184:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
185:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
186:   MatAssemblyBegin(prejac,MAT_FINAL_ASSEMBLY);
187:   MatAssemblyEnd(prejac,MAT_FINAL_ASSEMBLY);

189:   VecRestoreArrayRead(x,&xx);
190:   return 0;
191: }