Actual source code: ex59.c

petsc-3.6.4 2016-04-12
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;


 43:   PetscInitialize(&argc,&argv,(char*)0,help);
 44:   PetscOptionsGetInt(NULL,"-n",&n,NULL);
 45:   PetscOptionsGetBool(NULL,"-second_order",&second_order,NULL);
 46:   h    = 1.0/(n-1);

 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:      Create nonlinear solver context
 50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 52:   SNESCreate(PETSC_COMM_WORLD,&snes);

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:      Create vector data structures; set function evaluation routine
 56:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 58:   VecCreate(PETSC_COMM_SELF,&x);
 59:   VecSetSizes(x,PETSC_DECIDE,n);
 60:   VecSetFromOptions(x);
 61:   VecDuplicate(x,&r);
 62:   VecDuplicate(x,&F);

 64:   SNESSetFunction(snes,r,FormFunction,(void*)F);

 66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:      Create matrix data structures; set Jacobian evaluation routine
 68:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 70:   MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&J);

 72:   /*
 73:      Note that in this case we create separate matrices for the Jacobian
 74:      and preconditioner matrix.  Both of these are computed in the
 75:      routine FormJacobian()
 76:   */
 77:   /*  SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0); */
 78:   SNESSetJacobian(snes,J,J,FormJacobian,0);
 79:   /*  SNESSetJacobian(snes,J,JPrec,FormJacobian,0); */

 81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 82:      Customize nonlinear solver; set runtime options
 83:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 85:   SNESSetFromOptions(snes);

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:      Initialize application:
 89:      Store right-hand-side of PDE and exact solution
 90:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 92:   /* set right hand side and initial guess to be exact solution of continuim problem */
 93: #define SQR(x) ((x)*(x))
 94:   xp = 0.0;
 95:   for (i=0; i<n; i++)
 96:   {
 97:     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.);
 98:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
 99:     v2   = a*xp + (b-a)*PetscPowScalar(xp,k);
100:     VecSetValues(x,1,&i,&v2,INSERT_VALUES);
101:     xp  += h;
102:   }

104:   /* perturb initial guess */
105:   VecGetArray(x,&xx);
106:   for (i=0; i<n; i++) {
107:     v2   = xx[i]*sperturb;
108:     VecSetValues(x,1,&i,&v2,INSERT_VALUES);
109:   }
110:   VecRestoreArray(x,&xx);


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

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

122:   VecDestroy(&x);     VecDestroy(&r);
123:   VecDestroy(&F);     MatDestroy(&J);
124:   SNESDestroy(&snes);
125:   PetscFinalize();

127:   return 0;
128: }

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

137:   VecGetArrayRead(x,&xx);
138:   VecGetArray(f,&ff);
139:   VecGetArray((Vec)dummy,&FF);
140:   VecGetSize(x,&n);
141:   d    = (PetscReal)(n - 1); d2 = d*d;

143:   if (second_order) ff[0] = d*(0.5*d*(-xx[2] + 4.*xx[1] - 3.*xx[0]) - X0DOT);
144:   else ff[0] = d*(d*(xx[1] - xx[0]) - X0DOT);

146:   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];

148:   ff[n-1] = d*d*(xx[n-1] - X1);
149:   VecRestoreArrayRead(x,&xx);
150:   VecRestoreArray(f,&ff);
151:   VecRestoreArray((Vec)dummy,&FF);
152:   return 0;
153: }

155: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat prejac,void *dummy)
156: {
157:   const PetscScalar *xx;
158:   PetscScalar       A[3],d,d2;
159:   PetscInt          i,n,j[3];
160:   PetscErrorCode    ierr;

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

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

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

186:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
187:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
188:   MatAssemblyBegin(prejac,MAT_FINAL_ASSEMBLY);
189:   MatAssemblyEnd(prejac,MAT_FINAL_ASSEMBLY);

191:   VecRestoreArrayRead(x,&xx);
192:   return 0;
193: }