Actual source code: ex2.c

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

      4: /*
      5:    Include "petscdraw.h" so that we can use PETSc drawing routines.
      6:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
      7:    file automatically includes:
      8:      petscsys.h       - base PETSc routines   petscvec.h - vectors
      9:      petscmat.h - matrices
     10:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     11:      petscviewer.h - viewers               petscpc.h  - preconditioners
     12:      petscksp.h   - linear solvers
     13: */

     15: #include <petscsnes.h>

     17: /*
     18:    User-defined routines
     19: */
     20: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
     21: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
     22: extern PetscErrorCode FormInitialGuess(Vec);
     23: extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);

     25: /*
     26:    User-defined context for monitoring
     27: */
     28: typedef struct {
     29:   PetscViewer viewer;
     30: } MonitorCtx;

     32: int main(int argc, char **argv)
     33: {
     34:   SNES        snes;       /* SNES context */
     35:   Vec         x, r, F, U; /* vectors */
     36:   Mat         J;          /* Jacobian matrix */
     37:   MonitorCtx  monP;       /* monitoring context */
     38:   PetscInt    its, n = 5, i, maxit, maxf;
     39:   PetscMPIInt size;
     40:   PetscScalar h, xp, v, none = -1.0;
     41:   PetscReal   abstol, rtol, stol, norm;

     43:   PetscFunctionBeginUser;
     44:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     45:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     46:   PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
     47:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
     48:   h = 1.0 / (n - 1);

     50:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     51:      Create nonlinear solver context
     52:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     54:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

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

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

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

     78:   PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
     79:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
     80:   PetscCall(MatSetFromOptions(J));
     81:   PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));

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

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

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

    100:   /*
    101:      Set an optional user-defined monitoring routine
    102:   */
    103:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer));
    104:   PetscCall(SNESMonitorSet(snes, Monitor, &monP, 0));

    106:   /*
    107:      Set names for some vectors to facilitate monitoring (optional)
    108:   */
    109:   PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
    110:   PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));

    112:   /*
    113:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    114:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    115:   */
    116:   PetscCall(SNESSetFromOptions(snes));

    118:   /*
    119:      Print parameters used for convergence testing (optional) ... just
    120:      to demonstrate this routine; this information is also printed with
    121:      the option -snes_view
    122:   */
    123:   PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
    124:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));

    126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    127:      Initialize application:
    128:      Store right-hand side of PDE and exact solution
    129:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

    140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    141:      Evaluate initial guess; then solve nonlinear system
    142:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    143:   /*
    144:      Note: The user should initialize the vector, x, with the initial guess
    145:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    146:      to employ an initial guess of zero, the user should explicitly set
    147:      this vector to zero by calling VecSet().
    148:   */
    149:   PetscCall(FormInitialGuess(x));
    150:   PetscCall(SNESSolve(snes, NULL, x));
    151:   PetscCall(SNESGetIterationNumber(snes, &its));
    152:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));

    154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    155:      Check solution and clean up
    156:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    158:   /*
    159:      Check the error
    160:   */
    161:   PetscCall(VecAXPY(x, none, U));
    162:   PetscCall(VecNorm(x, NORM_2, &norm));
    163:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));

    165:   /*
    166:      Free work space.  All PETSc objects should be destroyed when they
    167:      are no longer needed.
    168:   */
    169:   PetscCall(VecDestroy(&x));
    170:   PetscCall(VecDestroy(&r));
    171:   PetscCall(VecDestroy(&U));
    172:   PetscCall(VecDestroy(&F));
    173:   PetscCall(MatDestroy(&J));
    174:   PetscCall(SNESDestroy(&snes));
    175:   PetscCall(PetscViewerDestroy(&monP.viewer));
    176:   PetscCall(PetscFinalize());
    177:   return 0;
    178: }
    179: /* ------------------------------------------------------------------- */
    180: /*
    181:    FormInitialGuess - Computes initial guess.

    183:    Input/Output Parameter:
    184: .  x - the solution vector
    185: */
    186: PetscErrorCode FormInitialGuess(Vec x)
    187: {
    188:   PetscFunctionBeginUser;
    189:   PetscCall(VecSet(x, 0.5));
    190:   PetscFunctionReturn(PETSC_SUCCESS);
    191: }
    192: /* ------------------------------------------------------------------- */
    193: /*
    194:    FormFunction - Evaluates nonlinear function, F(x).

    196:    Input Parameters:
    197: .  snes - the SNES context
    198: .  x - input vector
    199: .  ctx - optional user-defined context, as set by SNESSetFunction()

    201:    Output Parameter:
    202: .  f - function vector

    204:    Note:
    205:    The user-defined context can contain any application-specific data
    206:    needed for the function evaluation (such as various parameters, work
    207:    vectors, and grid information).  In this program the context is just
    208:    a vector containing the right-hand side of the discretized PDE.
    209:  */

    211: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
    212: {
    213:   Vec                g = (Vec)ctx;
    214:   const PetscScalar *xx, *gg;
    215:   PetscScalar       *ff, d;
    216:   PetscInt           i, n;

    218:   PetscFunctionBeginUser;
    219:   /*
    220:      Get pointers to vector data.
    221:        - For default PETSc vectors, VecGetArray() returns a pointer to
    222:          the data array.  Otherwise, the routine is implementation dependent.
    223:        - You MUST call VecRestoreArray() when you no longer need access to
    224:          the array.
    225:   */
    226:   PetscCall(VecGetArrayRead(x, &xx));
    227:   PetscCall(VecGetArray(f, &ff));
    228:   PetscCall(VecGetArrayRead(g, &gg));

    230:   /*
    231:      Compute function
    232:   */
    233:   PetscCall(VecGetSize(x, &n));
    234:   d     = (PetscReal)(n - 1);
    235:   d     = d * d;
    236:   ff[0] = xx[0];
    237:   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];
    238:   ff[n - 1] = xx[n - 1] - 1.0;

    240:   /*
    241:      Restore vectors
    242:   */
    243:   PetscCall(VecRestoreArrayRead(x, &xx));
    244:   PetscCall(VecRestoreArray(f, &ff));
    245:   PetscCall(VecRestoreArrayRead(g, &gg));
    246:   PetscFunctionReturn(PETSC_SUCCESS);
    247: }
    248: /* ------------------------------------------------------------------- */
    249: /*
    250:    FormJacobian - Evaluates Jacobian matrix.

    252:    Input Parameters:
    253: .  snes - the SNES context
    254: .  x - input vector
    255: .  dummy - optional user-defined context (not used here)

    257:    Output Parameters:
    258: .  jac - Jacobian matrix
    259: .  B - optionally different preconditioning matrix

    261: */

    263: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
    264: {
    265:   const PetscScalar *xx;
    266:   PetscScalar        A[3], d;
    267:   PetscInt           i, n, j[3];

    269:   PetscFunctionBeginUser;
    270:   /*
    271:      Get pointer to vector data
    272:   */
    273:   PetscCall(VecGetArrayRead(x, &xx));

    275:   /*
    276:      Compute Jacobian entries and insert into matrix.
    277:       - Note that in this case we set all elements for a particular
    278:         row at once.
    279:   */
    280:   PetscCall(VecGetSize(x, &n));
    281:   d = (PetscReal)(n - 1);
    282:   d = d * d;

    284:   /*
    285:      Interior grid points
    286:   */
    287:   for (i = 1; i < n - 1; i++) {
    288:     j[0] = i - 1;
    289:     j[1] = i;
    290:     j[2] = i + 1;
    291:     A[0] = A[2] = d;
    292:     A[1]        = -2.0 * d + 2.0 * xx[i];
    293:     PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
    294:   }

    296:   /*
    297:      Boundary points
    298:   */
    299:   i    = 0;
    300:   A[0] = 1.0;

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

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

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

    309:   /*
    310:      Restore vector
    311:   */
    312:   PetscCall(VecRestoreArrayRead(x, &xx));

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

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

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

    347:   PetscFunctionBeginUser;
    348:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm));
    349:   PetscCall(SNESGetConvergedReason(snes, &reason));
    350:   PetscCall(SNESGetSolution(snes, &x));
    351:   PetscCall(VecView(x, monP->viewer));
    352:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  converged = %s\n", SNESConvergedReasons[reason]));
    353:   PetscFunctionReturn(PETSC_SUCCESS);
    354: }

    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:
    371:       suffix: 4
    372:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view
    373:       requires: !single

    375:    test:
    376:       suffix: 5
    377:       filter: grep -v atol | sed -e "s/CONVERGED_ITS/DIVERGED_MAX_IT/g" | sed -e "s/CONVERGED_FNORM_RELATIVE/DIVERGED_MAX_IT/g"
    378:       args: -nox -snes_type {{newtonls newtontr ncg ngmres qn anderson nrichardson ms ksponly ksptransposeonly vinewtonrsls vinewtonssls fas ms}} -snes_max_it 1
    379:       requires: !single

    381: TEST*/