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