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: Notes:
24: Use `PETSC_CURRENT` to retain the value for any parameter
26: All parameters must be non-negative
28: Developer Note:
29: Why can't the values set with `SNESSetTolerances()` be used?
31: .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetTolerances()`, `SNESSetTrustRegionTolerance()`
32: @*/
33: PetscErrorCode SNESNGSSetTolerances(SNES snes, PetscReal abstol, PetscReal rtol, PetscReal stol, PetscInt maxit)
34: {
35: SNES_NGS *gs = (SNES_NGS *)snes->data;
37: PetscFunctionBegin;
40: if (abstol != (PetscReal)PETSC_CURRENT) {
41: PetscCheck(abstol >= 0.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)abstol);
42: gs->abstol = abstol;
43: }
44: if (rtol != (PetscReal)PETSC_CURRENT) {
45: 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);
46: gs->rtol = rtol;
47: }
48: if (stol != (PetscReal)PETSC_CURRENT) {
49: PetscCheck(stol >= 0.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Step tolerance %g must be non-negative", (double)stol);
50: gs->stol = stol;
51: }
52: if (maxit != PETSC_CURRENT) {
53: PetscCheck(maxit >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of iterations %" PetscInt_FMT " must be non-negative", maxit);
54: gs->max_its = maxit;
55: }
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: /*@
60: SNESNGSGetTolerances - Gets various parameters used in convergence tests for nonlinear Gauss-Seidel `SNESNCG`
62: Not Collective
64: Input Parameters:
65: + snes - the `SNES` context
66: . atol - absolute convergence tolerance
67: . rtol - relative convergence tolerance
68: . stol - convergence tolerance in terms of the norm
69: of the change in the solution between steps
70: - maxit - maximum number of iterations
72: Level: intermediate
74: Note:
75: The user can specify `NULL` for any parameter that is not needed.
77: .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetTolerances()`
78: @*/
79: PetscErrorCode SNESNGSGetTolerances(SNES snes, PetscReal *atol, PetscReal *rtol, PetscReal *stol, PetscInt *maxit)
80: {
81: SNES_NGS *gs = (SNES_NGS *)snes->data;
83: PetscFunctionBegin;
85: if (atol) *atol = gs->abstol;
86: if (rtol) *rtol = gs->rtol;
87: if (stol) *stol = gs->stol;
88: if (maxit) *maxit = gs->max_its;
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: /*@
93: SNESNGSSetSweeps - Sets the number of sweeps of nonlinear GS to use in `SNESNCG`
95: Logically Collective
97: Input Parameters:
98: + snes - the `SNES` context
99: - sweeps - the number of sweeps of nonlinear GS to perform.
101: Options Database Key:
102: . -snes_ngs_sweeps <n> - Number of sweeps of nonlinear GS to apply
104: Level: intermediate
106: .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetNGS()`, `SNESGetNGS()`, `SNESSetNPC()`, `SNESNGSGetSweeps()`
107: @*/
108: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
109: {
110: SNES_NGS *gs = (SNES_NGS *)snes->data;
112: PetscFunctionBegin;
114: gs->sweeps = sweeps;
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /*@
119: SNESNGSGetSweeps - Gets the number of sweeps nonlinear GS will use in `SNESNCG`
121: Input Parameter:
122: . snes - the `SNES` context
124: Output Parameter:
125: . sweeps - the number of sweeps of nonlinear GS to perform.
127: Level: intermediate
129: .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetNGS()`, `SNESGetNGS()`, `SNESSetNPC()`, `SNESNGSSetSweeps()`
130: @*/
131: PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt *sweeps)
132: {
133: SNES_NGS *gs = (SNES_NGS *)snes->data;
135: PetscFunctionBegin;
137: *sweeps = gs->sweeps;
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode SNESReset_NGS(SNES snes)
142: {
143: SNES_NGS *gs = (SNES_NGS *)snes->data;
145: PetscFunctionBegin;
146: PetscCall(ISColoringDestroy(&gs->coloring));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode SNESDestroy_NGS(SNES snes)
151: {
152: PetscFunctionBegin;
153: PetscCall(SNESReset_NGS(snes));
154: PetscCall(PetscFree(snes->data));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: static PetscErrorCode SNESSetUp_NGS(SNES snes)
159: {
160: PetscErrorCode (*f)(SNES, Vec, Vec, void *);
162: PetscFunctionBegin;
163: PetscCall(SNESGetNGS(snes, &f, NULL));
164: if (!f) PetscCall(SNESSetNGS(snes, SNESComputeNGSDefaultSecant, NULL));
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: static PetscErrorCode SNESSetFromOptions_NGS(SNES snes, PetscOptionItems *PetscOptionsObject)
169: {
170: SNES_NGS *gs = (SNES_NGS *)snes->data;
171: PetscInt sweeps, max_its = PETSC_CURRENT;
172: PetscReal rtol = PETSC_CURRENT, atol = PETSC_CURRENT, stol = PETSC_CURRENT;
173: PetscBool flg, flg1, flg2, flg3;
175: PetscFunctionBegin;
176: PetscOptionsHeadBegin(PetscOptionsObject, "SNES GS options");
177: /* GS Options */
178: PetscCall(PetscOptionsInt("-snes_ngs_sweeps", "Number of sweeps of GS to apply", "SNESComputeGS", gs->sweeps, &sweeps, &flg));
179: if (flg) PetscCall(SNESNGSSetSweeps(snes, sweeps));
180: PetscCall(PetscOptionsReal("-snes_ngs_atol", "Absolute residual tolerance for GS iteration", "SNESComputeGS", gs->abstol, &atol, &flg));
181: PetscCall(PetscOptionsReal("-snes_ngs_rtol", "Relative residual tolerance for GS iteration", "SNESComputeGS", gs->rtol, &rtol, &flg1));
182: PetscCall(PetscOptionsReal("-snes_ngs_stol", "Absolute update tolerance for GS iteration", "SNESComputeGS", gs->stol, &stol, &flg2));
183: PetscCall(PetscOptionsInt("-snes_ngs_max_it", "Maximum number of sweeps of GS to apply", "SNESComputeGS", gs->max_its, &max_its, &flg3));
184: if (flg || flg1 || flg2 || flg3) PetscCall(SNESNGSSetTolerances(snes, atol, rtol, stol, max_its));
185: flg = PETSC_FALSE;
186: PetscCall(PetscOptionsBool("-snes_ngs_secant", "Use finite difference secant approximation with coloring", "", flg, &flg, NULL));
187: if (flg) {
188: PetscCall(SNESSetNGS(snes, SNESComputeNGSDefaultSecant, NULL));
189: PetscCall(PetscInfo(snes, "Setting default finite difference secant approximation with coloring\n"));
190: }
191: PetscCall(PetscOptionsReal("-snes_ngs_secant_h", "Differencing parameter for secant search", "", gs->h, &gs->h, NULL));
192: 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));
194: PetscOptionsHeadEnd();
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: static PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
199: {
200: PetscErrorCode (*f)(SNES, Vec, Vec, void *);
201: SNES_NGS *gs = (SNES_NGS *)snes->data;
202: PetscBool iascii;
204: PetscFunctionBegin;
205: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
206: if (iascii) {
207: PetscCall(DMSNESGetNGS(snes->dm, &f, NULL));
208: if (f == SNESComputeNGSDefaultSecant) PetscCall(PetscViewerASCIIPrintf(viewer, " Use finite difference secant approximation with coloring with h = %g \n", (double)gs->h));
209: }
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: static PetscErrorCode SNESSolve_NGS(SNES snes)
214: {
215: Vec F;
216: Vec X;
217: Vec B;
218: PetscInt i;
219: PetscReal fnorm;
220: SNESNormSchedule normschedule;
222: PetscFunctionBegin;
223: 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);
225: PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
226: X = snes->vec_sol;
227: F = snes->vec_func;
228: B = snes->vec_rhs;
230: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
231: snes->iter = 0;
232: snes->norm = 0.;
233: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
234: snes->reason = SNES_CONVERGED_ITERATING;
236: PetscCall(SNESGetNormSchedule(snes, &normschedule));
237: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
238: /* compute the initial function and preconditioned update delX */
239: if (!snes->vec_func_init_set) {
240: PetscCall(SNESComputeFunction(snes, X, F));
241: } else snes->vec_func_init_set = PETSC_FALSE;
243: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
244: SNESCheckFunctionNorm(snes, fnorm);
245: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
246: snes->iter = 0;
247: snes->norm = fnorm;
248: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
249: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
251: /* test convergence */
252: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
253: PetscCall(SNESMonitor(snes, 0, snes->norm));
254: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
255: } else {
256: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
257: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
258: }
260: /* Call general purpose update function */
261: PetscTryTypeMethod(snes, update, snes->iter);
263: for (i = 0; i < snes->max_its; i++) {
264: PetscCall(SNESComputeNGS(snes, B, X));
265: /* only compute norms if requested or about to exit due to maximum iterations */
266: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
267: PetscCall(SNESComputeFunction(snes, X, F));
268: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
269: SNESCheckFunctionNorm(snes, fnorm);
270: }
271: /* Monitor convergence */
272: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
273: snes->iter = i + 1;
274: snes->norm = fnorm;
275: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
276: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, snes->iter));
277: /* Test for convergence */
278: PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
279: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
280: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
281: /* Call general purpose update function */
282: PetscTryTypeMethod(snes, update, snes->iter);
283: }
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*MC
288: SNESNGS - Either calls the user-provided Gauss-Seidel solution routine provided with `SNESSetNGS()` or does a finite difference secant approximation
289: using coloring {cite}`bruneknepleysmithtu15`.
291: Level: advanced
293: Options Database Keys:
294: + -snes_ngs_sweeps <n> - Number of sweeps of nonlinear GS to apply
295: . -snes_ngs_atol <atol> - Absolute residual tolerance for nonlinear GS iteration
296: . -snes_ngs_rtol <rtol> - Relative residual tolerance for nonlinear GS iteration
297: . -snes_ngs_stol <stol> - Absolute update tolerance for nonlinear GS iteration
298: . -snes_ngs_max_it <maxit> - Maximum number of sweeps of nonlinea GS to apply
299: . -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine,
300: this is used by default if no user provided Gauss-Seidel routine is available.
301: Requires either that a `DM` that can compute a coloring
302: is available or a Jacobian sparse matrix is provided (from which to get the coloring).
303: . -snes_ngs_secant_h <h> - Differencing parameter for secant approximation
304: . -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a `DM` is available.
305: - -snes_norm_schedule <none, always, initialonly, finalonly, initialfinalonly> - how often the residual norms are computed
307: Notes:
308: the Gauss-Seidel smoother is inherited through composition. If a solver has been created with `SNESGetNPC()`, it will have
309: its parent's Gauss-Seidel routine associated with it.
311: By default this routine computes the solution norm at each iteration, this can be time consuming, you can turn this off with `SNESSetNormSchedule()`
312: or `-snes_norm_schedule none`
314: .seealso: [](ch_snes), `SNESNCG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESSetNGS()`, `SNESType`, `SNESNGSSetSweeps()`, `SNESNGSSetTolerances()`,
315: `SNESSetNormSchedule()`, `SNESNGSGetTolerances()`, `SNESNGSSetSweeps()`
316: M*/
318: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
319: {
320: SNES_NGS *gs;
322: PetscFunctionBegin;
323: snes->ops->destroy = SNESDestroy_NGS;
324: snes->ops->setup = SNESSetUp_NGS;
325: snes->ops->setfromoptions = SNESSetFromOptions_NGS;
326: snes->ops->view = SNESView_NGS;
327: snes->ops->solve = SNESSolve_NGS;
328: snes->ops->reset = SNESReset_NGS;
330: snes->usesksp = PETSC_FALSE;
331: snes->usesnpc = PETSC_FALSE;
332: snes->alwayscomputesfinalresidual = PETSC_FALSE;
334: PetscCall(SNESParametersInitialize(snes));
335: PetscObjectParameterSetDefault(snes, max_funcs, 10000);
336: PetscObjectParameterSetDefault(snes, max_its, 10000);
338: PetscCall(PetscNew(&gs));
339: gs->sweeps = 1;
340: gs->rtol = 1e-5;
341: gs->abstol = PETSC_MACHINE_EPSILON;
342: gs->stol = 1000 * PETSC_MACHINE_EPSILON;
343: gs->max_its = 50;
344: gs->h = PETSC_SQRT_MACHINE_EPSILON;
345: snes->data = (void *)gs;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }