Actual source code: ex2.c

      2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
      3: This example employs a user-defined monitoring routine.\n\n";

      5: /*T
      6:    Concepts: SNES^basic uniprocessor example
      7:    Concepts: SNES^setting a user-defined monitoring routine
      8:    Processors: 1
      9: T*/

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

     22:  #include petscsnes.h

     24: /*
     25:    User-defined routines
     26: */

     32: /*
     33:    User-defined context for monitoring
     34: */
     35: typedef struct {
     36:    PetscViewer viewer;
     37: } MonitorCtx;

     41: int main(int argc,char **argv)
     42: {
     43:   SNES           snes;                   /* SNES context */
     44:   Vec            x,r,F,U;             /* vectors */
     45:   Mat            J;                      /* Jacobian matrix */
     46:   MonitorCtx     monP;                   /* monitoring context */
     48:   PetscInt       its,n = 5,i,maxit,maxf;
     49:   PetscMPIInt    size;
     50:   PetscScalar    h,xp,v,none = -1.0;
     51:   PetscReal      abstol,rtol,stol,norm;

     53:   PetscInitialize(&argc,&argv,(char *)0,help);
     54:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     55:   if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
     56:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
     57:   h = 1.0/(n-1);

     59:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     60:      Create nonlinear solver context
     61:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     63:   SNESCreate(PETSC_COMM_WORLD,&snes);

     65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     66:      Create vector data structures; set function evaluation routine
     67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     68:   /*
     69:      Note that we form 1 vector from scratch and then duplicate as needed.
     70:   */
     71:   VecCreate(PETSC_COMM_WORLD,&x);
     72:   VecSetSizes(x,PETSC_DECIDE,n);
     73:   VecSetFromOptions(x);
     74:   VecDuplicate(x,&r);
     75:   VecDuplicate(x,&F);
     76:   VecDuplicate(x,&U);

     78:   /*
     79:      Set function evaluation routine and vector
     80:   */
     81:   SNESSetFunction(snes,r,FormFunction,(void*)F);

     84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     85:      Create matrix data structure; set Jacobian evaluation routine
     86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     88:   MatCreate(PETSC_COMM_WORLD,&J);
     89:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
     90:   MatSetFromOptions(J);

     92:   /*
     93:      Set Jacobian matrix data structure and default Jacobian evaluation
     94:      routine. User can override with:
     95:      -snes_fd : default finite differencing approximation of Jacobian
     96:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
     97:                 (unless user explicitly sets preconditioner)
     98:      -snes_mf_operator : form preconditioning matrix as set by the user,
     99:                          but use matrix-free approx for Jacobian-vector
    100:                          products within Newton-Krylov method
    101:   */

    103:   SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL);

    105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    106:      Customize nonlinear solver; set runtime options
    107:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    109:   /*
    110:      Set an optional user-defined monitoring routine
    111:   */
    112:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
    113:   SNESSetMonitor(snes,Monitor,&monP,0);

    115:   /*
    116:      Set names for some vectors to facilitate monitoring (optional)
    117:   */
    118:   PetscObjectSetName((PetscObject)x,"Approximate Solution");
    119:   PetscObjectSetName((PetscObject)U,"Exact Solution");

    121:   /*
    122:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    123:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    124:   */
    125:   SNESSetFromOptions(snes);

    127:   /*
    128:      Print parameters used for convergence testing (optional) ... just
    129:      to demonstrate this routine; this information is also printed with
    130:      the option -snes_view
    131:   */
    132:   SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
    133:   PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",abstol,rtol,stol,maxit,maxf);

    135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    136:      Initialize application:
    137:      Store right-hand-side of PDE and exact solution
    138:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    140:   xp = 0.0;
    141:   for (i=0; i<n; i++) {
    142:     v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
    143:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
    144:     v= xp*xp*xp;
    145:     VecSetValues(U,1,&i,&v,INSERT_VALUES);
    146:     xp += h;
    147:   }

    149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    150:      Evaluate initial guess; then solve nonlinear system
    151:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    152:   /*
    153:      Note: The user should initialize the vector, x, with the initial guess
    154:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    155:      to employ an initial guess of zero, the user should explicitly set
    156:      this vector to zero by calling VecSet().
    157:   */
    158:   FormInitialGuess(x);
    159:   SNESSolve(snes,PETSC_NULL,x);
    160:   SNESGetIterationNumber(snes,&its);
    161:   PetscPrintf(PETSC_COMM_WORLD,"number of Newton iterations = %D\n\n",its);

    163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    164:      Check solution and clean up
    165:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    167:   /*
    168:      Check the error
    169:   */
    170:   VecAXPY(x,none,U);
    171:   VecNorm(x,NORM_2,&norm);
    172:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);

    175:   /*
    176:      Free work space.  All PETSc objects should be destroyed when they
    177:      are no longer needed.
    178:   */
    179:   VecDestroy(x);  VecDestroy(r);
    180:   VecDestroy(U);  VecDestroy(F);
    181:   MatDestroy(J);  SNESDestroy(snes);
    182:   PetscViewerDestroy(monP.viewer);
    183:   PetscFinalize();

    185:   return 0;
    186: }
    187: /* ------------------------------------------------------------------- */
    190: /*
    191:    FormInitialGuess - Computes initial guess.

    193:    Input/Output Parameter:
    194: .  x - the solution vector
    195: */
    196: PetscErrorCode FormInitialGuess(Vec x)
    197: {
    199:    PetscScalar    pfive = .50;
    200:    VecSet(x,pfive);
    201:    return 0;
    202: }
    203: /* ------------------------------------------------------------------- */
    206: /*
    207:    FormFunction - Evaluates nonlinear function, F(x).

    209:    Input Parameters:
    210: .  snes - the SNES context
    211: .  x - input vector
    212: .  ctx - optional user-defined context, as set by SNESSetFunction()

    214:    Output Parameter:
    215: .  f - function vector

    217:    Note:
    218:    The user-defined context can contain any application-specific data
    219:    needed for the function evaluation (such as various parameters, work
    220:    vectors, and grid information).  In this program the context is just
    221:    a vector containing the right-hand-side of the discretized PDE.
    222:  */

    224: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
    225: {
    226:    Vec            g = (Vec)ctx;
    227:    PetscScalar    *xx,*ff,*gg,d;
    229:    PetscInt       i,n;

    231:   /*
    232:      Get pointers to vector data.
    233:        - For default PETSc vectors, VecGetArray() returns a pointer to
    234:          the data array.  Otherwise, the routine is implementation dependent.
    235:        - You MUST call VecRestoreArray() when you no longer need access to
    236:          the array.
    237:   */
    238:    VecGetArray(x,&xx);
    239:    VecGetArray(f,&ff);
    240:    VecGetArray(g,&gg);

    242:   /*
    243:      Compute function
    244:   */
    245:    VecGetSize(x,&n);
    246:    d = (PetscReal)(n - 1); d = d*d;
    247:    ff[0]   = xx[0];
    248:    for (i=1; i<n-1; i++) {
    249:      ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
    250:    }
    251:    ff[n-1] = xx[n-1] - 1.0;

    253:   /*
    254:      Restore vectors
    255:   */
    256:   VecRestoreArray(x,&xx);
    257:   VecRestoreArray(f,&ff);
    258:   VecRestoreArray(g,&gg);
    259:   return 0;
    260: }
    261: /* ------------------------------------------------------------------- */
    264: /*
    265:    FormJacobian - Evaluates Jacobian matrix.

    267:    Input Parameters:
    268: .  snes - the SNES context
    269: .  x - input vector
    270: .  dummy - optional user-defined context (not used here)

    272:    Output Parameters:
    273: .  jac - Jacobian matrix
    274: .  B - optionally different preconditioning matrix
    275: .  flag - flag indicating matrix structure
    276: */

    278: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)
    279: {
    280:   PetscScalar    *xx,A[3],d;
    282:   PetscInt       i,n,j[3];

    284:   /*
    285:      Get pointer to vector data
    286:   */
    287:   VecGetArray(x,&xx);

    289:   /*
    290:      Compute Jacobian entries and insert into matrix.
    291:       - Note that in this case we set all elements for a particular
    292:         row at once.
    293:   */
    294:   VecGetSize(x,&n);
    295:   d = (PetscReal)(n - 1); d = d*d;

    297:   /*
    298:      Interior grid points
    299:   */
    300:   for (i=1; i<n-1; i++) {
    301:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
    302:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
    303:     MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
    304:   }

    306:   /*
    307:      Boundary points
    308:   */
    309:   i = 0;   A[0] = 1.0;
    310:   MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
    311:   i = n-1; A[0] = 1.0;
    312:   MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);

    314:   /*
    315:      Restore vector
    316:   */
    317:   VecRestoreArray(x,&xx);

    319:   /*
    320:      Assemble matrix
    321:   */
    322:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    323:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

    325:   return 0;
    326: }
    327: /* ------------------------------------------------------------------- */
    330: /*
    331:    Monitor - User-defined monitoring routine that views the
    332:    current iterate with an x-window plot.

    334:    Input Parameters:
    335:    snes - the SNES context
    336:    its - iteration number
    337:    norm - 2-norm function value (may be estimated)
    338:    ctx - optional user-defined context for private data for the
    339:          monitor routine, as set by SNESSetMonitor()

    341:    Note:
    342:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
    343:    such as -nox to deactivate all x-window output.
    344:  */
    345: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
    346: {
    348:   MonitorCtx     *monP = (MonitorCtx*) ctx;
    349:   Vec            x;

    351:   PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,fnorm);
    352:   SNESGetSolution(snes,&x);
    353:   VecView(x,monP->viewer);
    354:   return 0;
    355: }