Actual source code: snesgs.c

petsc-master 2019-11-16
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: .seealso: SNESSetTrustRegionTolerance()
 25: @*/
 26: PetscErrorCode  SNESNGSSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit)
 27: {
 28:   SNES_NGS *gs = (SNES_NGS*)snes->data;


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

 52: /*@
 53:    SNESNGSGetTolerances - Gets various parameters used in convergence tests.

 55:    Not Collective

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

 65:    Notes:
 66:    The user can specify NULL for any parameter that is not needed.

 68:    Level: intermediate

 70: .seealso: SNESSetTolerances()
 71: @*/
 72: PetscErrorCode  SNESNGSGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit)
 73: {
 74:   SNES_NGS *gs = (SNES_NGS*)snes->data;

 78:   if (atol) *atol = gs->abstol;
 79:   if (rtol) *rtol = gs->rtol;
 80:   if (stol) *stol = gs->stol;
 81:   if (maxit) *maxit = gs->max_its;
 82:   return(0);
 83: }

 85: /*@
 86:    SNESNGSSetSweeps - Sets the number of sweeps of GS to use.

 88:    Input Parameters:
 89: +  snes   - the SNES context
 90: -  sweeps  - the number of sweeps of GS to perform.

 92:    Level: intermediate

 94: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSGetSweeps()
 95: @*/

 97: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
 98: {
 99:   SNES_NGS *gs = (SNES_NGS*)snes->data;

103:   gs->sweeps = sweeps;
104:   return(0);
105: }

107: /*@
108:    SNESNGSGetSweeps - Gets the number of sweeps GS will use.

110:    Input Parameters:
111: .  snes   - the SNES context

113:    Output Parameters:
114: .  sweeps  - the number of sweeps of GS to perform.

116:    Level: intermediate

118: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSSetSweeps()
119: @*/
120: PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps)
121: {
122:   SNES_NGS *gs = (SNES_NGS*)snes->data;

126:   *sweeps = gs->sweeps;
127:   return(0);
128: }

130: PetscErrorCode SNESReset_NGS(SNES snes)
131: {
132:   SNES_NGS       *gs = (SNES_NGS*)snes->data;

136:   ISColoringDestroy(&gs->coloring);
137:   return(0);
138: }

140: PetscErrorCode SNESDestroy_NGS(SNES snes)
141: {

145:   SNESReset_NGS(snes);
146:   PetscFree(snes->data);
147:   return(0);
148: }

150: PetscErrorCode SNESSetUp_NGS(SNES snes)
151: {
153:   PetscErrorCode (*f)(SNES,Vec,Vec,void*);

156:   SNESGetNGS(snes,&f,NULL);
157:   if (!f) {
158:     SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
159:   }
160:   return(0);
161: }

163: PetscErrorCode SNESSetFromOptions_NGS(PetscOptionItems *PetscOptionsObject,SNES snes)
164: {
165:   SNES_NGS       *gs = (SNES_NGS*)snes->data;
167:   PetscInt       sweeps,max_its=PETSC_DEFAULT;
168:   PetscReal      rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT;
169:   PetscBool      flg,flg1,flg2,flg3;

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

194:   PetscOptionsTail();
195:   return(0);
196: }

198: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
199: {
201:   PetscErrorCode (*f)(SNES,Vec,Vec,void*);
202:   SNES_NGS       *gs = (SNES_NGS*)snes->data;
203:   PetscBool      iascii;

206:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
207:   if (iascii) {
208:     DMSNESGetNGS(snes->dm,&f,NULL);
209:     if (f == SNESComputeNGSDefaultSecant) {
210:       PetscViewerASCIIPrintf(viewer,"  Use finite difference secant approximation with coloring with h = %g \n",(double)gs->h);
211:     }
212:   }
213:   return(0);
214: }

