Actual source code: ex59.c

      1: static const char help[] = "Tries to solve u`` + u^{2} = f for an easy case and an impossible case.\n\n";

      3: /*
      4:        This example was contributed by Peter Graf to show how SNES fails when given a nonlinear problem with no solution.

      6:        Run with -n 14 to see fail to converge and -n 15 to see convergence

      8:        The option -second_order uses a different discretization of the Neumann boundary condition and always converges

     10: */

     12: #include <petscsnes.h>

     14: PetscBool second_order = PETSC_FALSE;
     15: #define X0DOT -2.0
     16: #define X1    5.0
     17: #define KPOW  2.0
     18: const PetscScalar sperturb = 1.1;

     20: /*
     21:    User-defined routines
     22: */
     23: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
     24: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);

     26: int main(int argc, char **argv)
     27: {
     28:   SNES              snes;    /* SNES context */
     29:   Vec               x, r, F; /* vectors */
     30:   Mat               J;       /* Jacobian */
     31:   PetscInt          it, n = 11, i;
     32:   PetscReal         h, xp = 0.0;
     33:   PetscScalar       v;
     34:   const PetscScalar a = X0DOT;
     35:   const PetscScalar b = X1;
     36:   const PetscScalar k = KPOW;
     37:   PetscScalar       v2;
     38:   PetscScalar      *xx;

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

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

     50:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

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

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

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

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

     68:   PetscCall(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 matrix used to compute the preconditioner.  Both of these are computed in the
     73:      routine FormJacobian()
     74:   */
     75:   /*  PetscCall(SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0)); */
     76:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, 0));
     77:   /*  PetscCall(SNESSetJacobian(snes,J,JPrec,FormJacobian,0)); */

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

     83:   PetscCall(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 continuum problem */
     91: #define SQR(x) ((x) * (x))
     92:   xp = 0.0;
     93:   for (i = 0; i < n; i++) {
     94:     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.);
     95:     PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
     96:     v2 = a * xp + (b - a) * PetscPowScalar(xp, k);
     97:     PetscCall(VecSetValues(x, 1, &i, &v2, INSERT_VALUES));
     98:     xp += h;
     99:   }

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

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

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

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

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

    133:   PetscFunctionBeginUser;
    134:   PetscCall(VecGetArrayRead(x, &xx));
    135:   PetscCall(VecGetArray(f, &ff));
    136:   PetscCall(VecGetArray((Vec)dummy, &FF));
    137:   PetscCall(VecGetSize(x, &n));
    138:   d  = (PetscReal)(n - 1);
    139:   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:   PetscCall(VecRestoreArrayRead(x, &xx));
    148:   PetscCall(VecRestoreArray(f, &ff));
    149:   PetscCall(VecRestoreArray((Vec)dummy, &FF));
    150:   PetscFunctionReturn(PETSC_SUCCESS);
    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];

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

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

    191:   i    = n - 1;
    192:   A[0] = d * d;
    193:   PetscCall(MatSetValues(prejac, 1, &i, 1, &i, &A[0], INSERT_VALUES));

    195:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    196:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
    197:   PetscCall(MatAssemblyBegin(prejac, MAT_FINAL_ASSEMBLY));
    198:   PetscCall(MatAssemblyEnd(prejac, MAT_FINAL_ASSEMBLY));

    200:   PetscCall(VecRestoreArrayRead(x, &xx));
    201:   PetscFunctionReturn(PETSC_SUCCESS);
    202: }

    204: /*TEST

    206:    test:
    207:       args: -n 14 -snes_monitor_short -snes_converged_reason
    208:       requires: !single

    210:    test:
    211:       suffix: 2
    212:       args: -n 15 -snes_monitor_short -snes_converged_reason
    213:       requires: !single

    215:    test:
    216:       suffix: 3
    217:       args: -n 14 -second_order -snes_monitor_short -snes_converged_reason
    218:       requires: !single

    220: TEST*/