Actual source code: snesrichardson.c
1: #include <../src/snes/impls/richardson/snesrichardsonimpl.h>
3: PetscErrorCode SNESReset_NRichardson(SNES snes)
4: {
6: return(0);
7: }
9: /*
10: SNESDestroy_NRichardson - Destroys the private SNES_NRichardson context that was created with SNESCreate_NRichardson().
12: Input Parameter:
13: . snes - the SNES context
15: Application Interface Routine: SNESDestroy()
16: */
17: PetscErrorCode SNESDestroy_NRichardson(SNES snes)
18: {
22: SNESReset_NRichardson(snes);
23: PetscFree(snes->data);
24: return(0);
25: }
27: /*
28: SNESSetUp_NRichardson - Sets up the internal data structures for the later use
29: of the SNESNRICHARDSON nonlinear solver.
31: Input Parameters:
32: + snes - the SNES context
33: - x - the solution vector
35: Application Interface Routine: SNESSetUp()
36: */
37: PetscErrorCode SNESSetUp_NRichardson(SNES snes)
38: {
40: if (snes->npcside== PC_RIGHT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"NRichardson only supports left preconditioning");
41: if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
42: return(0);
43: }
45: /*
46: SNESSetFromOptions_NRichardson - Sets various parameters for the SNESNEWTONLS method.
48: Input Parameter:
49: . snes - the SNES context
51: Application Interface Routine: SNESSetFromOptions()
52: */
53: static PetscErrorCode SNESSetFromOptions_NRichardson(PetscOptionItems *PetscOptionsObject,SNES snes)
54: {
58: PetscOptionsHead(PetscOptionsObject,"SNES Richardson options");
59: PetscOptionsTail();
60: return(0);
61: }
63: /*
64: SNESView_NRichardson - Prints info from the SNESRichardson data structure.
66: Input Parameters:
67: + SNES - the SNES context
68: - viewer - visualization context
70: Application Interface Routine: SNESView()
71: */
72: static PetscErrorCode SNESView_NRichardson(SNES snes, PetscViewer viewer)
73: {
74: PetscBool iascii;
78: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
79: if (iascii) {
80: }
81: return(0);
82: }
84: /*
85: SNESSolve_NRichardson - Solves a nonlinear system with the Richardson method.
87: Input Parameters:
88: . snes - the SNES context
90: Output Parameter:
91: . outits - number of iterations until termination
93: Application Interface Routine: SNESSolve()
94: */
95: PetscErrorCode SNESSolve_NRichardson(SNES snes)
96: {
97: Vec X, Y, F;
98: PetscReal xnorm, fnorm, ynorm;
99: PetscInt maxits, i;
100: PetscErrorCode ierr;
101: SNESLineSearchReason lsresult;
102: SNESConvergedReason reason;
105: if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
107: snes->reason = SNES_CONVERGED_ITERATING;
109: maxits = snes->max_its; /* maximum number of iterations */
110: X = snes->vec_sol; /* X^n */
111: Y = snes->vec_sol_update; /* \tilde X */
112: F = snes->vec_func; /* residual vector */
114: PetscObjectSAWsTakeAccess((PetscObject)snes);
115: snes->iter = 0;
116: snes->norm = 0.;
117: PetscObjectSAWsGrantAccess((PetscObject)snes);
119: if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
120: SNESApplyNPC(snes,X,NULL,F);
121: SNESGetConvergedReason(snes->npc,&reason);
122: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
123: snes->reason = SNES_DIVERGED_INNER;
124: return(0);
125: }
126: VecNorm(F,NORM_2,&fnorm);
127: } else {
128: if (!snes->vec_func_init_set) {
129: SNESComputeFunction(snes,X,F);
130: } else snes->vec_func_init_set = PETSC_FALSE;
132: VecNorm(F,NORM_2,&fnorm);
133: SNESCheckFunctionNorm(snes,fnorm);
134: }
135: if (snes->npc && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
136: SNESApplyNPC(snes,X,F,Y);
137: SNESGetConvergedReason(snes->npc,&reason);
138: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
139: snes->reason = SNES_DIVERGED_INNER;
140: return(0);
141: }
142: } else {
143: VecCopy(F,Y);
144: }
146: PetscObjectSAWsTakeAccess((PetscObject)snes);
147: snes->norm = fnorm;
148: PetscObjectSAWsGrantAccess((PetscObject)snes);
149: SNESLogConvergenceHistory(snes,fnorm,0);
150: SNESMonitor(snes,0,fnorm);
152: /* test convergence */
153: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
154: if (snes->reason) return(0);
156: /* Call general purpose update function */
157: if (snes->ops->update) {
158: (*snes->ops->update)(snes, snes->iter);
159: }
161: /* set parameter for default relative tolerance convergence test */
162: snes->ttol = fnorm*snes->rtol;
163: /* test convergence */
164: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
165: if (snes->reason) return(0);
167: for (i = 1; i < maxits+1; i++) {
168: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
169: SNESLineSearchGetReason(snes->linesearch, &lsresult);
170: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
171: if (lsresult) {
172: if (++snes->numFailures >= snes->maxFailures) {
173: snes->reason = SNES_DIVERGED_LINE_SEARCH;
174: break;
175: }
176: }
177: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
178: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
179: break;
180: }
182: /* Monitor convergence */
183: PetscObjectSAWsTakeAccess((PetscObject)snes);
184: snes->iter = i;
185: snes->norm = fnorm;
186: snes->xnorm = xnorm;
187: snes->ynorm = ynorm;
188: PetscObjectSAWsGrantAccess((PetscObject)snes);
189: SNESLogConvergenceHistory(snes,snes->norm,0);
190: SNESMonitor(snes,snes->iter,snes->norm);
191: /* Test for convergence */
192: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
193: if (snes->reason) break;
195: /* Call general purpose update function */
196: if (snes->ops->update) {
197: (*snes->ops->update)(snes, snes->iter);
198: }
200: if (snes->npc) {
201: if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
202: SNESApplyNPC(snes,X,NULL,Y);
203: VecNorm(F,NORM_2,&fnorm);
204: VecCopy(Y,F);
205: } else {
206: SNESApplyNPC(snes,X,F,Y);
207: }
208: SNESGetConvergedReason(snes->npc,&reason);
209: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
210: snes->reason = SNES_DIVERGED_INNER;
211: return(0);
212: }
213: } else {
214: VecCopy(F,Y);
215: }
216: }
217: if (i == maxits+1) {
218: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
219: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
220: }
221: return(0);
222: }
224: /*MC
225: SNESNRICHARDSON - Richardson nonlinear solver that uses successive substitutions, also sometimes known as Picard iteration.
227: Level: beginner
229: Options Database:
230: + -snes_linesearch_type <l2,cp,basic> Line search type.
231: - -snes_linesearch_damping<1.0> Damping for the line search.
233: Notes:
234: If no inner nonlinear preconditioner is provided then solves F(x) - b = 0 using x^{n+1} = x^{n} - lambda
235: (F(x^n) - b) where lambda is obtained either SNESLineSearchSetDamping(), -snes_damping or a line search. If
236: an inner nonlinear preconditioner is provided (either with -npc_snes_type or SNESSetNPC()) then the inner
237: solver is called an initial solution x^n and the nonlinear Richardson uses x^{n+1} = x^{n} + lambda d^{n}
238: where d^{n} = \hat{x}^{n} - x^{n} where \hat{x}^{n} is the solution returned from the inner solver.
240: The update, especially without inner nonlinear preconditioner, may be ill-scaled. If using the basic
241: linesearch, one may have to scale the update with -snes_linesearch_damping
243: This uses no derivative information thus will be much slower then Newton's method obtained with -snes_type ls
245: Only supports left non-linear preconditioning.
247: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESNGMRES, SNESQN, SNESNCG
248: M*/
249: PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes)
250: {
251: PetscErrorCode ierr;
252: SNES_NRichardson *neP;
253: SNESLineSearch linesearch;
256: snes->ops->destroy = SNESDestroy_NRichardson;
257: snes->ops->setup = SNESSetUp_NRichardson;
258: snes->ops->setfromoptions = SNESSetFromOptions_NRichardson;
259: snes->ops->view = SNESView_NRichardson;
260: snes->ops->solve = SNESSolve_NRichardson;
261: snes->ops->reset = SNESReset_NRichardson;
263: snes->usesksp = PETSC_FALSE;
264: snes->usesnpc = PETSC_TRUE;
266: snes->npcside= PC_LEFT;
268: SNESGetLineSearch(snes, &linesearch);
269: if (!((PetscObject)linesearch)->type_name) {
270: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
271: }
273: snes->alwayscomputesfinalresidual = PETSC_TRUE;
275: PetscNewLog(snes,&neP);
276: snes->data = (void*) neP;
278: if (!snes->tolerancesset) {
279: snes->max_funcs = 30000;
280: snes->max_its = 10000;
281: snes->stol = 1e-20;
282: }
283: return(0);
284: }