Actual source code: snesgs.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
  1:  #include <../src/snes/impls/gs/gsimpl.h>

  3: /*@
  4:    SNESNGSSetTolerances - Sets various parameters used in convergence tests.

  6:    Logically Collective on SNES

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


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

 22:    Level: intermediate

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

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


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

 54: /*@
 55:    SNESNGSGetTolerances - Gets various parameters used in convergence tests.

 57:    Not Collective

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

 67:    Notes:
 68:    The user can specify NULL for any parameter that is not needed.

 70:    Level: intermediate

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

 74: .seealso: SNESSetTolerances()
 75: @*/
 76: PetscErrorCode  SNESNGSGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit)
 77: {
 78:   SNES_NGS *gs = (SNES_NGS*)snes->data;

 82:   if (atol) *atol = gs->abstol;
 83:   if (rtol) *rtol = gs->rtol;
 84:   if (stol) *stol = gs->stol;
 85:   if (maxit) *maxit = gs->max_its;
 86:   return(0);
 87: }

 89: /*@
 90:    SNESNGSSetSweeps - Sets the number of sweeps of GS to use.

 92:    Input Parameters:
 93: +  snes   - the SNES context
 94: -  sweeps  - the number of sweeps of GS to perform.

 96:    Level: intermediate

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

100: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetPC(), SNESNGSGetSweeps()
101: @*/

103: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
104: {
105:   SNES_NGS *gs = (SNES_NGS*)snes->data;

109:   gs->sweeps = sweeps;
110:   return(0);
111: }

113: /*@
114:    SNESNGSGetSweeps - Gets the number of sweeps GS will use.

116:    Input Parameters:
117: .  snes   - the SNES context

119:    Output Parameters:
120: .  sweeps  - the number of sweeps of GS to perform.

122:    Level: intermediate

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

126: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetPC(), SNESNGSSetSweeps()
127: @*/
128: PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps)
129: {
130:   SNES_NGS *gs = (SNES_NGS*)snes->data;

134:   *sweeps = gs->sweeps;
135:   return(0);
136: }

138: PetscErrorCode SNESReset_NGS(SNES snes)
139: {
140:   SNES_NGS       *gs = (SNES_NGS*)snes->data;

144:   ISColoringDestroy(&gs->coloring);
145:   return(0);
146: }

148: PetscErrorCode SNESDestroy_NGS(SNES snes)
149: {

153:   SNESReset_NGS(snes);
154:   PetscFree(snes->data);
155:   return(0);
156: }

158: PetscErrorCode SNESSetUp_NGS(SNES snes)
159: {
161:   PetscErrorCode (*f)(SNES,Vec,Vec,void*);

164:   SNESGetNGS(snes,&f,NULL);
165:   if (!f) {
166:     SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
167:   }
168:   return(0);
169: }

171: PetscErrorCode SNESSetFromOptions_NGS(PetscOptionItems *PetscOptionsObject,SNES snes)
172: {
173:   SNES_NGS       *gs = (SNES_NGS*)snes->data;
175:   PetscInt       sweeps,max_its=PETSC_DEFAULT;
176:   PetscReal      rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT;
177:   PetscBool      flg,flg1,flg2,flg3;

180:   PetscOptionsHead(PetscOptionsObject,"SNES GS options");
181:   /* GS Options */
182:   PetscOptionsInt("-snes_ngs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",gs->sweeps,&sweeps,&flg);
183:   if (flg) {
184:     SNESNGSSetSweeps(snes,sweeps);
185:   }
186:   PetscOptionsReal("-snes_ngs_atol","Absolute residual tolerance for GS iteration","SNESComputeGS",gs->abstol,&atol,&flg);
187:   PetscOptionsReal("-snes_ngs_rtol","Relative residual tolerance for GS iteration","SNESComputeGS",gs->rtol,&rtol,&flg1);
188:   PetscOptionsReal("-snes_ngs_stol","Absolute update tolerance for GS iteration","SNESComputeGS",gs->stol,&stol,&flg2);
189:   PetscOptionsInt("-snes_ngs_max_it","Maximum number of sweeps of GS to apply","SNESComputeGS",gs->max_its,&max_its,&flg3);
190:   if (flg || flg1 || flg2 || flg3) {
191:     SNESNGSSetTolerances(snes,atol,rtol,stol,max_its);
192:   }
193:   flg  = PETSC_FALSE;
194:   PetscOptionsBool("-snes_ngs_secant","Use finite difference secant approximation with coloring","",flg,&flg,NULL);
195:   if (flg) {
196:     SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
197:     PetscInfo(snes,"Setting default finite difference secant approximation with coloring\n");
198:   }
199:   PetscOptionsReal("-snes_ngs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);
200:   PetscOptionsBool("-snes_ngs_secant_mat_coloring","Use the graph coloring of the Jacobian for the secant GS","",gs->secant_mat,&gs->secant_mat,&flg);

202:   PetscOptionsTail();
203:   return(0);
204: }

206: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
207: {
209:   PetscErrorCode (*f)(SNES,Vec,Vec,void*);
210:   SNES_NGS       *gs = (SNES_NGS*)snes->data;
211:   PetscBool      iascii;

214:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
215:   if (iascii) {
216:     DMSNESGetNGS(snes->dm,&f,NULL);
217:     if (f == SNESComputeNGSDefaultSecant) {
218:       PetscViewerASCIIPrintf(viewer,"  Use finite difference secant approximation with coloring with h = %g \n",(double)gs->h);
219:     }
220:   }
221:   return(0);
222: }

224: PetscErrorCode SNESSolve_NGS(SNES snes)
225: {
226:   Vec              F;
227:   Vec              X;
228:   Vec              B;
229:   PetscInt         i;
230:   PetscReal        fnorm;
231:   PetscErrorCode   ierr;
232:   SNESNormSchedule normschedule;


236:   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);

238:   PetscCitationsRegister(SNESCitation,&SNEScite);
239:   X = snes->vec_sol;
240:   F = snes->vec_func;
241:   B = snes->vec_rhs;

243:   PetscObjectSAWsTakeAccess((PetscObject)snes);
244:   snes->iter   = 0;
245:   snes->norm   = 0.;
246:   PetscObjectSAWsGrantAccess((PetscObject)snes);
247:   snes->reason = SNES_CONVERGED_ITERATING;

249:   SNESGetNormSchedule(snes, &normschedule);
250:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
251:     /* compute the initial function and preconditioned update delX */
252:     if (!snes->vec_func_init_set) {
253:       SNESComputeFunction(snes,X,F);
254:     } else snes->vec_func_init_set = PETSC_FALSE;

256:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
257:     SNESCheckFunctionNorm(snes,fnorm);
258:     PetscObjectSAWsTakeAccess((PetscObject)snes);
259:     snes->iter = 0;
260:     snes->norm = fnorm;
261:     PetscObjectSAWsGrantAccess((PetscObject)snes);
262:     SNESLogConvergenceHistory(snes,snes->norm,0);
263:     SNESMonitor(snes,0,snes->norm);

265:     /* test convergence */
266:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
267:     if (snes->reason) return(0);
268:   } else {
269:     PetscObjectSAWsGrantAccess((PetscObject)snes);
270:     SNESLogConvergenceHistory(snes,snes->norm,0);
271:   }

