Actual source code: snesrichardson.c

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

  3: PetscErrorCode SNESReset_NRichardson(SNES snes)
  4: {
  5:   return 0;
  6: }

  8: /*
  9:   SNESDestroy_NRichardson - Destroys the private SNES_NRichardson context that was created with SNESCreate_NRichardson().

 11:   Input Parameter:
 12: . snes - the SNES context

 14:   Application Interface Routine: SNESDestroy()
 15: */
 16: PetscErrorCode SNESDestroy_NRichardson(SNES snes)
 17: {
 18:   SNESReset_NRichardson(snes);
 19:   PetscFree(snes->data);
 20:   return 0;
 21: }

 23: /*
 24:    SNESSetUp_NRichardson - Sets up the internal data structures for the later use
 25:    of the SNESNRICHARDSON nonlinear solver.

 27:    Input Parameters:
 28: +  snes - the SNES context
 29: -  x - the solution vector

 31:    Application Interface Routine: SNESSetUp()
 32:  */
 33: PetscErrorCode SNESSetUp_NRichardson(SNES snes)
 34: {
 36:   if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
 37:   return 0;
 38: }

 40: /*
 41:   SNESSetFromOptions_NRichardson - Sets various parameters for the SNESNEWTONLS method.

 43:   Input Parameter:
 44: . snes - the SNES context

 46:   Application Interface Routine: SNESSetFromOptions()
 47: */
 48: static PetscErrorCode SNESSetFromOptions_NRichardson(SNES snes, PetscOptionItems *PetscOptionsObject)
 49: {
 50:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES Richardson options");
 51:   PetscOptionsHeadEnd();
 52:   return 0;
 53: }

 55: /*
 56:   SNESView_NRichardson - Prints info from the SNESRichardson data structure.

 58:   Input Parameters:
 59: + SNES - the SNES context
 60: - viewer - visualization context

 62:   Application Interface Routine: SNESView()
 63: */
 64: static PetscErrorCode SNESView_NRichardson(SNES snes, PetscViewer viewer)
 65: {
 66:   PetscBool iascii;

 68:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 69:   if (iascii) { }
 70:   return 0;
 71: }

 73: /*
 74:   SNESSolve_NRichardson - Solves a nonlinear system with the Richardson method.

 76:   Input Parameters:
 77: . snes - the SNES context

 79:   Output Parameter:
 80: . outits - number of iterations until termination

 82:   Application Interface Routine: SNESSolve()
 83: */
 84: PetscErrorCode SNESSolve_NRichardson(SNES snes)
 85: {
 86:   Vec                  X, Y, F;
 87:   PetscReal            xnorm, fnorm, ynorm;
 88:   PetscInt             maxits, i;
 89:   SNESLineSearchReason lsresult;
 90:   SNESConvergedReason  reason;


 94:   snes->reason = SNES_CONVERGED_ITERATING;

 96:   maxits = snes->max_its;        /* maximum number of iterations */
 97:   X      = snes->vec_sol;        /* X^n */
 98:   Y      = snes->vec_sol_update; /* \tilde X */
 99:   F      = snes->vec_func;       /* residual vector */

101:   PetscObjectSAWsTakeAccess((PetscObject)snes);
102:   snes->iter = 0;
103:   snes->norm = 0.;
104:   PetscObjectSAWsGrantAccess((PetscObject)snes);

106:   if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
107:     SNESApplyNPC(snes, X, NULL, F);
108:     SNESGetConvergedReason(snes->npc, &reason);
109:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
110:       snes->reason = SNES_DIVERGED_INNER;
111:       return 0;
112:     }
113:     VecNorm(F, NORM_2, &fnorm);
114:   } else {
115:     if (!snes->vec_func_init_set) {
116:       SNESComputeFunction(snes, X, F);
117:     } else snes->vec_func_init_set = PETSC_FALSE;

119:     VecNorm(F, NORM_2, &fnorm);
120:     SNESCheckFunctionNorm(snes, fnorm);
121:   }
122:   if (snes->npc && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
123:     SNESApplyNPC(snes, X, F, Y);
124:     SNESGetConvergedReason(snes->npc, &reason);
125:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
126:       snes->reason = SNES_DIVERGED_INNER;
127:       return 0;
128:     }
129:   } else {
130:     VecCopy(F, Y);
131:   }

