Actual source code: snesgs.c

  1: #include <../src/snes/impls/gs/gsimpl.h>

  3: /*@
  4:    SNESNGSSetTolerances - Sets various parameters used in convergence tests for nonlinear Gauss-Seidel `SNESNCG`

  6:    Logically Collective

  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

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

 21:    Level: intermediate

 23: .seealso: `SNESNCG`, `SNESSetTrustRegionTolerance()`
 24: @*/
 25: PetscErrorCode SNESNGSSetTolerances(SNES snes, PetscReal abstol, PetscReal rtol, PetscReal stol, PetscInt maxit)
 26: {
 27:   SNES_NGS *gs = (SNES_NGS *)snes->data;


 31:   if (abstol != PETSC_DEFAULT) {
 33:     gs->abstol = abstol;
 34:   }
 35:   if (rtol != PETSC_DEFAULT) {
 37:     gs->rtol = rtol;
 38:   }
 39:   if (stol != PETSC_DEFAULT) {
 41:     gs->stol = stol;
 42:   }
 43:   if (maxit != PETSC_DEFAULT) {
 45:     gs->max_its = maxit;
 46:   }
 47:   return 0;
 48: }

 50: /*@
 51:    SNESNGSGetTolerances - Gets various parameters used in convergence tests for nonlinear Gauss-Seidel `SNESNCG`

 53:    Not Collective

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

 63:    Note:
 64:    The user can specify NULL for any parameter that is not needed.

 66:    Level: intermediate

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

 75:   if (atol) *atol = gs->abstol;
 76:   if (rtol) *rtol = gs->rtol;
 77:   if (stol) *stol = gs->stol;
 78:   if (maxit) *maxit = gs->max_its;
 79:   return 0;
 80: }

 82: /*@
 83:    SNESNGSSetSweeps - Sets the number of sweeps of nonlinear GS to use in `SNESNCG`

 85:    Input Parameters:
 86: +  snes   - the `SNES` context
 87: -  sweeps  - the number of sweeps of nonlinear GS to perform.

 89:   Options Database Key:
 90: .   -snes_ngs_sweeps <n> - Number of sweeps of nonlinear GS to apply

 92:    Level: intermediate

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

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

102:   gs->sweeps = sweeps;
103:   return 0;
104: }

106: /*@
107:    SNESNGSGetSweeps - Gets the number of sweeps nonlinear GS will use in `SNESNCG`

109:    Input Parameters:
110: .  snes   - the `SNES` context

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

115:    Level: intermediate

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

124:   *sweeps = gs->sweeps;
125:   return 0;
126: }

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

132:   ISColoringDestroy(&gs->coloring);
133:   return 0;
134: }

136: PetscErrorCode SNESDestroy_NGS(SNES snes)
137: {
138:   SNESReset_NGS(snes);
139:   PetscFree(snes->data);
140:   return 0;
141: }

143: PetscErrorCode SNESSetUp_NGS(SNES snes)
144: {
145:   PetscErrorCode (*f)(SNES, Vec, Vec, void *);

147:   SNESGetNGS(snes, &f, NULL);
148:   if (!f) SNESSetNGS(snes, SNESComputeNGSDefaultSecant, NULL);
149:   return 0;
150: }

152: PetscErrorCode SNESSetFromOptions_NGS(SNES snes, PetscOptionItems *PetscOptionsObject)
153: {
154:   SNES_NGS *gs = (SNES_NGS *)snes->data;
155:   PetscInt  sweeps, max_its = PETSC_DEFAULT;
156:   PetscReal rtol = PETSC_DEFAULT, atol = PETSC_DEFAULT, stol = PETSC_DEFAULT;
157:   PetscBool flg, flg1, flg2, flg3;

159:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES GS options");
160:   /* GS Options */
161:   PetscOptionsInt("-snes_ngs_sweeps", "Number of sweeps of GS to apply", "SNESComputeGS", gs->sweeps, &sweeps, &flg);
162:   if (flg) SNESNGSSetSweeps(snes, sweeps);
163:   PetscOptionsReal("-snes_ngs_atol", "Absolute residual tolerance for GS iteration", "SNESComputeGS", gs->abstol, &atol, &flg);
164:   PetscOptionsReal("-snes_ngs_rtol", "Relative residual tolerance for GS iteration", "SNESComputeGS", gs->rtol, &rtol, &flg1);
165:   PetscOptionsReal("-snes_ngs_stol", "Absolute update tolerance for GS iteration", "SNESComputeGS", gs->stol, &stol, &flg2);
166:   PetscOptionsInt("-snes_ngs_max_it", "Maximum number of sweeps of GS to apply", "SNESComputeGS", gs->max_its, &max_its, &flg3);
167:   if (flg || flg1 || flg2 || flg3) SNESNGSSetTolerances(snes, atol, rtol, stol, max_its);
168:   flg = PETSC_FALSE;
169:   PetscOptionsBool("-snes_ngs_secant", "Use finite difference secant approximation with coloring", "", flg, &flg, NULL);
170:   if (flg) {
171:     SNESSetNGS(snes, SNESComputeNGSDefaultSecant, NULL);
172:     PetscInfo(snes, "Setting default finite difference secant approximation with coloring\n");
173:   }
174:   PetscOptionsReal("-snes_ngs_secant_h", "Differencing parameter for secant search", "", gs->h, &gs->h, NULL);
175:   PetscOptionsBool("-snes_ngs_secant_mat_coloring", "Use the graph coloring of the Jacobian for the secant GS", "", gs->secant_mat, &gs->secant_mat, &flg);