273:   /* Call general purpose update function */
274:   if (snes->ops->update) {
275:     (*snes->ops->update)(snes, snes->iter);
276:   }

278:   for (i = 0; i < snes->max_its; i++) {
279:     SNESComputeNGS(snes, B, X);
280:     /* only compute norms if requested or about to exit due to maximum iterations */
281:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
282:       SNESComputeFunction(snes,X,F);
283:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
284:       SNESCheckFunctionNorm(snes,fnorm);
285:       /* Monitor convergence */
286:       PetscObjectSAWsTakeAccess((PetscObject)snes);
287:       snes->iter = i+1;
288:       snes->norm = fnorm;
289:       PetscObjectSAWsGrantAccess((PetscObject)snes);
290:       SNESLogConvergenceHistory(snes,snes->norm,0);
291:       SNESMonitor(snes,snes->iter,snes->norm);
292:     }
293:     /* Test for convergence */
294:     if (normschedule == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
295:     if (snes->reason) return(0);
296:     /* Call general purpose update function */
297:     if (snes->ops->update) {
298:       (*snes->ops->update)(snes, snes->iter);
299:     }
300:   }
301:   if (normschedule == SNES_NORM_ALWAYS) {
302:     if (i == snes->max_its) {
303:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
304:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
305:     }
306:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
307:   return(0);
308: }

310: /*MC
311:   SNESNGS - Either calls the user-provided solution routine provided with SNESSetNGS() or does a finite difference secant approximation
312:             using coloring.

314:    Level: advanced

316:   Options Database:
317: +   -snes_ngs_sweeps - Number of sweeps of GS to apply
318: .    -snes_ngs_atol - Absolute residual tolerance for GS iteration
319: .    -snes_ngs_rtol - Relative residual tolerance for GS iteration
320: .    -snes_ngs_stol - Absolute update tolerance for GS iteration
321: .    -snes_ngs_max_it - Maximum number of sweeps of GS to apply
322: .    -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine, this is
323:                         used by default if no user provided Gauss-Seidel routine is available. Requires either that a DM that can compute a coloring
324:                         is available or a Jacobian sparse matrix is provided (from which to get the coloring).
325: .    -snes_ngs_secant_h - Differencing parameter for secant approximation
326: -    -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a DM is available.


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

333:    References:
334: .  1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
335:    SIAM Review, 57(4), 2015

337: .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types), SNESNGSSetSweeps(), SNESNGSSetTolerances()
338: M*/

340: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
341: {
342:   SNES_NGS        *gs;

346:   snes->ops->destroy        = SNESDestroy_NGS;
347:   snes->ops->setup          = SNESSetUp_NGS;
348:   snes->ops->setfromoptions = SNESSetFromOptions_NGS;
349:   snes->ops->view           = SNESView_NGS;
350:   snes->ops->solve          = SNESSolve_NGS;
351:   snes->ops->reset          = SNESReset_NGS;

353:   snes->usesksp = PETSC_FALSE;
354:   snes->usesnpc = PETSC_FALSE;

356:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

358:   if (!snes->tolerancesset) {
359:     snes->max_its   = 10000;
360:     snes->max_funcs = 10000;
361:   }

363:   PetscNewLog(snes,&gs);

365:   gs->sweeps  = 1;
366:   gs->rtol    = 1e-5;
367:   gs->abstol  = 1e-15;
368:   gs->stol    = 1e-12;
369:   gs->max_its = 50;
370:   gs->h       = 1e-8;

372:   snes->data = (void*) gs;
373:   return(0);
374: }