133:   PetscObjectSAWsTakeAccess((PetscObject)snes);
134:   snes->norm = fnorm;
135:   PetscObjectSAWsGrantAccess((PetscObject)snes);
136:   SNESLogConvergenceHistory(snes, fnorm, 0);
137:   SNESMonitor(snes, 0, fnorm);

139:   /* test convergence */
140:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
141:   if (snes->reason) return 0;

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

146:   /* set parameter for default relative tolerance convergence test */
147:   snes->ttol = fnorm * snes->rtol;
148:   /* test convergence */
149:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
150:   if (snes->reason) return 0;

152:   for (i = 1; i < maxits + 1; i++) {
153:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
154:     SNESLineSearchGetReason(snes->linesearch, &lsresult);
155:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
156:     if (lsresult) {
157:       if (++snes->numFailures >= snes->maxFailures) {
158:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
159:         break;
160:       }
161:     }
162:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
163:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
164:       break;
165:     }

167:     /* Monitor convergence */
168:     PetscObjectSAWsTakeAccess((PetscObject)snes);
169:     snes->iter  = i;
170:     snes->norm  = fnorm;
171:     snes->xnorm = xnorm;
172:     snes->ynorm = ynorm;
173:     PetscObjectSAWsGrantAccess((PetscObject)snes);
174:     SNESLogConvergenceHistory(snes, snes->norm, 0);
175:     SNESMonitor(snes, snes->iter, snes->norm);
176:     /* Test for convergence */
177:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
178:     if (snes->reason) break;

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

183:     if (snes->npc) {
184:       if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
185:         SNESApplyNPC(snes, X, NULL, Y);
186:         VecNorm(F, NORM_2, &fnorm);
187:         VecCopy(Y, F);
188:       } else {
189:         SNESApplyNPC(snes, X, F, Y);
190:       }
191:       SNESGetConvergedReason(snes->npc, &reason);
192:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
193:         snes->reason = SNES_DIVERGED_INNER;
194:         return 0;
195:       }
196:     } else {
197:       VecCopy(F, Y);
198:     }
199:   }
200:   if (i == maxits + 1) {
201:     PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
202:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
203:   }
204:   return 0;
205: }

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

210:    Options Database Keys:
211: +  -snes_linesearch_type <l2,cp,basic> - Line search type.
212: -  -snes_linesearch_damping<1.0> - Damping for the line search.

214:    Level: beginner

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

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

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

228:    Only supports left non-linear preconditioning.

230: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESNCG`
231: M*/
232: PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes)
233: {
234:   SNES_NRichardson *neP;
235:   SNESLineSearch    linesearch;

237:   snes->ops->destroy        = SNESDestroy_NRichardson;
238:   snes->ops->setup          = SNESSetUp_NRichardson;
239:   snes->ops->setfromoptions = SNESSetFromOptions_NRichardson;
240:   snes->ops->view           = SNESView_NRichardson;
241:   snes->ops->solve          = SNESSolve_NRichardson;
242:   snes->ops->reset          = SNESReset_NRichardson;

244:   snes->usesksp = PETSC_FALSE;
245:   snes->usesnpc = PETSC_TRUE;

247:   snes->npcside = PC_LEFT;

249:   SNESGetLineSearch(snes, &linesearch);
250:   if (!((PetscObject)linesearch)->type_name) SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);

252:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

254:   PetscNew(&neP);
255:   snes->data = (void *)neP;

257:   if (!snes->tolerancesset) {
258:     snes->max_funcs = 30000;
259:     snes->max_its   = 10000;
260:     snes->stol      = 1e-20;
261:   }
262:   return 0;
263: }