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

      2: static char help[] = "Newton's method for a two-variable system, sequential.\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: /*
     20:    User-defined routines
     21: */
     22: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     23: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
     24: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     25: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);

     29: int main(int argc,char **argv)
     30: {
     31:   SNES           snes;         /* nonlinear solver context */
     32:   KSP            ksp;          /* linear solver context */
     33:   PC             pc;           /* preconditioner context */
     34:   Vec            x,r;          /* solution, residual vectors */
     35:   Mat            J;            /* Jacobian matrix */
     37:   PetscInt       its;
     38:   PetscMPIInt    size;
     39:   PetscScalar    pfive = .5,*xx;
     40:   PetscBool      flg;

     42:   PetscInitialize(&argc,&argv,(char*)0,help);
     43:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     44:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");

     46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     47:      Create nonlinear solver context
     48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     49:   SNESCreate(PETSC_COMM_WORLD,&snes);

     51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     52:      Create matrix and vector data structures; set corresponding routines
     53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     54:   /*
     55:      Create vectors for solution and nonlinear function
     56:   */
     57:   VecCreate(PETSC_COMM_WORLD,&x);
     58:   VecSetSizes(x,PETSC_DECIDE,2);
     59:   VecSetFromOptions(x);
     60:   VecDuplicate(x,&r);

     62:   /*
     63:      Create Jacobian matrix data structure
     64:   */
     65:   MatCreate(PETSC_COMM_WORLD,&J);
     66:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
     67:   MatSetFromOptions(J);
     68:   MatSetUp(J);

     70:   PetscOptionsHasName(NULL,"-hard",&flg);
     71:   if (!flg) {
     72:     /*
     73:      Set function evaluation routine and vector.
     74:     */
     75:     SNESSetFunction(snes,r,FormFunction1,NULL);

     77:     /*
     78:      Set Jacobian matrix data structure and Jacobian evaluation routine
     79:     */
     80:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
     81:   } else {
     82:     SNESSetFunction(snes,r,FormFunction2,NULL);
     83:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
     84:   }

     86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     87:      Customize nonlinear solver; set runtime options
     88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     89:   /*
     90:      Set linear solver defaults for this problem. By extracting the
     91:      KSP and PC contexts from the SNES context, we can then
     92:      directly call any KSP and PC routines to set various options.
     93:   */
     94:   SNESGetKSP(snes,&ksp);
     95:   KSPGetPC(ksp,&pc);
     96:   PCSetType(pc,PCNONE);
     97:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

     99:   /*
    100:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    101:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    102:      These options will override those specified above as long as
    103:      SNESSetFromOptions() is called _after_ any other customization
    104:      routines.
    105:   */
    106:   SNESSetFromOptions(snes);

    108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    109:      Evaluate initial guess; then solve nonlinear system
    110:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    111:   if (!flg) {
    112:     VecSet(x,pfive);
    113:   } else {
    114:     VecGetArray(x,&xx);
    115:     xx[0] = 2.0; xx[1] = 3.0;
    116:     VecRestoreArray(x,&xx);
    117:   }
    118:   /*
    119:      Note: The user should initialize the vector, x, with the initial guess
    120:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    121:      to employ an initial guess of zero, the user should explicitly set
    122:      this vector to zero by calling VecSet().
    123:   */

    125:   SNESSolve(snes,NULL,x);
    126:   SNESGetIterationNumber(snes,&its);
    127:   if (flg) {
    128:     Vec f;
    129:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    130:     SNESGetFunction(snes,&f,0,0);
    131:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
    132:   }

    134:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);

    136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    137:      Free work space.  All PETSc objects should be destroyed when they
    138:      are no longer needed.
    139:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    141:   VecDestroy(&x); VecDestroy(&r);
    142:   MatDestroy(&J); SNESDestroy(&snes);
    143:   PetscFinalize();
    144:   return 0;
    145: }
    146: /* ------------------------------------------------------------------- */
    149: /*
    150:    FormFunction1 - Evaluates nonlinear function, F(x).

    152:    Input Parameters:
    153: .  snes - the SNES context
    154: .  x    - input vector
    155: .  ctx  - optional user-defined context

    157:    Output Parameter:
    158: .  f - function vector
    159:  */
    160: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
    161: {
    162:   PetscErrorCode    ierr;
    163:   const PetscScalar *xx;
    164:   PetscScalar       *ff;

    166:   /*
    167:    Get pointers to vector data.
    168:       - For default PETSc vectors, VecGetArray() returns a pointer to
    169:         the data array.  Otherwise, the routine is implementation dependent.
    170:       - You MUST call VecRestoreArray() when you no longer need access to
    171:         the array.
    172:    */
    173:   VecGetArrayRead(x,&xx);
    174:   VecGetArray(f,&ff);

    176:   /* Compute function */
    177:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
    178:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

    180:   /* Restore vectors */
    181:   VecRestoreArrayRead(x,&xx);
    182:   VecRestoreArray(f,&ff);
    183:   return 0;
    184: }
    185: /* ------------------------------------------------------------------- */
    188: /*
    189:    FormJacobian1 - Evaluates Jacobian matrix.

    191:    Input Parameters:
    192: .  snes - the SNES context
    193: .  x - input vector
    194: .  dummy - optional user-defined context (not used here)

    196:    Output Parameters:
    197: .  jac - Jacobian matrix
    198: .  B - optionally different preconditioning matrix
    199: .  flag - flag indicating matrix structure
    200: */
    201: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    202: {
    203:   const PetscScalar *xx;
    204:   PetscScalar       A[4];
    205:   PetscErrorCode    ierr;
    206:   PetscInt          idx[2] = {0,1};

    208:   /*
    209:      Get pointer to vector data
    210:   */
    211:   VecGetArrayRead(x,&xx);

    213:   /*
    214:      Compute Jacobian entries and insert into matrix.
    215:       - Since this is such a small problem, we set all entries for
    216:         the matrix at once.
    217:   */
    218:   A[0]  = 2.0*xx[0] + xx[1]; A[1] = xx[0];
    219:   A[2]  = xx[1]; A[3] = xx[0] + 2.0*xx[1];
    220:   MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
    221:   *flag = SAME_NONZERO_PATTERN;

    223:   /*
    224:      Restore vector
    225:   */
    226:   VecRestoreArrayRead(x,&xx);

    228:   /*
    229:      Assemble matrix
    230:   */
    231:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
    232:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
    233:   if (*jac != *B) {
    234:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    235:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    236:   }
    237:   return 0;
    238: }

    240: /* ------------------------------------------------------------------- */
    243: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
    244: {
    245:   PetscErrorCode    ierr;
    246:   const PetscScalar *xx;
    247:   PetscScalar       *ff;

    249:   /*
    250:      Get pointers to vector data.
    251:        - For default PETSc vectors, VecGetArray() returns a pointer to
    252:          the data array.  Otherwise, the routine is implementation dependent.
    253:        - You MUST call VecRestoreArray() when you no longer need access to
    254:          the array.
    255:   */
    256:   VecGetArrayRead(x,&xx);
    257:   VecGetArray(f,&ff);

    259:   /*
    260:      Compute function
    261:   */
    262:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
    263:   ff[1] = xx[1];

    265:   /*
    266:      Restore vectors
    267:   */
    268:   VecRestoreArrayRead(x,&xx);
    269:   VecRestoreArray(f,&ff);
    270:   return 0;
    271: }
    272: /* ------------------------------------------------------------------- */
    275: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    276: {
    277:   const PetscScalar *xx;
    278:   PetscScalar       A[4];
    279:   PetscErrorCode    ierr;
    280:   PetscInt          idx[2] = {0,1};

    282:   /*
    283:      Get pointer to vector data
    284:   */
    285:   VecGetArrayRead(x,&xx);

    287:   /*
    288:      Compute Jacobian entries and insert into matrix.
    289:       - Since this is such a small problem, we set all entries for
    290:         the matrix at once.
    291:   */
    292:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
    293:   A[2]  = 0.0;                     A[3] = 1.0;
    294:   MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
    295:   *flag = SAME_NONZERO_PATTERN;

    297:   /*
    298:      Restore vector
    299:   */
    300:   VecRestoreArrayRead(x,&xx);

    302:   /*
    303:      Assemble matrix
    304:   */
    305:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
    306:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
    307:   if (*jac != *B) {
    308:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    309:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    310:   }
    311:   return 0;
    312: }