Actual source code: snesrichardson.c

      1: #include <../src/snes/impls/richardson/snesrichardsonimpl.h>

      3: static PetscErrorCode SNESDestroy_NRichardson(SNES snes)
      4: {
      5:   PetscFunctionBegin;
      6:   PetscCall(PetscFree(snes->data));
      7:   PetscFunctionReturn(PETSC_SUCCESS);
      8: }

     10: static PetscErrorCode SNESSetUp_NRichardson(SNES snes)
     11: {
     12:   PetscFunctionBegin;
     13:   PetscCheck(snes->npcside != PC_RIGHT, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "NRichardson only supports left preconditioning");
     14:   if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
     15:   PetscFunctionReturn(PETSC_SUCCESS);
     16: }

     18: static PetscErrorCode SNESSetFromOptions_NRichardson(SNES snes, PetscOptionItems PetscOptionsObject)
     19: {
     20:   PetscFunctionBegin;
     21:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES Richardson options");
     22:   PetscOptionsHeadEnd();
     23:   PetscFunctionReturn(PETSC_SUCCESS);
     24: }

     26: static PetscErrorCode SNESSolve_NRichardson(SNES snes)
     27: {
     28:   Vec                  X, Y, F;
     29:   PetscReal            xnorm, fnorm, ynorm;
     30:   PetscInt             maxits, i;
     31:   SNESLineSearchReason lsresult;
     32:   SNESConvergedReason  reason;

     34:   PetscFunctionBegin;
     35:   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

     37:   snes->reason = SNES_CONVERGED_ITERATING;

     39:   maxits = snes->max_its;        /* maximum number of iterations */
     40:   X      = snes->vec_sol;        /* X^n */
     41:   Y      = snes->vec_sol_update; /* \tilde X */
     42:   F      = snes->vec_func;       /* residual vector */

     44:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
     45:   snes->iter = 0;
     46:   snes->norm = 0.;
     47:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

     49:   if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
     50:     PetscCall(SNESApplyNPC(snes, X, NULL, F));
     51:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
     52:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
     53:       snes->reason = SNES_DIVERGED_INNER;
     54:       PetscFunctionReturn(PETSC_SUCCESS);
     55:     }
     56:     PetscCall(VecNorm(F, NORM_2, &fnorm));
     57:   } else {
     58:     if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
     59:     else snes->vec_func_init_set = PETSC_FALSE;

     61:     PetscCall(VecNorm(F, NORM_2, &fnorm));
     62:     SNESCheckFunctionNorm(snes, fnorm);
     63:   }
     64:   if (snes->npc && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
     65:     PetscCall(SNESApplyNPC(snes, X, F, Y));
     66:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
     67:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
     68:       snes->reason = SNES_DIVERGED_INNER;
     69:       PetscFunctionReturn(PETSC_SUCCESS);
     70:     }
     71:   } else {
     72:     PetscCall(VecCopy(F, Y));
     73:   }

     75:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
     76:   snes->norm = fnorm;
     77:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
     78:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));

     80:   /* test convergence */
     81:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
     82:   PetscCall(SNESMonitor(snes, 0, fnorm));
     83:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

     85:   /* Call general purpose update function */
     86:   PetscTryTypeMethod(snes, update, snes->iter);

     88:   for (i = 1; i < maxits + 1; i++) {
     89:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y));
     90:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsresult));
     91:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
     92:     if (lsresult) {
     93:       if (++snes->numFailures >= snes->maxFailures) {
     94:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
     95:         break;
     96:       }
     97:     }
     98:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
     99:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    100:       break;
    101:     }

    103:     /* Monitor convergence */
    104:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    105:     snes->iter  = i;
    106:     snes->norm  = fnorm;
    107:     snes->xnorm = xnorm;
    108:     snes->ynorm = ynorm;
    109:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    110:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
    111:     /* Test for convergence */
    112:     PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
    113:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
    114:     if (snes->reason) break;

    116:     /* Call general purpose update function */
    117:     PetscTryTypeMethod(snes, update, snes->iter);

    119:     if (snes->npc) {
    120:       if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
    121:         PetscCall(SNESApplyNPC(snes, X, NULL, Y));
    122:         PetscCall(VecNorm(F, NORM_2, &fnorm));
    123:         PetscCall(VecCopy(Y, F));
    124:       } else {
    125:         PetscCall(SNESApplyNPC(snes, X, F, Y));
    126:       }
    127:       PetscCall(SNESGetConvergedReason(snes->npc, &reason));
    128:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    129:         snes->reason = SNES_DIVERGED_INNER;
    130:         PetscFunctionReturn(PETSC_SUCCESS);
    131:       }
    132:     } else {
    133:       PetscCall(VecCopy(F, Y));
    134:     }
    135:   }
    136:   if (i == maxits + 1) {
    137:     PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
    138:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
    139:   }
    140:   PetscFunctionReturn(PETSC_SUCCESS);
    141: }

    143: /*MC
    144:    SNESNRICHARDSON - Richardson nonlinear solver that uses successive substitutions, also sometimes known as Picard iteration.

    146:    Options Database Keys:
    147: +  -snes_linesearch_type <l2,cp,basic> - Line search type.
    148: -  -snes_linesearch_damping <1.0>      - Damping for the line search.

    150:    Level: beginner

    152:    Notes:
    153:    If no inner nonlinear preconditioner is provided then solves $F(x) - b = 0 $ using $x^{n+1} = x^{n} - \lambda
    154:    (F(x^n) - b) $ where $ \lambda$ is obtained with either `SNESLineSearchSetDamping()`, `-snes_damping` or a line search.  If
    155:    an inner nonlinear preconditioner is provided (either with `-npc_snes_typ`e or `SNESSetNPC()`) then the inner
    156:    solver is called on the initial solution $x^n$ and the nonlinear Richardson uses $ x^{n+1} = x^{n} + \lambda d^{n}$
    157:    where $d^{n} = \hat{x}^{n} - x^{n} $ where $\hat{x}^{n} $ is the solution returned from the inner solver.

    159:    The update, especially without inner nonlinear preconditioner, may be ill-scaled.  If using the basic
    160:    linesearch, one may have to scale the update with `-snes_linesearch_damping`

    162:    This uses no derivative information provided with `SNESSetJacobian()` thus it will be much slower than Newton's method obtained with `-snes_type ls`

    164:    Only supports left non-linear preconditioning.

    166: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESNCG`,
    167:           `SNESLineSearchSetDamping()`
    168: M*/
    169: PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes)
    170: {
    171:   SNES_NRichardson *neP;
    172:   SNESLineSearch    linesearch;

    174:   PetscFunctionBegin;
    175:   snes->ops->destroy        = SNESDestroy_NRichardson;
    176:   snes->ops->setup          = SNESSetUp_NRichardson;
    177:   snes->ops->setfromoptions = SNESSetFromOptions_NRichardson;
    178:   snes->ops->solve          = SNESSolve_NRichardson;

    180:   snes->usesksp = PETSC_FALSE;
    181:   snes->usesnpc = PETSC_TRUE;

    183:   snes->npcside = PC_LEFT;

    185:   PetscCall(SNESGetLineSearch(snes, &linesearch));
    186:   if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT));

    188:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    190:   PetscCall(SNESParametersInitialize(snes));
    191:   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
    192:   PetscObjectParameterSetDefault(snes, max_its, 10000);
    193:   PetscObjectParameterSetDefault(snes, stol, 1e-20);

    195:   PetscCall(PetscNew(&neP));
    196:   snes->data = (void *)neP;
    197:   PetscFunctionReturn(PETSC_SUCCESS);
    198: }