216: PetscErrorCode SNESSolve_NGS(SNES snes)
217: {
218:   Vec              F;
219:   Vec              X;
220:   Vec              B;
221:   PetscInt         i;
222:   PetscReal        fnorm;
223:   PetscErrorCode   ierr;
224:   SNESNormSchedule normschedule;


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

230:   PetscCitationsRegister(SNESCitation,&SNEScite);
231:   X = snes->vec_sol;
232:   F = snes->vec_func;
233:   B = snes->vec_rhs;

235:   PetscObjectSAWsTakeAccess((PetscObject)snes);
236:   snes->iter   = 0;
237:   snes->norm   = 0.;
238:   PetscObjectSAWsGrantAccess((PetscObject)snes);
239:   snes->reason = SNES_CONVERGED_ITERATING;

241:   SNESGetNormSchedule(snes, &normschedule);
242:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
243:     /* compute the initial function and preconditioned update delX */
244:     if (!snes->vec_func_init_set) {
245:       SNESComputeFunction(snes,X,F);
246:     } else snes->vec_func_init_set = PETSC_FALSE;

248:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
249:     SNESCheckFunctionNorm(snes,fnorm);
250:     PetscObjectSAWsTakeAccess((PetscObject)snes);
251:     snes->iter = 0;
252:     snes->norm = fnorm;
253:     PetscObjectSAWsGrantAccess((PetscObject)snes);
254:     SNESLogConvergenceHistory(snes,snes->norm,0);
255:     SNESMonitor(snes,0,snes->norm);

257:     /* test convergence */
258:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
259:     if (snes->reason) return(0);
260:   } else {
261:     PetscObjectSAWsGrantAccess((PetscObject)snes);
262:     SNESLogConvergenceHistory(snes,snes->norm,0);
263:   }

265:   /* Call general purpose update function */
266:   if (snes->ops->update) {
267:     (*snes->ops->update)(snes, snes->iter);
268:   }

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

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

306:    Level: advanced

308:   Options Database:
309: +   -snes_ngs_sweeps <n> - Number of sweeps of GS to apply
310: .    -snes_ngs_atol <atol> - Absolute residual tolerance for GS iteration
311: .    -snes_ngs_rtol <rtol> - Relative residual tolerance for GS iteration
312: .    -snes_ngs_stol <stol> - Absolute update tolerance for GS iteration
313: .    -snes_ngs_max_it <maxit> - Maximum number of sweeps of GS to apply
314: .    -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine, this is
315:                         used by default if no user provided Gauss-Seidel routine is available. Requires either that a DM that can compute a coloring
316:                         is available or a Jacobian sparse matrix is provided (from which to get the coloring).
317: .    -snes_ngs_secant_h <h> - Differencing parameter for secant approximation
318: .    -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a DM is available.
319: -    -snes_norm_schedule <none, always, initialonly, finalonly, initalfinalonly> - how often the residual norms are computed

321:   Notes:
322:   the Gauss-Seidel smoother is inherited through composition.  If a solver has been created with SNESGetNPC(), it will have
323:   its parent's Gauss-Seidel routine associated with it.

325:   By default this routine computes the solution norm at each iteration, this can be time consuming, you can turn this off with SNESSetNormSchedule() or -snes_norm_schedule
326:    References:
327: .  1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
328:    SIAM Review, 57(4), 2015

330: .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types), SNESNGSSetSweeps(), SNESNGSSetTolerances(),
331:           SNESSetNormSchedule()
332: M*/

334: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
335: {
336:   SNES_NGS        *gs;

340:   snes->ops->destroy        = SNESDestroy_NGS;
341:   snes->ops->setup          = SNESSetUp_NGS;
342:   snes->ops->setfromoptions = SNESSetFromOptions_NGS;
343:   snes->ops->view           = SNESView_NGS;
344:   snes->ops->solve          = SNESSolve_NGS;
345:   snes->ops->reset          = SNESReset_NGS;

347:   snes->usesksp = PETSC_FALSE;
348:   snes->usesnpc = PETSC_FALSE;

350:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

352:   if (!snes->tolerancesset) {
353:     snes->max_its   = 10000;
354:     snes->max_funcs = 10000;
355:   }

357:   PetscNewLog(snes,&gs);

359:   gs->sweeps  = 1;
360:   gs->rtol    = 1e-5;
361:   gs->abstol  = PETSC_MACHINE_EPSILON;
362:   gs->stol    = 1000*PETSC_MACHINE_EPSILON;
363:   gs->max_its = 50;
364:   gs->h       = PETSC_SQRT_MACHINE_EPSILON;

366:   snes->data = (void*) gs;
367:   return(0);
368: }