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;

 29:   PetscFunctionBegin;

 32:   if (abstol != (PetscReal)PETSC_DEFAULT) {
 33:     PetscCheck(abstol >= 0.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)abstol);
 34:     gs->abstol = abstol;
 35:   }
 36:   if (rtol != (PetscReal)PETSC_DEFAULT) {
 37:     PetscCheck(rtol >= 0.0 && rtol < 1.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative and less than 1.0", (double)rtol);
 38:     gs->rtol = rtol;
 39:   }
 40:   if (stol != (PetscReal)PETSC_DEFAULT) {
 41:     PetscCheck(stol >= 0.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Step tolerance %g must be non-negative", (double)stol);
 42:     gs->stol = stol;
 43:   }
 44:   if (maxit != PETSC_DEFAULT) {
 45:     PetscCheck(maxit >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of iterations %" PetscInt_FMT " must be non-negative", maxit);
 46:     gs->max_its = maxit;
 47:   }
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

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

 54:    Not Collective

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

 64:    Level: intermediate

 66:    Note:
 67:    The user can specify NULL for any parameter that is not needed.

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

 75:   PetscFunctionBegin;
 77:   if (atol) *atol = gs->abstol;
 78:   if (rtol) *rtol = gs->rtol;
 79:   if (stol) *stol = gs->stol;
 80:   if (maxit) *maxit = gs->max_its;
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

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

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

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

 94:    Level: intermediate

 96: .seealso: `SNESNCG`, `SNESSetNGS()`, `SNESGetNGS()`, `SNESSetNPC()`, `SNESNGSGetSweeps()`
 97: @*/

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

103:   PetscFunctionBegin;
105:   gs->sweeps = sweeps;
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

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

112:    Input Parameter:
113: .  snes   - the `SNES` context

115:    Output Parameter:
116: .  sweeps  - the number of sweeps of nonlinear GS to perform.

118:    Level: intermediate

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

126:   PetscFunctionBegin;
128:   *sweeps = gs->sweeps;
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }

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

136:   PetscFunctionBegin;
137:   PetscCall(ISColoringDestroy(&gs->coloring));
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: PetscErrorCode SNESDestroy_NGS(SNES snes)
142: {
143:   PetscFunctionBegin;
144:   PetscCall(SNESReset_NGS(snes));
145:   PetscCall(PetscFree(snes->data));
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }

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

153:   PetscFunctionBegin;
154:   PetscCall(SNESGetNGS(snes, &f, NULL));
155:   if (!f) PetscCall(SNESSetNGS(snes, SNESComputeNGSDefaultSecant, NULL));
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

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

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

185:   PetscOptionsHeadEnd();
186:   PetscFunctionReturn(PETSC_SUCCESS);
187: }

189: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
190: {
191:   PetscErrorCode (*f)(SNES, Vec, Vec, void *);
192:   SNES_NGS *gs = (SNES_NGS *)snes->data;
193:   PetscBool iascii;

195:   PetscFunctionBegin;
196:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
197:   if (iascii) {
198:     PetscCall(DMSNESGetNGS(snes->dm, &f, NULL));
199:     if (f == SNESComputeNGSDefaultSecant) PetscCall(PetscViewerASCIIPrintf(viewer, "  Use finite difference secant approximation with coloring with h = %g \n", (double)gs->h));
200:   }
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: PetscErrorCode SNESSolve_NGS(SNES snes)
205: {
206:   Vec              F;
207:   Vec              X;
208:   Vec              B;
209:   PetscInt         i;
210:   PetscReal        fnorm;
211:   SNESNormSchedule normschedule;

213:   PetscFunctionBegin;

215:   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

217:   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
218:   X = snes->vec_sol;
219:   F = snes->vec_func;
220:   B = snes->vec_rhs;

222:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
223:   snes->iter = 0;
224:   snes->norm = 0.;
225:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
226:   snes->reason = SNES_CONVERGED_ITERATING;

228:   PetscCall(SNESGetNormSchedule(snes, &normschedule));
229:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
230:     /* compute the initial function and preconditioned update delX */
231:     if (!snes->vec_func_init_set) {
232:       PetscCall(SNESComputeFunction(snes, X, F));
233:     } else snes->vec_func_init_set = PETSC_FALSE;

235:     PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
236:     SNESCheckFunctionNorm(snes, fnorm);
237:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
238:     snes->iter = 0;
239:     snes->norm = fnorm;
240:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
241:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
242:     PetscCall(SNESMonitor(snes, 0, snes->norm));

244:     /* test convergence */
245:     PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
246:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
247:   } else {
248:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
249:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
250:   }

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

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

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

289:    Level: advanced

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

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

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

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

315: .seealso: `SNESNCG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESSetNGS()`, `SNESType`, `SNESNGSSetSweeps()`, `SNESNGSSetTolerances()`,
316:           `SNESSetNormSchedule()`
317: M*/

319: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
320: {
321:   SNES_NGS *gs;

323:   PetscFunctionBegin;
324:   snes->ops->destroy        = SNESDestroy_NGS;
325:   snes->ops->setup          = SNESSetUp_NGS;
326:   snes->ops->setfromoptions = SNESSetFromOptions_NGS;
327:   snes->ops->view           = SNESView_NGS;
328:   snes->ops->solve          = SNESSolve_NGS;
329:   snes->ops->reset          = SNESReset_NGS;

331:   snes->usesksp = PETSC_FALSE;
332:   snes->usesnpc = PETSC_FALSE;

334:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

336:   if (!snes->tolerancesset) {
337:     snes->max_its   = 10000;
338:     snes->max_funcs = 10000;
339:   }

341:   PetscCall(PetscNew(&gs));

343:   gs->sweeps  = 1;
344:   gs->rtol    = 1e-5;
345:   gs->abstol  = PETSC_MACHINE_EPSILON;
346:   gs->stol    = 1000 * PETSC_MACHINE_EPSILON;
347:   gs->max_its = 50;
348:   gs->h       = PETSC_SQRT_MACHINE_EPSILON;

350:   snes->data = (void *)gs;
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }