Actual source code: ex6.c

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

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

     11: #include <petscsnes.h>

     13: /*
     14:    User-defined routines
     15: */
     16: PETSC_EXTERN PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
     17: PETSC_EXTERN PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
     18: PETSC_EXTERN PetscErrorCode FormInitialGuess(Vec);
     19: PETSC_EXTERN PetscErrorCode MySNESConvergedReasonView(SNES,void*);
     20: PETSC_EXTERN PetscErrorCode MyKSPConvergedReasonView(KSP,void*);

     22: /*
     23:    User-defined context for monitoring
     24: */
     25: typedef struct {
     26:   PetscViewer viewer;
     27: } ReasonViewCtx;

     29: int main(int argc,char **argv)
     30: {
     31:   SNES           snes;                /* SNES context */
     32:   KSP            ksp;                 /* KSP context */
     33:   Vec            x,r,F,U;             /* vectors */
     34:   Mat            J;                   /* Jacobian matrix */
     35:   ReasonViewCtx     monP;             /* monitoring context */
     37:   PetscInt       its,n = 5,i;
     38:   PetscMPIInt    size;
     39:   PetscScalar    h,xp,v;
     40:   MPI_Comm       comm;

     42:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     43:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     44:   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
     45:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
     46:   h    = 1.0/(n-1);
     47:   comm = PETSC_COMM_WORLD;
     48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     49:      Create nonlinear solver context
     50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     52:   SNESCreate(comm,&snes);

     54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     55:      Create vector data structures; set function evaluation routine
     56:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     57:   /*
     58:      Note that we form 1 vector from scratch and then duplicate as needed.
     59:   */
     60:   VecCreate(comm,&x);
     61:   VecSetSizes(x,PETSC_DECIDE,n);
     62:   VecSetFromOptions(x);
     63:   VecDuplicate(x,&r);
     64:   VecDuplicate(x,&F);
     65:   VecDuplicate(x,&U);

     67:   /*
     68:      Set function evaluation routine and vector
     69:   */
     70:   SNESSetFunction(snes,r,FormFunction,(void*)F);

     73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     74:      Create matrix data structure; set Jacobian evaluation routine
     75:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     77:   MatCreate(comm,&J);
     78:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
     79:   MatSetFromOptions(J);
     80:   MatSeqAIJSetPreallocation(J,3,NULL);

     82:   /*
     83:      Set Jacobian matrix data structure and default Jacobian evaluation
     84:      routine. User can override with:
     85:      -snes_fd : default finite differencing approximation of Jacobian
     86:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
     87:                 (unless user explicitly sets preconditioner)
     88:      -snes_mf_operator : form preconditioning matrix as set by the user,
     89:                          but use matrix-free approx for Jacobian-vector
     90:                          products within Newton-Krylov method
     91:   */

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

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:      Customize nonlinear solver; set runtime options
     97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     99:   /*
    100:      Set an optional user-defined reasonview routine
    101:   */
    102:   PetscViewerASCIIGetStdout(comm,&monP.viewer);
    103:   /* Just make sure we can not repeat addding the same function
    104:    * PETSc will be able to igore the repeated function
    105:    */
    106:   for (i=0; i<4; i++){
    107:     SNESConvergedReasonViewSet(snes,MySNESConvergedReasonView,&monP,0);
    108:   }
    109:   SNESGetKSP(snes,&ksp);
    110:   /* Just make sure we can not repeat addding the same function
    111:    * PETSc will be able to igore the repeated function
    112:    */
    113:   for (i=0; i<4; i++){
    114:     KSPConvergedReasonViewSet(ksp,MyKSPConvergedReasonView,&monP,0);
    115:   }
    116:   /*
    117:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    118:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    119:   */
    120:   SNESSetFromOptions(snes);

    122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    123:      Initialize application:
    124:      Store right-hand-side of PDE and exact solution
    125:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    127:   xp = 0.0;
    128:   for (i=0; i<n; i++) {
    129:     v    = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
    130:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
    131:     v    = xp*xp*xp;
    132:     VecSetValues(U,1,&i,&v,INSERT_VALUES);
    133:     xp  += h;
    134:   }

    136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    137:      Evaluate initial guess; then solve nonlinear system
    138:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    139:   /*
    140:      Note: The user should initialize the vector, x, with the initial guess
    141:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    142:      to employ an initial guess of zero, the user should explicitly set
    143:      this vector to zero by calling VecSet().
    144:   */
    145:   FormInitialGuess(x);
    146:   SNESSolve(snes,NULL,x);
    147:   SNESGetIterationNumber(snes,&its);

    149:   /*
    150:      Free work space.  All PETSc objects should be destroyed when they
    151:      are no longer needed.
    152:   */
    153:   VecDestroy(&x);  VecDestroy(&r);
    154:   VecDestroy(&U);  VecDestroy(&F);
    155:   MatDestroy(&J);  SNESDestroy(&snes);
    156:   /*PetscViewerDestroy(&monP.viewer);*/
    157:   PetscFinalize();
    158:   return ierr;
    159: }
    160: /* ------------------------------------------------------------------- */
    161: /*
    162:    FormInitialGuess - Computes initial guess.

    164:    Input/Output Parameter:
    165: .  x - the solution vector
    166: */
    167: PetscErrorCode FormInitialGuess(Vec x)
    168: {
    170:   PetscScalar    pfive = .50;
    171:   VecSet(x,pfive);
    172:   return 0;
    173: }
    174: /* ------------------------------------------------------------------- */
    175: /*
    176:    FormFunction - Evaluates nonlinear function, F(x).

    178:    Input Parameters:
    179: .  snes - the SNES context
    180: .  x - input vector
    181: .  ctx - optional user-defined context, as set by SNESSetFunction()

    183:    Output Parameter:
    184: .  f - function vector

    186:    Note:
    187:    The user-defined context can contain any application-specific data
    188:    needed for the function evaluation (such as various parameters, work
    189:    vectors, and grid information).  In this program the context is just
    190:    a vector containing the right-hand-side of the discretized PDE.
    191:  */

    193: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
    194: {
    195:   Vec               g = (Vec)ctx;
    196:   const PetscScalar *xx,*gg;
    197:   PetscScalar       *ff,d;
    198:   PetscErrorCode    ierr;
    199:   PetscInt          i,n;

    201:   /*
    202:      Get pointers to vector data.
    203:        - For default PETSc vectors, VecGetArray() returns a pointer to
    204:          the data array.  Otherwise, the routine is implementation dependent.
    205:        - You MUST call VecRestoreArray() when you no longer need access to
    206:          the array.
    207:   */
    208:   VecGetArrayRead(x,&xx);
    209:   VecGetArray(f,&ff);
    210:   VecGetArrayRead(g,&gg);

    212:   /*
    213:      Compute function
    214:   */
    215:   VecGetSize(x,&n);
    216:   d     = (PetscReal)(n - 1); d = d*d;
    217:   ff[0] = xx[0];
    218:   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];
    219:   ff[n-1] = xx[n-1] - 1.0;

    221:   /*
    222:      Restore vectors
    223:   */
    224:   VecRestoreArrayRead(x,&xx);
    225:   VecRestoreArray(f,&ff);
    226:   VecRestoreArrayRead(g,&gg);
    227:   return 0;
    228: }
    229: /* ------------------------------------------------------------------- */
    230: /*
    231:    FormJacobian - Evaluates Jacobian matrix.

    233:    Input Parameters:
    234: .  snes - the SNES context
    235: .  x - input vector
    236: .  dummy - optional user-defined context (not used here)

    238:    Output Parameters:
    239: .  jac - Jacobian matrix
    240: .  B - optionally different preconditioning matrix

    242: */

    244: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
    245: {
    246:   const PetscScalar *xx;
    247:   PetscScalar       A[3],d;
    248:   PetscErrorCode    ierr;
    249:   PetscInt          i,n,j[3];

    251:   /*
    252:      Get pointer to vector data
    253:   */
    254:   VecGetArrayRead(x,&xx);

    256:   /*
    257:      Compute Jacobian entries and insert into matrix.
    258:       - Note that in this case we set all elements for a particular
    259:         row at once.
    260:   */
    261:   VecGetSize(x,&n);
    262:   d    = (PetscReal)(n - 1); d = d*d;

    264:   /*
    265:      Interior grid points
    266:   */
    267:   for (i=1; i<n-1; i++) {
    268:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
    269:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
    270:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
    271:   }

    273:   /*
    274:      Boundary points
    275:   */
    276:   i = 0;   A[0] = 1.0;

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

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

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

    284:   /*
    285:      Restore vector
    286:   */
    287:   VecRestoreArrayRead(x,&xx);

    289:   /*
    290:      Assemble matrix
    291:   */
    292:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    293:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    294:   if (jac != B) {
    295:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    296:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
    297:   }
    298:   return 0;
    299: }

    301: PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx)
    302: {
    303:   PetscErrorCode        ierr;
    304:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
    305:   PetscViewer           viewer = monP->viewer;
    306:   SNESConvergedReason   reason;
    307:   const char            *strreason;

    309:   SNESGetConvergedReason(snes,&reason);
    310:   SNESGetConvergedReasonString(snes,&strreason);
    311:   PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n");
    312:   PetscViewerASCIIAddTab(viewer,1);
    313:   if (reason > 0) {
    314:     PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason);
    315:   } else if (reason <= 0) {
    316:     PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason);
    317:   }
    318:   PetscViewerASCIISubtractTab(viewer,1);
    319:   return 0;
    320: }

    322: PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx)
    323: {
    324:   PetscErrorCode        ierr;
    325:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
    326:   PetscViewer           viewer = monP->viewer;
    327:   KSPConvergedReason    reason;
    328:   const char            *reasonstr;

    330:   KSPGetConvergedReason(ksp,&reason);
    331:   KSPGetConvergedReasonString(ksp,&reasonstr);
    332:   PetscViewerASCIIAddTab(viewer,2);
    333:   PetscViewerASCIIPrintf(viewer,"Customized KSP converged reason view\n");
    334:   PetscViewerASCIIAddTab(viewer,1);
    335:   if (reason > 0) {
    336:     PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",reasonstr);
    337:   } else if (reason <= 0) {
    338:     PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",reasonstr);
    339:   }
    340:   PetscViewerASCIISubtractTab(viewer,3);
    341:   return 0;
    342: }

    345: /*TEST

    347:    test:
    348:       suffix: 1
    349:       nsize: 1

    351:    test:
    352:       suffix: 2
    353:       nsize: 1
    354:       args: -ksp_converged_reason_view_cancel

    356:    test:
    357:       suffix: 3
    358:       nsize: 1
    359:       args: -ksp_converged_reason_view_cancel -ksp_converged_reason

    361:    test:
    362:       suffix: 4
    363:       nsize: 1
    364:       args: -snes_converged_reason_view_cancel

    366:    test:
    367:       suffix: 5
    368:       nsize: 1
    369:       args: -snes_converged_reason_view_cancel -snes_converged_reason

    371: TEST*/