Actual source code: snesrichardson.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  1: #include <../src/snes/impls/richardson/snesrichardsonimpl.h>


  6: PetscErrorCode SNESReset_NRichardson(SNES snes)
  7: {
  9:   return(0);
 10: }

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

 15:   Input Parameter:
 16: . snes - the SNES context

 18:   Application Interface Routine: SNESDestroy()
 19: */
 22: PetscErrorCode SNESDestroy_NRichardson(SNES snes)
 23: {

 27:   SNESReset_NRichardson(snes);
 28:   PetscFree(snes->data);
 29:   return(0);
 30: }

 32: /*
 33:    SNESSetUp_NRichardson - Sets up the internal data structures for the later use
 34:    of the SNESNRICHARDSON nonlinear solver.

 36:    Input Parameters:
 37: +  snes - the SNES context
 38: -  x - the solution vector

 40:    Application Interface Routine: SNESSetUp()
 41:  */
 44: PetscErrorCode SNESSetUp_NRichardson(SNES snes)
 45: {
 47:   if (snes->pcside == PC_RIGHT) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"NRichardson only supports left preconditioning");}
 48:   if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
 49:   return(0);
 50: }

 52: /*
 53:   SNESSetFromOptions_NRichardson - Sets various parameters for the SNESNEWTONLS method.

 55:   Input Parameter:
 56: . snes - the SNES context

 58:   Application Interface Routine: SNESSetFromOptions()
 59: */
 62: static PetscErrorCode SNESSetFromOptions_NRichardson(SNES snes)
 63: {
 65:   SNESLineSearch linesearch;

 68:   PetscOptionsHead("SNES Richardson options");
 69:   PetscOptionsTail();
 70:   if (!snes->linesearch) {
 71:     SNESGetLineSearch(snes, &linesearch);
 72:     SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
 73:   }
 74:   return(0);
 75: }

 77: /*
 78:   SNESView_NRichardson - Prints info from the SNESRichardson data structure.

 80:   Input Parameters:
 81: + SNES - the SNES context
 82: - viewer - visualization context

 84:   Application Interface Routine: SNESView()
 85: */
 88: static PetscErrorCode SNESView_NRichardson(SNES snes, PetscViewer viewer)
 89: {
 90:   PetscBool      iascii;

 94:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 95:   if (iascii) {
 96:   }
 97:   return(0);
 98: }

100: /*
101:   SNESSolve_NRichardson - Solves a nonlinear system with the Richardson method.

103:   Input Parameters:
104: . snes - the SNES context

106:   Output Parameter:
107: . outits - number of iterations until termination

109:   Application Interface Routine: SNESSolve()
110: */
113: PetscErrorCode SNESSolve_NRichardson(SNES snes)
114: {
115:   Vec                 X, Y, F;
116:   PetscReal           xnorm, fnorm, ynorm;
117:   PetscInt            maxits, i;
118:   PetscErrorCode      ierr;
119:   PetscBool           lsSuccess;
120:   SNESConvergedReason reason;

123:   snes->reason = SNES_CONVERGED_ITERATING;

125:   maxits = snes->max_its;        /* maximum number of iterations */
126:   X      = snes->vec_sol;        /* X^n */
127:   Y      = snes->vec_sol_update; /* \tilde X */
128:   F      = snes->vec_func;       /* residual vector */

130:   PetscObjectSAWsTakeAccess((PetscObject)snes);
131:   snes->iter = 0;
132:   snes->norm = 0.;
133:   PetscObjectSAWsGrantAccess((PetscObject)snes);

135:   if (snes->pc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
136:     SNESApplyNPC(snes,X,NULL,F);
137:     SNESGetConvergedReason(snes->pc,&reason);
138:     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
139:       snes->reason = SNES_DIVERGED_INNER;
140:       return(0);
141:     }
142:     VecNorm(F,NORM_2,&fnorm);
143:   } else {
144:     if (!snes->vec_func_init_set) {
145:       SNESComputeFunction(snes,X,F);
146:       if (snes->domainerror) {
147:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
148:         return(0);
149:       }
150:     } else snes->vec_func_init_set = PETSC_FALSE;

152:     VecNorm(F,NORM_2,&fnorm);
153:     if (PetscIsInfOrNanReal(fnorm)) {
154:       snes->reason = SNES_DIVERGED_FNORM_NAN;
155:       return(0);
156:     }
157:   }
158:   if (snes->pc && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
159:       SNESApplyNPC(snes,X,F,Y);
160:       SNESGetConvergedReason(snes->pc,&reason);
161:       if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
162:         snes->reason = SNES_DIVERGED_INNER;
163:         return(0);
164:       }
165:   } else {
166:     VecCopy(F,Y);
167:   }

