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*/

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

     24: #include <petscsnes.h>

     26: /*
     27:    User-defined routines
     28: */
     29: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
     30: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
     31: extern PetscErrorCode FormInitialGuess(Vec);
     32: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);

     34: /*
     35:    User-defined context for monitoring
     36: */
     37: typedef struct {
     38:   PetscViewer viewer;
     39: } 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);if (ierr) return ierr;
     54:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     55:   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
     56:   PetscOptionsGetInt(NULL,NULL,"-n",&n,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);
     91:   MatSeqAIJSetPreallocation(J,3,NULL);

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

    104:   SNESSetJacobian(snes,J,J,FormJacobian,NULL);

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

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

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

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

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

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

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

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

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

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

    176:   /*
    177:      Free work space.  All PETSc objects should be destroyed when they
    178:      are no longer needed.
    179:   */
    180:   VecDestroy(&x);  VecDestroy(&r);
    181:   VecDestroy(&U);  VecDestroy(&F);
    182:   MatDestroy(&J);  SNESDestroy(&snes);
    183:   PetscViewerDestroy(&monP.viewer);
    184:   PetscFinalize();
    185:   return ierr;
    186: }
    187: /* ------------------------------------------------------------------- */
    188: /*
    189:    FormInitialGuess - Computes initial guess.

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

    205:    Input Parameters:
    206: .  snes - the SNES context
    207: .  x - input vector
    208: .  ctx - optional user-defined context, as set by SNESSetFunction()

    210:    Output Parameter:
    211: .  f - function vector

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

    220: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
    221: {
    222:   Vec               g = (Vec)ctx;
    223:   const PetscScalar *xx,*gg;
    224:   PetscScalar       *ff,d;
    225:   PetscErrorCode    ierr;
    226:   PetscInt          i,n;

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

    239:   /*
    240:      Compute function
    241:   */
    242:   VecGetSize(x,&n);
    243:   d     = (PetscReal)(n - 1); d = d*d;
    244:   ff[0] = xx[0];
    245:   for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
    246:   ff[n-1] = xx[n-1] - 1.0;

    248:   /*
    249:      Restore vectors
    250:   */
    251:   VecRestoreArrayRead(x,&xx);
    252:   VecRestoreArray(f,&ff);
    253:   VecRestoreArrayRead(g,&gg);
    254:   return 0;
    255: }
    256: /* ------------------------------------------------------------------- */
    257: /*
    258:    FormJacobian - Evaluates Jacobian matrix.

    260:    Input Parameters:
    261: .  snes - the SNES context
    262: .  x - input vector
    263: .  dummy - optional user-defined context (not used here)

    265:    Output Parameters:
    266: .  jac - Jacobian matrix
    267: .  B - optionally different preconditioning matrix

    269: */

    271: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
    272: {
    273:   const PetscScalar *xx;
    274:   PetscScalar       A[3],d;
    275:   PetscErrorCode    ierr;
    276:   PetscInt          i,n,j[3];

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

    283:   /*
    284:      Compute Jacobian entries and insert into matrix.
    285:       - Note that in this case we set all elements for a particular
    286:         row at once.
    287:   */
    288:   VecGetSize(x,&n);
    289:   d    = (PetscReal)(n - 1); d = d*d;

    291:   /*
    292:      Interior grid points
    293:   */
    294:   for (i=1; i<n-1; i++) {
    295:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
    296:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
    297:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
    298:   }

    300:   /*
    301:      Boundary points
    302:   */
    303:   i = 0;   A[0] = 1.0;

    305:   MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);

    307:   i = n-1; A[0] = 1.0;

    309:   MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);

    311:   /*
    312:      Restore vector
    313:   */
    314:   VecRestoreArrayRead(x,&xx);

    316:   /*
    317:      Assemble matrix
    318:   */
    319:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    320:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    321:   if (jac != B) {
    322:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    323:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
    324:   }
    325:   return 0;
    326: }
    327: /* ------------------------------------------------------------------- */
    328: /*
    329:    Monitor - User-defined monitoring routine that views the
    330:    current iterate with an x-window plot.

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

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

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

    356: /*TEST

    358:    test:
    359:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

    361:    test:
    362:       suffix: 2
    363:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
    364:       requires: !single

    366:    test:
    367:       suffix: 3
    368:       args: -nox -malloc no -options_left no -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

    370: TEST*/