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:   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);if (ierr) return ierr;
     44:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
     45:   PetscOptionsGetBool(NULL,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();
    126:   return ierr;
    127: }

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

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

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

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

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

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

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

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

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

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

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

    195: /*TEST

    197:    test:
    198:       args: -n 14 -snes_monitor_short -snes_converged_reason
    199:       requires: !single

    201:    test:
    202:       suffix: 2
    203:       args: -n 15 -snes_monitor_short -snes_converged_reason
    204:       requires: !single

    206:    test:
    207:       suffix: 3
    208:       args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
    209:       requires: !single

    211: TEST*/