169:   PetscObjectSAWsTakeAccess((PetscObject)snes);
170:   snes->norm = fnorm;
171:   PetscObjectSAWsGrantAccess((PetscObject)snes);
172:   SNESLogConvergenceHistory(snes,fnorm,0);
173:   SNESMonitor(snes,0,fnorm);

175:   /* test convergence */
176:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
177:   if (snes->reason) return(0);

179:   /* Call general purpose update function */
180:   if (snes->ops->update) {
181:     (*snes->ops->update)(snes, snes->iter);
182:   }

184:   /* set parameter for default relative tolerance convergence test */
185:   snes->ttol = fnorm*snes->rtol;
186:   /* test convergence */
187:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
188:   if (snes->reason) return(0);

190:   for (i = 1; i < maxits+1; i++) {
191:     lsSuccess = PETSC_TRUE;

193:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
194:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
195:     SNESLineSearchGetSuccess(snes->linesearch, &lsSuccess);
196:     if (!lsSuccess) {
197:       if (++snes->numFailures >= snes->maxFailures) {
198:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
199:         break;
200:       }
201:     }
202:     if (snes->nfuncs >= snes->max_funcs) {
203:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
204:       break;
205:     }
206:     if (snes->domainerror) {
207:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
208:       return(0);
209:     }

211:     /* Monitor convergence */
212:     PetscObjectSAWsTakeAccess((PetscObject)snes);
213:     snes->iter = i;
214:     snes->norm = fnorm;
215:     PetscObjectSAWsGrantAccess((PetscObject)snes);
216:     SNESLogConvergenceHistory(snes,snes->norm,0);
217:     SNESMonitor(snes,snes->iter,snes->norm);
218:     /* Test for convergence */
219:     (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
220:     if (snes->reason) break;

222:     /* Call general purpose update function */
223:     if (snes->ops->update) {
224:       (*snes->ops->update)(snes, snes->iter);
225:     }

227:     if (snes->pc) {
228:       if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
229:         SNESApplyNPC(snes,X,NULL,Y);
230:         VecNorm(F,NORM_2,&fnorm);
231:         VecCopy(Y,F);
232:       } else {
233:         SNESApplyNPC(snes,X,F,Y);
234:       }
235:       SNESGetConvergedReason(snes->pc,&reason);
236:       if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
237:         snes->reason = SNES_DIVERGED_INNER;
238:         return(0);
239:       }
240:     } else {
241:       VecCopy(F,Y);
242:     }
243:   }
244:   if (i == maxits+1) {
245:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
246:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
247:   }
248:   return(0);
249: }

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

254:   Level: beginner

256:   Options Database:
257: +   -snes_linesearch_type <l2,cp,basic> Line search type.
258: -   -snes_linesearch_damping<1.0> Damping for the line search.

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

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

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

271: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESNGMRES, SNESQN, SNESNCG
272: M*/
275: PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes)
276: {
277:   PetscErrorCode   ierr;
278:   SNES_NRichardson *neP;

281:   snes->ops->destroy        = SNESDestroy_NRichardson;
282:   snes->ops->setup          = SNESSetUp_NRichardson;
283:   snes->ops->setfromoptions = SNESSetFromOptions_NRichardson;
284:   snes->ops->view           = SNESView_NRichardson;
285:   snes->ops->solve          = SNESSolve_NRichardson;
286:   snes->ops->reset          = SNESReset_NRichardson;

288:   snes->usesksp = PETSC_FALSE;
289:   snes->usespc  = PETSC_TRUE;

291:   snes->pcside = PC_LEFT;

293:   PetscNewLog(snes,&neP);
294:   snes->data = (void*) neP;

296:   if (!snes->tolerancesset) {
297:     snes->max_funcs = 30000;
298:     snes->max_its   = 10000;
299:     snes->stol      = 1e-20;
300:   }
301:   return(0);
302: }