Actual source code: ex6.c

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

      4: #include <petscsnes.h>

      6: /*
      7:    User-defined routines
      8: */
      9: PETSC_EXTERN PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
     10: PETSC_EXTERN PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
     11: PETSC_EXTERN PetscErrorCode FormInitialGuess(Vec);
     12: PETSC_EXTERN PetscErrorCode MySNESConvergedReasonView(SNES, void *);
     13: PETSC_EXTERN PetscErrorCode MyKSPConvergedReasonView(KSP, void *);

     15: /*
     16:    User-defined context for monitoring
     17: */
     18: typedef struct {
     19:   PetscViewer viewer;
     20: } ReasonViewCtx;

     22: int main(int argc, char **argv)
     23: {
     24:   SNES          snes;       /* SNES context */
     25:   KSP           ksp;        /* KSP context */
     26:   Vec           x, r, F, U; /* vectors */
     27:   Mat           J;          /* Jacobian matrix */
     28:   ReasonViewCtx monP;       /* monitoring context */
     29:   PetscInt      its, n = 5, i;
     30:   PetscMPIInt   size;
     31:   PetscScalar   h, xp, v;
     32:   MPI_Comm      comm;

     34:   PetscFunctionBeginUser;
     35:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     36:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     37:   PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
     38:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
     39:   h    = 1.0 / (n - 1);
     40:   comm = PETSC_COMM_WORLD;
     41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     42:      Create nonlinear solver context
     43:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     45:   PetscCall(SNESCreate(comm, &snes));

     47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     48:      Create vector data structures; set function evaluation routine
     49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     50:   /*
     51:      Note that we form 1 vector from scratch and then duplicate as needed.
     52:   */
     53:   PetscCall(VecCreate(comm, &x));
     54:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
     55:   PetscCall(VecSetFromOptions(x));
     56:   PetscCall(VecDuplicate(x, &r));
     57:   PetscCall(VecDuplicate(x, &F));
     58:   PetscCall(VecDuplicate(x, &U));

     60:   /*
     61:      Set function evaluation routine and vector
     62:   */
     63:   PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));

     65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     66:      Create matrix data structure; set Jacobian evaluation routine
     67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     69:   PetscCall(MatCreate(comm, &J));
     70:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
     71:   PetscCall(MatSetFromOptions(J));
     72:   PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));

     74:   /*
     75:      Set Jacobian matrix data structure and default Jacobian evaluation
     76:      routine. User can override with:
     77:      -snes_fd : default finite differencing approximation of Jacobian
     78:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
     79:                 (unless user explicitly sets preconditioner)
     80:      -snes_mf_operator : form preconditioning matrix as set by the user,
     81:                          but use matrix-free approx for Jacobian-vector
     82:                          products within Newton-Krylov method
     83:   */

     85:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));

     87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     88:      Customize nonlinear solver; set runtime options
     89:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     91:   /*
     92:      Set an optional user-defined reasonview routine
     93:   */
     94:   PetscCall(PetscViewerASCIIGetStdout(comm, &monP.viewer));
     95:   /* Just make sure we can not repeat adding the same function
     96:    * PETSc will be able to ignore the repeated function
     97:    */
     98:   for (i = 0; i < 4; i++) PetscCall(SNESConvergedReasonViewSet(snes, MySNESConvergedReasonView, &monP, 0));
     99:   PetscCall(SNESGetKSP(snes, &ksp));
    100:   /* Just make sure we can not repeat adding the same function
    101:    * PETSc will be able to ignore the repeated function
    102:    */
    103:   for (i = 0; i < 4; i++) PetscCall(KSPConvergedReasonViewSet(ksp, MyKSPConvergedReasonView, &monP, 0));
    104:   /*
    105:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    106:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    107:   */
    108:   PetscCall(SNESSetFromOptions(snes));

    110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    111:      Initialize application:
    112:      Store right-hand side of PDE and exact solution
    113:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    115:   xp = 0.0;
    116:   for (i = 0; i < n; i++) {
    117:     v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
    118:     PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
    119:     v = xp * xp * xp;
    120:     PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES));
    121:     xp += h;
    122:   }

    124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    125:      Evaluate initial guess; then solve nonlinear system
    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:   */
    133:   PetscCall(FormInitialGuess(x));
    134:   PetscCall(SNESSolve(snes, NULL, x));
    135:   PetscCall(SNESGetIterationNumber(snes, &its));

    137:   /*
    138:      Free work space.  All PETSc objects should be destroyed when they
    139:      are no longer needed.
    140:   */
    141:   PetscCall(VecDestroy(&x));
    142:   PetscCall(VecDestroy(&r));
    143:   PetscCall(VecDestroy(&U));
    144:   PetscCall(VecDestroy(&F));
    145:   PetscCall(MatDestroy(&J));
    146:   PetscCall(SNESDestroy(&snes));
    147:   PetscCall(PetscFinalize());
    148:   return 0;
    149: }

    151: /*
    152:    FormInitialGuess - Computes initial guess.

    154:    Input/Output Parameter:
    155: .  x - the solution vector
    156: */
    157: PetscErrorCode FormInitialGuess(Vec x)
    158: {
    159:   PetscScalar pfive = .50;

    161:   PetscFunctionBeginUser;
    162:   PetscCall(VecSet(x, pfive));
    163:   PetscFunctionReturn(PETSC_SUCCESS);
    164: }

    166: /*
    167:    FormFunction - Evaluates nonlinear function, F(x).

    169:    Input Parameters:
    170: .  snes - the SNES context
    171: .  x - input vector
    172: .  ctx - optional user-defined context, as set by SNESSetFunction()

    174:    Output Parameter:
    175: .  f - function vector

    177:    Note:
    178:    The user-defined context can contain any application-specific data
    179:    needed for the function evaluation (such as various parameters, work
    180:    vectors, and grid information).  In this program the context is just
    181:    a vector containing the right-hand side of the discretized PDE.
    182:  */

    184: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
    185: {
    186:   Vec                g = (Vec)ctx;
    187:   const PetscScalar *xx, *gg;
    188:   PetscScalar       *ff, d;
    189:   PetscInt           i, n;

    191:   PetscFunctionBeginUser;
    192:   /*
    193:      Get pointers to vector data.
    194:        - For default PETSc vectors, VecGetArray() returns a pointer to
    195:          the data array.  Otherwise, the routine is implementation dependent.
    196:        - You MUST call VecRestoreArray() when you no longer need access to
    197:          the array.
    198:   */
    199:   PetscCall(VecGetArrayRead(x, &xx));
    200:   PetscCall(VecGetArray(f, &ff));
    201:   PetscCall(VecGetArrayRead(g, &gg));

    203:   /*
    204:      Compute function
    205:   */
    206:   PetscCall(VecGetSize(x, &n));
    207:   d     = (PetscReal)(n - 1);
    208:   d     = d * d;
    209:   ff[0] = xx[0];
    210:   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];
    211:   ff[n - 1] = xx[n - 1] - 1.0;

    213:   /*
    214:      Restore vectors
    215:   */
    216:   PetscCall(VecRestoreArrayRead(x, &xx));
    217:   PetscCall(VecRestoreArray(f, &ff));
    218:   PetscCall(VecRestoreArrayRead(g, &gg));
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }
    221: /* ------------------------------------------------------------------- */
    222: /*
    223:    FormJacobian - Evaluates Jacobian matrix.

    225:    Input Parameters:
    226: .  snes - the SNES context
    227: .  x - input vector
    228: .  dummy - optional user-defined context (not used here)

    230:    Output Parameters:
    231: .  jac - Jacobian matrix
    232: .  B - optionally different preconditioning matrix

    234: */

    236: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
    237: {
    238:   const PetscScalar *xx;
    239:   PetscScalar        A[3], d;
    240:   PetscInt           i, n, j[3];

    242:   PetscFunctionBeginUser;
    243:   /*
    244:      Get pointer to vector data
    245:   */
    246:   PetscCall(VecGetArrayRead(x, &xx));

    248:   /*
    249:      Compute Jacobian entries and insert into matrix.
    250:       - Note that in this case we set all elements for a particular
    251:         row at once.
    252:   */
    253:   PetscCall(VecGetSize(x, &n));
    254:   d = (PetscReal)(n - 1);
    255:   d = d * d;

    257:   /*
    258:      Interior grid points
    259:   */
    260:   for (i = 1; i < n - 1; i++) {
    261:     j[0] = i - 1;
    262:     j[1] = i;
    263:     j[2] = i + 1;
    264:     A[0] = A[2] = d;
    265:     A[1]        = -2.0 * d + 2.0 * xx[i];
    266:     PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
    267:   }

    269:   /*
    270:      Boundary points
    271:   */
    272:   i    = 0;
    273:   A[0] = 1.0;

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

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

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

    282:   /*
    283:      Restore vector
    284:   */
    285:   PetscCall(VecRestoreArrayRead(x, &xx));

    287:   /*
    288:      Assemble matrix
    289:   */
    290:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    291:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    292:   if (jac != B) {
    293:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    294:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
    295:   }
    296:   PetscFunctionReturn(PETSC_SUCCESS);
    297: }

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

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

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

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

    342: /*TEST

    344:    test:
    345:       suffix: 1
    346:       nsize: 1
    347:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"

    349:    test:
    350:       suffix: 2
    351:       nsize: 1
    352:       args: -ksp_converged_reason_view_cancel
    353:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"

    355:    test:
    356:       suffix: 3
    357:       nsize: 1
    358:       args: -ksp_converged_reason_view_cancel -ksp_converged_reason
    359:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"

    361:    test:
    362:       suffix: 4
    363:       nsize: 1
    364:       args: -snes_converged_reason_view_cancel
    365:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"

    367:    test:
    368:       suffix: 5
    369:       nsize: 1
    370:       args: -snes_converged_reason_view_cancel -snes_converged_reason
    371:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"

    373: TEST*/