Actual source code: ex42.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n";

      4: /*T
      5:    Concepts: SNES^basic example
      6: T*/

      8: /*
      9:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
     10:    file automatically includes:
     11:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     12:      petscmat.h - matrices
     13:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     14:      petscviewer.h - viewers               petscpc.h  - preconditioners
     15:      petscksp.h   - linear solvers
     16: */
     17: #include <petscsnes.h>

     19: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     20: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);

     24: int main(int argc,char **argv)
     25: {
     26:   SNES                snes;    /* nonlinear solver context */
     27:   Vec                 x,r;     /* solution, residual vectors */
     28:   Mat                 J;       /* Jacobian matrix */
     29:   PetscErrorCode      ierr;
     30:   PetscInt            its;
     31:   PetscScalar         *xx;
     32:   SNESConvergedReason reason;

     34:   PetscInitialize(&argc,&argv,(char*)0,help);

     36:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     37:      Create nonlinear solver context
     38:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     39:   SNESCreate(PETSC_COMM_WORLD,&snes);

     41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     42:      Create matrix and vector data structures; set corresponding routines
     43:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     44:   /*
     45:      Create vectors for solution and nonlinear function
     46:   */
     47:   VecCreate(PETSC_COMM_WORLD,&x);
     48:   VecSetSizes(x,PETSC_DECIDE,2);
     49:   VecSetFromOptions(x);
     50:   VecDuplicate(x,&r);

     52:   /*
     53:      Create Jacobian matrix data structure
     54:   */
     55:   MatCreate(PETSC_COMM_WORLD,&J);
     56:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
     57:   MatSetFromOptions(J);
     58:   MatSetUp(J);

     60:   /*
     61:      Set function evaluation routine and vector.
     62:   */
     63:   SNESSetFunction(snes,r,FormFunction1,NULL);

     65:   /*
     66:      Set Jacobian matrix data structure and Jacobian evaluation routine
     67:   */
     68:   SNESSetJacobian(snes,J,J,FormJacobian1,NULL);

     70:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     71:      Customize nonlinear solver; set runtime options
     72:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     73:   SNESSetFromOptions(snes);

     75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76:      Evaluate initial guess; then solve nonlinear system
     77:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     78:   VecGetArray(x,&xx);
     79:   xx[0] = -1.2; xx[1] = 1.0;
     80:   VecRestoreArray(x,&xx);

     82:   /*
     83:      Note: The user should initialize the vector, x, with the initial guess
     84:      for the nonlinear solver prior to calling SNESSolve().  In particular,
     85:      to employ an initial guess of zero, the user should explicitly set
     86:      this vector to zero by calling VecSet().
     87:   */

     89:   SNESSolve(snes,NULL,x);
     90:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
     91:   SNESGetIterationNumber(snes,&its);
     92:   SNESGetConvergedReason(snes,&reason);
     93:   /*
     94:      Some systems computes a residual that is identically zero, thus converging
     95:      due to FNORM_ABS, others converge due to FNORM_RELATIVE.  Here, we only
     96:      report the reason if the iteration did not converge so that the tests are
     97:      reproducible.
     98:   */
     99:   PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],its);

    101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    102:      Free work space.  All PETSc objects should be destroyed when they
    103:      are no longer needed.
    104:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    106:   VecDestroy(&x); VecDestroy(&r);
    107:   MatDestroy(&J); SNESDestroy(&snes);
    108:   PetscFinalize();
    109:   return 0;
    110: }
    111: /* ------------------------------------------------------------------- */
    114: /*
    115:    FormFunction1 - Evaluates nonlinear function, F(x).

    117:    Input Parameters:
    118: .  snes - the SNES context
    119: .  x    - input vector
    120: .  ctx  - optional user-defined context

    122:    Output Parameter:
    123: .  f - function vector
    124:  */
    125: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
    126: {
    128:   PetscScalar    *xx,*ff;

    130:   /*
    131:     Get pointers to vector data.
    132:     - For default PETSc vectors, VecGetArray() returns a pointer to
    133:     the data array.  Otherwise, the routine is implementation dependent.
    134:     - You MUST call VecRestoreArray() when you no longer need access to
    135:     the array.
    136:   */
    137:   VecGetArray(x,&xx);
    138:   VecGetArray(f,&ff);

    140:   /* Compute function */
    141:   ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
    142:   ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];

    144:   /* Restore vectors */
    145:   VecRestoreArray(x,&xx);
    146:   VecRestoreArray(f,&ff);
    147:   return 0;
    148: }
    149: /* ------------------------------------------------------------------- */
    152: /*
    153:    FormJacobian1 - Evaluates Jacobian matrix.

    155:    Input Parameters:
    156: .  snes - the SNES context
    157: .  x - input vector
    158: .  dummy - optional user-defined context (not used here)

    160:    Output Parameters:
    161: .  jac - Jacobian matrix
    162: .  B - optionally different preconditioning matrix
    163: .  flag - flag indicating matrix structure
    164: */
    165: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    166: {
    167:   PetscScalar    *xx,A[4];
    169:   PetscInt       idx[2] = {0,1};

    171:   /*
    172:      Get pointer to vector data
    173:   */
    174:   VecGetArray(x,&xx);

    176:   /*
    177:      Compute Jacobian entries and insert into matrix.
    178:       - Since this is such a small problem, we set all entries for
    179:         the matrix at once.
    180:   */
    181:   A[0]  = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1];
    182:   A[1]  = -400.0*xx[0];
    183:   A[2]  = -400.0*xx[0];
    184:   A[3]  = 200;
    185:   MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
    186:   *flag = SAME_NONZERO_PATTERN;

    188:   /*
    189:      Restore vector
    190:   */
    191:   VecRestoreArray(x,&xx);

    193:   /*
    194:      Assemble matrix
    195:   */
    196:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
    197:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
    198:   if (*jac != *B) {
    199:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    200:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    201:   }
    202:   return 0;
    203: }