177:   PetscOptionsHeadEnd();
178:   return 0;
179: }

181: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
182: {
183:   PetscErrorCode (*f)(SNES, Vec, Vec, void *);
184:   SNES_NGS *gs = (SNES_NGS *)snes->data;
185:   PetscBool iascii;

187:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
188:   if (iascii) {
189:     DMSNESGetNGS(snes->dm, &f, NULL);
190:     if (f == SNESComputeNGSDefaultSecant) PetscViewerASCIIPrintf(viewer, "  Use finite difference secant approximation with coloring with h = %g \n", (double)gs->h);
191:   }
192:   return 0;
193: }

195: PetscErrorCode SNESSolve_NGS(SNES snes)
196: {
197:   Vec              F;
198:   Vec              X;
199:   Vec              B;
200:   PetscInt         i;
201:   PetscReal        fnorm;
202:   SNESNormSchedule normschedule;



207:   PetscCitationsRegister(SNESCitation, &SNEScite);
208:   X = snes->vec_sol;
209:   F = snes->vec_func;
210:   B = snes->vec_rhs;

212:   PetscObjectSAWsTakeAccess((PetscObject)snes);
213:   snes->iter = 0;
214:   snes->norm = 0.;
215:   PetscObjectSAWsGrantAccess((PetscObject)snes);
216:   snes->reason = SNES_CONVERGED_ITERATING;

218:   SNESGetNormSchedule(snes, &normschedule);
219:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
220:     /* compute the initial function and preconditioned update delX */
221:     if (!snes->vec_func_init_set) {
222:       SNESComputeFunction(snes, X, F);
223:     } else snes->vec_func_init_set = PETSC_FALSE;

225:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
226:     SNESCheckFunctionNorm(snes, fnorm);
227:     PetscObjectSAWsTakeAccess((PetscObject)snes);
228:     snes->iter = 0;
229:     snes->norm = fnorm;
230:     PetscObjectSAWsGrantAccess((PetscObject)snes);
231:     SNESLogConvergenceHistory(snes, snes->norm, 0);
232:     SNESMonitor(snes, 0, snes->norm);

234:     /* test convergence */
235:     PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
236:     if (snes->reason) return 0;
237:   } else {
238:     PetscObjectSAWsGrantAccess((PetscObject)snes);
239:     SNESLogConvergenceHistory(snes, snes->norm, 0);
240:   }

242:   /* Call general purpose update function */
243:   PetscTryTypeMethod(snes, update, snes->iter);

