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: }