Actual source code: snesgs.c

petsc-3.5.0 2014-06-30
Report Typos and Errors
  1: #include <../src/snes/impls/gs/gsimpl.h>      /*I "petscsnes.h"  I*/

  5: /*@
  6:    SNESNGSSetTolerances - Sets various parameters used in convergence tests.

  8:    Logically Collective on SNES

 10:    Input Parameters:
 11: +  snes - the SNES context
 12: .  abstol - absolute convergence tolerance
 13: .  rtol - relative convergence tolerance
 14: .  stol -  convergence tolerance in terms of the norm of the change in the solution between steps,  || delta x || < stol*|| x ||
 15: -  maxit - maximum number of iterations


 18:    Options Database Keys:
 19: +    -snes_ngs_atol <abstol> - Sets abstol
 20: .    -snes_ngs_rtol <rtol> - Sets rtol
 21: .    -snes_ngs_stol <stol> - Sets stol
 22: -    -snes_max_it <maxit> - Sets maxit

 24:    Level: intermediate

 26: .keywords: SNES, nonlinear, gauss-seidel, set, convergence, tolerances

 28: .seealso: SNESSetTrustRegionTolerance()
 29: @*/
 30: PetscErrorCode  SNESNGSSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit)
 31: {
 32:   SNES_NGS *gs = (SNES_NGS*)snes->data;


 37:   if (abstol != PETSC_DEFAULT) {
 38:     if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
 39:     gs->abstol = abstol;
 40:   }
 41:   if (rtol != PETSC_DEFAULT) {
 42:     if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol);
 43:     gs->rtol = rtol;
 44:   }
 45:   if (stol != PETSC_DEFAULT) {
 46:     if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol);
 47:     gs->stol = stol;
 48:   }
 49:   if (maxit != PETSC_DEFAULT) {
 50:     if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
 51:     gs->max_its = maxit;
 52:   }
 53:   return(0);
 54: }

 58: /*@
 59:    SNESNGSGetTolerances - Gets various parameters used in convergence tests.

 61:    Not Collective

 63:    Input Parameters:
 64: +  snes - the SNES context
 65: .  atol - absolute convergence tolerance
 66: .  rtol - relative convergence tolerance
 67: .  stol -  convergence tolerance in terms of the norm
 68:            of the change in the solution between steps
 69: -  maxit - maximum number of iterations

 71:    Notes:
 72:    The user can specify NULL for any parameter that is not needed.

 74:    Level: intermediate

 76: .keywords: SNES, nonlinear, get, convergence, tolerances

 78: .seealso: SNESSetTolerances()
 79: @*/
 80: PetscErrorCode  SNESNGSGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit)
 81: {
 82:   SNES_NGS *gs = (SNES_NGS*)snes->data;

 86:   if (atol) *atol = gs->abstol;
 87:   if (rtol) *rtol = gs->rtol;
 88:   if (stol) *stol = gs->stol;
 89:   if (maxit) *maxit = gs->max_its;
 90:   return(0);
 91: }

 95: /*@
 96:    SNESNGSSetSweeps - Sets the number of sweeps of GS to use.

 98:    Input Parameters:
 99: +  snes   - the SNES context
100: -  sweeps  - the number of sweeps of GS to perform.

102:    Level: intermediate

104: .keywords: SNES, nonlinear, set, Gauss-Siedel

106: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetPC(), SNESNGSGetSweeps()
107: @*/

109: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
110: {
111:   SNES_NGS *gs = (SNES_NGS*)snes->data;

115:   gs->sweeps = sweeps;
116:   return(0);
117: }

121: /*@
122:    SNESNGSGetSweeps - Gets the number of sweeps GS will use.

124:    Input Parameters:
125: .  snes   - the SNES context

127:    Output Parameters:
128: .  sweeps  - the number of sweeps of GS to perform.

130:    Level: intermediate

132: .keywords: SNES, nonlinear, set, Gauss-Siedel

134: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetPC(), SNESNGSSetSweeps()
135: @*/
136: PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps)
137: {
138:   SNES_NGS *gs = (SNES_NGS*)snes->data;

142:   *sweeps = gs->sweeps;
143:   return(0);
144: }


