Actual source code: snesrichardson.c

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

  3: PetscErrorCode SNESReset_NRichardson(SNES snes)
  4: {
  5:   PetscFunctionBegin;
  6:   PetscFunctionReturn(PETSC_SUCCESS);
  7: }

  9: PetscErrorCode SNESDestroy_NRichardson(SNES snes)
 10: {
 11:   PetscFunctionBegin;
 12:   PetscCall(SNESReset_NRichardson(snes));
 13:   PetscCall(PetscFree(snes->data));
 14:   PetscFunctionReturn(PETSC_SUCCESS);
 15: }

 17: PetscErrorCode SNESSetUp_NRichardson(SNES snes)
 18: {
 19:   PetscFunctionBegin;
 20:   PetscCheck(snes->npcside != PC_RIGHT, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "NRichardson only supports left preconditioning");
 21:   if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
 22:   PetscFunctionReturn(PETSC_SUCCESS);
 23: }

 25: static PetscErrorCode SNESSetFromOptions_NRichardson(SNES snes, PetscOptionItems *PetscOptionsObject)
 26: {
 27:   PetscFunctionBegin;
 28:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES Richardson options");
 29:   PetscOptionsHeadEnd();
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: static PetscErrorCode SNESView_NRichardson(SNES snes, PetscViewer viewer)
 34: {
 35:   PetscBool iascii;

 37:   PetscFunctionBegin;
 38:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 39:   if (iascii) { }
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: PetscErrorCode SNESSolve_NRichardson(SNES snes)
 44: {
 45:   Vec                  X, Y, F;
 46:   PetscReal            xnorm, fnorm, ynorm;
 47:   PetscInt             maxits, i;
 48:   SNESLineSearchReason lsresult;
 49:   SNESConvergedReason  reason;

 51:   PetscFunctionBegin;
 52:   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);

 54:   snes->reason = SNES_CONVERGED_ITERATING;

 56:   maxits = snes->max_its;        /* maximum number of iterations */
 57:   X      = snes->vec_sol;        /* X^n */
 58:   Y      = snes->vec_sol_update; /* \tilde X */
 59:   F      = snes->vec_func;       /* residual vector */

 61:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
 62:   snes->iter = 0;
 63:   snes->norm = 0.;
 64:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

 66:   if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
 67:     PetscCall(SNESApplyNPC(snes, X, NULL, F));
 68:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
 69:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 70:       snes->reason = SNES_DIVERGED_INNER;
 71:       PetscFunctionReturn(PETSC_SUCCESS);
 72:     }
 73:     PetscCall(VecNorm(F, NORM_2, &fnorm));
 74:   } else {
 75:     if (!snes->vec_func_init_set) {
 76:       PetscCall(SNESComputeFunction(snes, X, F));
 77:     } else snes->vec_func_init_set = PETSC_FALSE;

 79:     PetscCall(VecNorm(F, NORM_2, &fnorm));
 80:     SNESCheckFunctionNorm(snes, fnorm);
 81:   }
 82:   if (snes->npc && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
 83:     PetscCall(SNESApplyNPC(snes, X, F, Y));
 84:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
 85:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 86:       snes->reason = SNES_DIVERGED_INNER;
 87:       PetscFunctionReturn(PETSC_SUCCESS);
 88:     }
 89:   } else {
 90:     PetscCall(VecCopy(F, Y));
 91:   }

 93:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
 94:   snes->norm = fnorm;
 95:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
 96:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
 97:   PetscCall(SNESMonitor(snes, 0, fnorm));

 99:   /* test convergence */
100:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
101:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

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

106:   /* set parameter for default relative tolerance convergence test */
107:   snes->ttol = fnorm * snes->rtol;
108:   /* test convergence */
109:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
110:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

112:   for (i = 1; i < maxits + 1; i++) {
113:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y));
114:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsresult));
115:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
116:     if (lsresult) {
117:       if (++snes->numFailures >= snes->maxFailures) {
118:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
119:         break;
120:       }
121:     }
122:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
123:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
124:       break;
125:     }

127:     /* Monitor convergence */
128:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
129:     snes->iter  = i;
130:     snes->norm  = fnorm;
131:     snes->xnorm = xnorm;
132:     snes->ynorm = ynorm;
133:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
134:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
135:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
136:     /* Test for convergence */
137:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
138:     if (snes->reason) break;

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

143:     if (snes->npc) {
144:       if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
145:         PetscCall(SNESApplyNPC(snes, X, NULL, Y));
146:         PetscCall(VecNorm(F, NORM_2, &fnorm));
147:         PetscCall(VecCopy(Y, F));
148:       } else {
149:         PetscCall(SNESApplyNPC(snes, X, F, Y));
150:       }
151:       PetscCall(SNESGetConvergedReason(snes->npc, &reason));
152:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
153:         snes->reason = SNES_DIVERGED_INNER;
154:         PetscFunctionReturn(PETSC_SUCCESS);
155:       }
156:     } else {
157:       PetscCall(VecCopy(F, Y));
158:     }
159:   }
160:   if (i == maxits + 1) {
161:     PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
162:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
163:   }
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

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

170:    Options Database Keys:
171: +  -snes_linesearch_type <l2,cp,basic> - Line search type.
172: -  -snes_linesearch_damping<1.0> - Damping for the line search.

174:    Level: beginner

176:    Notes:
177:    If no inner nonlinear preconditioner is provided then solves F(x) - b = 0 using x^{n+1} = x^{n} - lambda
178:    (F(x^n) - b) where lambda is obtained either `SNESLineSearchSetDamping()`, -snes_damping or a line search.  If
179:    an inner nonlinear preconditioner is provided (either with -npc_snes_type or `SNESSetNPC()`) then the inner
180:    solver is called an initial solution x^n and the nonlinear Richardson uses x^{n+1} = x^{n} + lambda d^{n}
181:    where d^{n} = \hat{x}^{n} - x^{n} where \hat{x}^{n} is the solution returned from the inner solver.

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

186:    This uses no derivative information thus will be much slower then Newton's method obtained with -snes_type ls

188:    Only supports left non-linear preconditioning.

190: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESNCG`
191: M*/
192: PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes)
193: {
194:   SNES_NRichardson *neP;
195:   SNESLineSearch    linesearch;

197:   PetscFunctionBegin;
198:   snes->ops->destroy        = SNESDestroy_NRichardson;
199:   snes->ops->setup          = SNESSetUp_NRichardson;
200:   snes->ops->setfromoptions = SNESSetFromOptions_NRichardson;
201:   snes->ops->view           = SNESView_NRichardson;
202:   snes->ops->solve          = SNESSolve_NRichardson;
203:   snes->ops->reset          = SNESReset_NRichardson;

205:   snes->usesksp = PETSC_FALSE;
206:   snes->usesnpc = PETSC_TRUE;

208:   snes->npcside = PC_LEFT;

210:   PetscCall(SNESGetLineSearch(snes, &linesearch));
211:   if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHL2));

213:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

215:   PetscCall(PetscNew(&neP));
216:   snes->data = (void *)neP;

218:   if (!snes->tolerancesset) {
219:     snes->max_funcs = 30000;
220:     snes->max_its   = 10000;
221:     snes->stol      = 1e-20;
222:   }
223:   PetscFunctionReturn(PETSC_SUCCESS);
224: }