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: [](ch_snes), `SNES`, `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: [](ch_snes), `SNES`, `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:   Logically Collective

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

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

 96:   Level: intermediate

 98: .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetNGS()`, `SNESGetNGS()`, `SNESSetNPC()`, `SNESNGSGetSweeps()`
 99: @*/
100: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
101: {
102:   SNES_NGS *gs = (SNES_NGS *)snes->data;

104:   PetscFunctionBegin;
106:   gs->sweeps = sweeps;
107:   PetscFunctionReturn(PETSC_SUCCESS);
108: }

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

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

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

119:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

214:   PetscFunctionBegin;

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

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

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

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

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

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

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

256:   for (i = 0; i < snes->max_its; i++) {
257:     PetscCall(SNESComputeNGS(snes, B, X));
258:     /* only compute norms if requested or about to exit due to maximum iterations */
259:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
260:       PetscCall(SNESComputeFunction(snes, X, F));
261:       PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
262:       SNESCheckFunctionNorm(snes, fnorm);
263:     }
264:     /* Monitor convergence */
265:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
266:     snes->iter = i + 1;
267:     snes->norm = fnorm;
268:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
269:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, snes->iter));
270:     /* Test for convergence */
271:     PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
272:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
273:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
274:     /* Call general purpose update function */
275:     PetscTryTypeMethod(snes, update, snes->iter);
276:   }
277:   PetscFunctionReturn(PETSC_SUCCESS);
278: }

280: /*MC
281:   SNESNGS - Either calls the user-provided Gauss-Seidel solution routine provided with `SNESSetNGS()` or does a finite difference secant approximation
282:             using coloring {cite}`bruneknepleysmithtu15`.

284:    Level: advanced

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

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

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

307: .seealso: [](ch_snes), `SNESNCG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESSetNGS()`, `SNESType`, `SNESNGSSetSweeps()`, `SNESNGSSetTolerances()`,
308:           `SNESSetNormSchedule()`, `SNESNGSGetTolerances()`, `SNESNGSSetSweeps()`
309: M*/

311: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
312: {
313:   SNES_NGS *gs;

315:   PetscFunctionBegin;
316:   snes->ops->destroy        = SNESDestroy_NGS;
317:   snes->ops->setup          = SNESSetUp_NGS;
318:   snes->ops->setfromoptions = SNESSetFromOptions_NGS;
319:   snes->ops->view           = SNESView_NGS;
320:   snes->ops->solve          = SNESSolve_NGS;
321:   snes->ops->reset          = SNESReset_NGS;

323:   snes->usesksp = PETSC_FALSE;
324:   snes->usesnpc = PETSC_FALSE;

326:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

328:   if (!snes->tolerancesset) {
329:     snes->max_its   = 10000;
330:     snes->max_funcs = 10000;
331:   }

333:   PetscCall(PetscNew(&gs));

335:   gs->sweeps  = 1;
336:   gs->rtol    = 1e-5;
337:   gs->abstol  = PETSC_MACHINE_EPSILON;
338:   gs->stol    = 1000 * PETSC_MACHINE_EPSILON;
339:   gs->max_its = 50;
340:   gs->h       = PETSC_SQRT_MACHINE_EPSILON;

342:   snes->data = (void *)gs;
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }