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:   PetscInt          it,n = 11,i;
 33:   PetscReal         h,xp = 0.0;
 34:   PetscScalar       v;
 35:   const PetscScalar a = X0DOT;
 36:   const PetscScalar b = X1;
 37:   const PetscScalar k = KPOW;
 38:   PetscScalar       v2;
 39:   PetscScalar       *xx;

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

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 47:      Create nonlinear solver context
 48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 50:   SNESCreate(PETSC_COMM_WORLD,&snes);

 52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53:      Create vector data structures; set function evaluation routine
 54:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

 64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:      Create matrix data structures; set Jacobian evaluation routine
 66:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

 79:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80:      Customize nonlinear solver; set runtime options
 81:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 83:   SNESSetFromOptions(snes);

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

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

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

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

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

119:   VecDestroy(&x));     PetscCall(VecDestroy(&r);
120:   VecDestroy(&F));     PetscCall(MatDestroy(&J);
121:   SNESDestroy(&snes);
122:   PetscFinalize();
123:   return 0;
124: }

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

132:   VecGetArrayRead(x,&xx);
133:   VecGetArray(f,&ff);
134:   VecGetArray((Vec)dummy,&FF);
135:   VecGetSize(x,&n);
136:   d    = (PetscReal)(n - 1); d2 = d*d;

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

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

143:   ff[n-1] = d*d*(xx[n-1] - X1);
144:   VecRestoreArrayRead(x,&xx);
145:   VecRestoreArray(f,&ff);
146:   VecRestoreArray((Vec)dummy,&FF);
147:   return 0;
148: }

150: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat prejac,void *dummy)
151: {
152:   const PetscScalar *xx;
153:   PetscScalar       A[3],d,d2;
154:   PetscInt          i,n,j[3];

156:   VecGetSize(x,&n);
157:   VecGetArrayRead(x,&xx);
158:   d    = (PetscReal)(n - 1); d2 = d*d;

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

176:   i    = n-1;
177:   A[0] = d*d;
178:   MatSetValues(prejac,1,&i,1,&i,&A[0],INSERT_VALUES);

180:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
181:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
182:   MatAssemblyBegin(prejac,MAT_FINAL_ASSEMBLY);
183:   MatAssemblyEnd(prejac,MAT_FINAL_ASSEMBLY);

185:   VecRestoreArrayRead(x,&xx);
186:   return 0;
187: }

189: /*TEST

191:    test:
192:       args: -n 14 -snes_monitor_short -snes_converged_reason
193:       requires: !single

195:    test:
196:       suffix: 2
197:       args: -n 15 -snes_monitor_short -snes_converged_reason
198:       requires: !single

200:    test:
201:       suffix: 3
202:       args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
203:       requires: !single

205: TEST*/