Actual source code: ex1.c

      2: static char help[] = "Newton's method for a two-variable system, sequential.\n\n";

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

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

This examples solves either
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{x^2_0 + x_0 x_1 - 3}{x_0 x_1 + x^2_1 - 6}
\end{equation}
or if the {\tt -hard} options is given
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
\end{equation}

     29: #include <petscsnes.h>

     31: /*
     32:    User-defined routines
     33: */
     34: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
     35: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
     36: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
     37: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);

     39: int main(int argc,char **argv)
     40: {
     41:   SNES           snes;         /* nonlinear solver context */
     42:   KSP            ksp;          /* linear solver context */
     43:   PC             pc;           /* preconditioner context */
     44:   Vec            x,r;          /* solution, residual vectors */
     45:   Mat            J;            /* Jacobian matrix */
     47:   PetscMPIInt    size;
     48:   PetscScalar    pfive = .5,*xx;
     49:   PetscBool      flg;

     51:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     52:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     53:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");

     55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     56:      Create nonlinear solver context
     57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     58:   SNESCreate(PETSC_COMM_WORLD,&snes);

     60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     61:      Create matrix and vector data structures; set corresponding routines
     62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     63:   /*
     64:      Create vectors for solution and nonlinear function
     65:   */
     66:   VecCreate(PETSC_COMM_WORLD,&x);
     67:   VecSetSizes(x,PETSC_DECIDE,2);
     68:   VecSetFromOptions(x);
     69:   VecDuplicate(x,&r);

     71:   /*
     72:      Create Jacobian matrix data structure
     73:   */
     74:   MatCreate(PETSC_COMM_WORLD,&J);
     75:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
     76:   MatSetFromOptions(J);
     77:   MatSetUp(J);

     79:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
     80:   if (!flg) {
     81:     /*
     82:      Set function evaluation routine and vector.
     83:     */
     84:     SNESSetFunction(snes,r,FormFunction1,NULL);

     86:     /*
     87:      Set Jacobian matrix data structure and Jacobian evaluation routine
     88:     */
     89:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
     90:   } else {
     91:     SNESSetFunction(snes,r,FormFunction2,NULL);
     92:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
     93:   }

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:      Customize nonlinear solver; set runtime options
     97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     98:   /*
     99:      Set linear solver defaults for this problem. By extracting the
    100:      KSP and PC contexts from the SNES context, we can then
    101:      directly call any KSP and PC routines to set various options.
    102:   */
    103:   SNESGetKSP(snes,&ksp);
    104:   KSPGetPC(ksp,&pc);
    105:   PCSetType(pc,PCNONE);
    106:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

    108:   /*
    109:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    110:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    111:      These options will override those specified above as long as
    112:      SNESSetFromOptions() is called _after_ any other customization
    113:      routines.
    114:   */
    115:   SNESSetFromOptions(snes);

    117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    118:      Evaluate initial guess; then solve nonlinear system
    119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    120:   if (!flg) {
    121:     VecSet(x,pfive);
    122:   } else {
    123:     VecGetArray(x,&xx);
    124:     xx[0] = 2.0; xx[1] = 3.0;
    125:     VecRestoreArray(x,&xx);
    126:   }
    127:   /*
    128:      Note: The user should initialize the vector, x, with the initial guess
    129:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    130:      to employ an initial guess of zero, the user should explicitly set
    131:      this vector to zero by calling VecSet().
    132:   */

    134:   SNESSolve(snes,NULL,x);
    135:   if (flg) {
    136:     Vec f;
    137:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    138:     SNESGetFunction(snes,&f,0,0);
    139:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
    140:   }

    143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    144:      Free work space.  All PETSc objects should be destroyed when they
    145:      are no longer needed.
    146:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    148:   VecDestroy(&x); VecDestroy(&r);
    149:   MatDestroy(&J); SNESDestroy(&snes);
    150:   PetscFinalize();
    151:   return ierr;
    152: }
    153: /* ------------------------------------------------------------------- */
    154: /*
    155:    FormFunction1 - Evaluates nonlinear function, F(x).

    157:    Input Parameters:
    158: .  snes - the SNES context
    159: .  x    - input vector
    160: .  ctx  - optional user-defined context

    162:    Output Parameter:
    163: .  f - function vector
    164:  */
    165: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
    166: {
    167:   PetscErrorCode    ierr;
    168:   const PetscScalar *xx;
    169:   PetscScalar       *ff;

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

    181:   /* Compute function */
    182:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
    183:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

    185:   /* Restore vectors */
    186:   VecRestoreArrayRead(x,&xx);
    187:   VecRestoreArray(f,&ff);
    188:   return 0;
    189: }
    190: /* ------------------------------------------------------------------- */
    191: /*
    192:    FormJacobian1 - Evaluates Jacobian matrix.

    194:    Input Parameters:
    195: .  snes - the SNES context
    196: .  x - input vector
    197: .  dummy - optional user-defined context (not used here)

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

    211:   /*
    212:      Get pointer to vector data
    213:   */
    214:   VecGetArrayRead(x,&xx);

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

    225:   /*
    226:      Restore vector
    227:   */
    228:   VecRestoreArrayRead(x,&xx);

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

    242: /* ------------------------------------------------------------------- */
    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: /* ------------------------------------------------------------------- */
    273: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
    274: {
    275:   const PetscScalar *xx;
    276:   PetscScalar       A[4];
    277:   PetscErrorCode    ierr;
    278:   PetscInt          idx[2] = {0,1};

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

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

    294:   /*
    295:      Restore vector
    296:   */
    297:   VecRestoreArrayRead(x,&xx);

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

    314: /*TEST

    316:    test:
    317:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
    318:       requires: !single

    320:    test:
    321:       suffix: 2
    322:       requires: !single
    323:       args:  -snes_monitor_short
    324:       output_file: output/ex1_1.out

    326:    test:
    327:       suffix: 3
    328:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab  -snes_monitor_short
    329:       requires: !single
    330:       output_file: output/ex1_1.out

    332:    test:
    333:       suffix: 4
    334:       args: -ksp_view_solution ascii:ex1_2_sol.tmp::append  -snes_monitor_short
    335:       requires: !single
    336:       output_file: output/ex1_1.out

    338:    test:
    339:       suffix: 5
    340:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append  -snes_monitor_short
    341:       requires: !single
    342:       output_file: output/ex1_1.out

    344:    test:
    345:       suffix: 6
    346:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append  -snes_monitor_short
    347:       requires: !single
    348:       output_file: output/ex1_1.out

    350:    test:
    351:       suffix: X
    352:       args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4
    353:       requires: !single x

    355: TEST*/