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: }