149: PetscErrorCode SNESDefaultApplyNGS(SNES snes,Vec X,Vec F,void *ctx)
150: {
152:   /* see if there's a coloring on the DM */
153:   return(0);
154: }

158: PetscErrorCode SNESReset_NGS(SNES snes)
159: {
161:   return(0);
162: }

166: PetscErrorCode SNESDestroy_NGS(SNES snes)
167: {

171:   SNESReset_NGS(snes);
172:   PetscFree(snes->data);
173:   return(0);
174: }

178: PetscErrorCode SNESSetUp_NGS(SNES snes)
179: {
181:   PetscErrorCode (*f)(SNES,Vec,Vec,void*);

184:   SNESGetNGS(snes,&f,NULL);
185:   if (!f) {
186:     SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
187:   }
188:   return(0);
189: }

193: PetscErrorCode SNESSetFromOptions_NGS(SNES snes)
194: {
195:   SNES_NGS       *gs = (SNES_NGS*)snes->data;
197:   PetscInt       sweeps,max_its=PETSC_DEFAULT;
198:   PetscReal      rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT;
199:   PetscBool      flg,flg1,flg2,flg3;

202:   PetscOptionsHead("SNES GS options");
203:   /* GS Options */
204:   PetscOptionsInt("-snes_ngs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",gs->sweeps,&sweeps,&flg);
205:   if (flg) {
206:     SNESNGSSetSweeps(snes,sweeps);
207:   }
208:   PetscOptionsReal("-snes_ngs_atol","Number of sweeps of GS to apply","SNESComputeGS",gs->abstol,&atol,&flg);
209:   PetscOptionsReal("-snes_ngs_rtol","Number of sweeps of GS to apply","SNESComputeGS",gs->rtol,&rtol,&flg1);
210:   PetscOptionsReal("-snes_ngs_stol","Number of sweeps of GS to apply","SNESComputeGS",gs->stol,&stol,&flg2);
211:   PetscOptionsInt("-snes_ngs_max_it","Number of sweeps of GS to apply","SNESComputeGS",gs->max_its,&max_its,&flg3);
212:   if (flg || flg1 || flg2 || flg3) {
213:     SNESNGSSetTolerances(snes,atol,rtol,stol,max_its);
214:   }
215:   flg  = PETSC_FALSE;
216:   PetscOptionsBool("-snes_ngs_secant","Use pointwise secant local Jacobian approximation","",flg,&flg,NULL);
217:   if (flg) {
218:     SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
219:     PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");
220:   }
221:   PetscOptionsReal("-snes_ngs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);
222:   PetscOptionsBool("-snes_ngs_secant_mat_coloring","Use the Jacobian coloring for the secant GS","",gs->secant_mat,&gs->secant_mat,&flg);

224:   PetscOptionsTail();
225:   return(0);
226: }

230: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
231: {
233:   return(0);
234: }

238: PetscErrorCode SNESSolve_NGS(SNES snes)
239: {
240:   Vec              F;
241:   Vec              X;
242:   Vec              B;
243:   PetscInt         i;
244:   PetscReal        fnorm;
245:   PetscErrorCode   ierr;
246:   SNESNormSchedule normschedule;

249:   PetscCitationsRegister(SNESCitation,&SNEScite);
250:   X = snes->vec_sol;
251:   F = snes->vec_func;
252:   B = snes->vec_rhs;

254:   PetscObjectSAWsTakeAccess((PetscObject)snes);
255:   snes->iter   = 0;
256:   snes->norm   = 0.;
257:   PetscObjectSAWsGrantAccess((PetscObject)snes);
258:   snes->reason = SNES_CONVERGED_ITERATING;

260:   SNESGetNormSchedule(snes, &normschedule);
261:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
262:     /* compute the initial function and preconditioned update delX */
263:     if (!snes->vec_func_init_set) {
264:       SNESComputeFunction(snes,X,F);
265:       if (snes->domainerror) {
266:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
267:         return(0);
268:       }
269:     } else snes->vec_func_init_set = PETSC_FALSE;

271:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
272:     if (PetscIsInfOrNanReal(fnorm)) {
273:       snes->reason = SNES_DIVERGED_FNORM_NAN;
274:       return(0);
275:     }
276:     PetscObjectSAWsTakeAccess((PetscObject)snes);
277:     snes->iter = 0;
278:     snes->norm = fnorm;
279:     PetscObjectSAWsGrantAccess((PetscObject)snes);
280:     SNESLogConvergenceHistory(snes,snes->norm,0);
281:     SNESMonitor(snes,0,snes->norm);

283:     /* test convergence */
284:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
285:     if (snes->reason) return(0);
286:   } else {
287:     PetscObjectSAWsGrantAccess((PetscObject)snes);
288:     SNESLogConvergenceHistory(snes,snes->norm,0);
289:     SNESMonitor(snes,0,snes->norm);
290:   }

292:   /* Call general purpose update function */
293:   if (snes->ops->update) {
294:     (*snes->ops->update)(snes, snes->iter);
295:   }

297:   for (i = 0; i < snes->max_its; i++) {
298:     SNESComputeNGS(snes, B, X);
299:     /* only compute norms if requested or about to exit due to maximum iterations */
300:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
301:       SNESComputeFunction(snes,X,F);
302:       if (snes->domainerror) {
303:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
304:         return(0);
305:       }
306:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
307:       if (PetscIsInfOrNanReal(fnorm)) {
308:         snes->reason = SNES_DIVERGED_FNORM_NAN;
309:         return(0);
310:       }
311:     }
312:     /* Monitor convergence */
313:     PetscObjectSAWsTakeAccess((PetscObject)snes);
314:     snes->iter = i+1;
315:     snes->norm = fnorm;
316:     PetscObjectSAWsGrantAccess((PetscObject)snes);
317:     SNESLogConvergenceHistory(snes,snes->norm,0);
318:     SNESMonitor(snes,snes->iter,snes->norm);
319:     /* Test for convergence */
320:     if (normschedule == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
321:     if (snes->reason) return(0);
322:     /* Call general purpose update function */
323:     if (snes->ops->update) {
324:       (*snes->ops->update)(snes, snes->iter);
325:     }
326:   }
327:   if (normschedule == SNES_NORM_ALWAYS) {
328:     if (i == snes->max_its) {
329:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
330:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
331:     }
332:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
333:   return(0);
334: }

336: /*MC
337:   SNESNGS - Just calls the user-provided solution routine provided with SNESSetNGS()

339:    Level: advanced

341:   Notes:
342:   the Gauss-Seidel smoother is inherited through composition.  If a solver has been created with SNESGetPC(), it will have
343:   its parent's Gauss-Seidel routine associated with it.

345: .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types)
346: M*/

350: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
351: {
352:   SNES_NGS        *gs;

356:   snes->ops->destroy        = SNESDestroy_NGS;
357:   snes->ops->setup          = SNESSetUp_NGS;
358:   snes->ops->setfromoptions = SNESSetFromOptions_NGS;
359:   snes->ops->view           = SNESView_NGS;
360:   snes->ops->solve          = SNESSolve_NGS;
361:   snes->ops->reset          = SNESReset_NGS;

363:   snes->usesksp = PETSC_FALSE;
364:   snes->usespc  = PETSC_FALSE;

366:   if (!snes->tolerancesset) {
367:     snes->max_its   = 10000;
368:     snes->max_funcs = 10000;
369:   }

371:   PetscNewLog(snes,&gs);

373:   gs->sweeps  = 1;
374:   gs->rtol    = 1e-5;
375:   gs->abstol  = 1e-15;
376:   gs->stol    = 1e-12;
377:   gs->max_its = 50;
378:   gs->h       = 1e-8;

380:   snes->data = (void*) gs;
381:   return(0);
382: }