245:   for (i = 0; i < snes->max_its; i++) {
246:     SNESComputeNGS(snes, B, X);
247:     /* only compute norms if requested or about to exit due to maximum iterations */
248:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
249:       SNESComputeFunction(snes, X, F);
250:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
251:       SNESCheckFunctionNorm(snes, fnorm);
252:       /* Monitor convergence */
253:       PetscObjectSAWsTakeAccess((PetscObject)snes);
254:       snes->iter = i + 1;
255:       snes->norm = fnorm;
256:       PetscObjectSAWsGrantAccess((PetscObject)snes);
257:       SNESLogConvergenceHistory(snes, snes->norm, 0);
258:       SNESMonitor(snes, snes->iter, snes->norm);
259:     }
260:     /* Test for convergence */
261:     if (normschedule == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
262:     if (snes->reason) return 0;
263:     /* Call general purpose update function */
264:     PetscTryTypeMethod(snes, update, snes->iter);
265:   }
266:   if (normschedule == SNES_NORM_ALWAYS) {
267:     if (i == snes->max_its) {
268:       PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its);
269:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
270:     }
271:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
272:   return 0;
273: }

275: /*MC
276:   SNESNGS - Either calls the user-provided solution routine provided with `SNESSetNGS()` or does a finite difference secant approximation
277:             using coloring.

279:    Level: advanced

281:   Options Database Keys:
282: +   -snes_ngs_sweeps <n> - Number of sweeps of nonlinear GS to apply
283: .    -snes_ngs_atol <atol> - Absolute residual tolerance for nonlinear GS iteration
284: .    -snes_ngs_rtol <rtol> - Relative residual tolerance for nonlinear GS iteration
285: .    -snes_ngs_stol <stol> - Absolute update tolerance for nonlinear GS iteration
286: .    -snes_ngs_max_it <maxit> - Maximum number of sweeps of nonlinea GS to apply
287: .    -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine, this is
288:                         used by default if no user provided Gauss-Seidel routine is available. Requires either that a `DM` that can compute a coloring
289:                         is available or a Jacobian sparse matrix is provided (from which to get the coloring).
290: .    -snes_ngs_secant_h <h> - Differencing parameter for secant approximation
291: .    -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a DM is available.
292: -    -snes_norm_schedule <none, always, initialonly, finalonly, initialfinalonly> - how often the residual norms are computed

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

298:   By default this routine computes the solution norm at each iteration, this can be time consuming, you can turn this off with `SNESSetNormSchedule()`
299:   or -snes_norm_schedule none

301:    References:
302: .  * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
303:    SIAM Review, 57(4), 2015

305: .seealso: `SNESNCG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESSetNGS()`, `SNESType`, `SNESNGSSetSweeps()`, `SNESNGSSetTolerances()`,
306:           `SNESSetNormSchedule()`
307: M*/

309: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
310: {
311:   SNES_NGS *gs;

313:   snes->ops->destroy        = SNESDestroy_NGS;
314:   snes->ops->setup          = SNESSetUp_NGS;
315:   snes->ops->setfromoptions = SNESSetFromOptions_NGS;
316:   snes->ops->view           = SNESView_NGS;
317:   snes->ops->solve          = SNESSolve_NGS;
318:   snes->ops->reset          = SNESReset_NGS;

320:   snes->usesksp = PETSC_FALSE;
321:   snes->usesnpc = PETSC_FALSE;

323:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

325:   if (!snes->tolerancesset) {
326:     snes->max_its   = 10000;
327:     snes->max_funcs = 10000;
328:   }

330:   PetscNew(&gs);

332:   gs->sweeps  = 1;
333:   gs->rtol    = 1e-5;
334:   gs->abstol  = PETSC_MACHINE_EPSILON;
335:   gs->stol    = 1000 * PETSC_MACHINE_EPSILON;
336:   gs->max_its = 50;
337:   gs->h       = PETSC_SQRT_MACHINE_EPSILON;

339:   snes->data = (void *)gs;
340:   return 0;
341: }