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(PetscOptionItems *PetscOptionsObject,SNES snes)
 49: {
 50:   PetscOptionsHead(PetscOptionsObject,"SNES Richardson options");
 51:   PetscOptionsTail();
 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:   }
 71:   return 0;
 72: }

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

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

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

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


 95:   snes->reason = SNES_CONVERGED_ITERATING;

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

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

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

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

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

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

144:   /* Call general purpose update function */
145:   if (snes->ops->update) {
146:     (*snes->ops->update)(snes, snes->iter);
147:   }

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

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

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

183:     /* Call general purpose update function */
184:     if (snes->ops->update) {
185:       (*snes->ops->update)(snes, snes->iter);
186:     }

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

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

215:   Level: beginner

217:   Options Database:
218: +   -snes_linesearch_type <l2,cp,basic> - Line search type.
219: -   -snes_linesearch_damping<1.0> - Damping for the line search.

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

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

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

233:      Only supports left non-linear preconditioning.

235: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESNGMRES, SNESQN, SNESNCG
236: M*/
237: PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes)
238: {
239:   SNES_NRichardson *neP;
240:   SNESLineSearch   linesearch;

242:   snes->ops->destroy        = SNESDestroy_NRichardson;
243:   snes->ops->setup          = SNESSetUp_NRichardson;
244:   snes->ops->setfromoptions = SNESSetFromOptions_NRichardson;
245:   snes->ops->view           = SNESView_NRichardson;
246:   snes->ops->solve          = SNESSolve_NRichardson;
247:   snes->ops->reset          = SNESReset_NRichardson;

249:   snes->usesksp = PETSC_FALSE;
250:   snes->usesnpc = PETSC_TRUE;

252:   snes->npcside= PC_LEFT;

254:   SNESGetLineSearch(snes, &linesearch);
255:   if (!((PetscObject)linesearch)->type_name) {
256:     SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
257:   }

259:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

261:   PetscNewLog(snes,&neP);
262:   snes->data = (void*) neP;

264:   if (!snes->tolerancesset) {
265:     snes->max_funcs = 30000;
266:     snes->max_its   = 10000;
267:     snes->stol      = 1e-20;
268:   }
269:   return 0;
270: }