1 #include <../src/snes/impls/gs/gsimpl.h> /*I "petscsnes.h" I*/
2
3 /*@
4 SNESNGSSetTolerances - Sets various parameters used in convergence tests for nonlinear Gauss-Seidel `SNESNCG`
5
6 Logically Collective
7
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
14
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
20
21 Level: intermediate
22
23 Notes:
24 Use `PETSC_CURRENT` to retain the value for any parameter
25
26 All parameters must be non-negative
27
28 Developer Note:
29 Why can't the values set with `SNESSetTolerances()` be used?
30
31 .seealso: [](ch_snes), `SNES`, `SNESNCG`
32 @*/
SNESNGSSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit)33 PetscErrorCode SNESNGSSetTolerances(SNES snes, PetscReal abstol, PetscReal rtol, PetscReal stol, PetscInt maxit)
34 {
35 SNES_NGS *gs = (SNES_NGS *)snes->data;
36
37 PetscFunctionBegin;
38 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
39
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 }
58
59 /*@
60 SNESNGSGetTolerances - Gets various parameters used in convergence tests for nonlinear Gauss-Seidel `SNESNCG`
61
62 Not Collective
63
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
71
72 Level: intermediate
73
74 Note:
75 The user can specify `NULL` for any parameter that is not needed.
76
77 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetTolerances()`
78 @*/
SNESNGSGetTolerances(SNES snes,PetscReal * atol,PetscReal * rtol,PetscReal * stol,PetscInt * maxit)79 PetscErrorCode SNESNGSGetTolerances(SNES snes, PetscReal *atol, PetscReal *rtol, PetscReal *stol, PetscInt *maxit)
80 {
81 SNES_NGS *gs = (SNES_NGS *)snes->data;
82
83 PetscFunctionBegin;
84 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
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 }
91
92 /*@
93 SNESNGSSetSweeps - Sets the number of sweeps of nonlinear GS to use in `SNESNCG`
94
95 Logically Collective
96
97 Input Parameters:
98 + snes - the `SNES` context
99 - sweeps - the number of sweeps of nonlinear GS to perform.
100
101 Options Database Key:
102 . -snes_ngs_sweeps <n> - Number of sweeps of nonlinear GS to apply
103
104 Level: intermediate
105
106 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetNGS()`, `SNESGetNGS()`, `SNESSetNPC()`, `SNESNGSGetSweeps()`
107 @*/
SNESNGSSetSweeps(SNES snes,PetscInt sweeps)108 PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
109 {
110 SNES_NGS *gs = (SNES_NGS *)snes->data;
111
112 PetscFunctionBegin;
113 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
114 gs->sweeps = sweeps;
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 /*@
119 SNESNGSGetSweeps - Gets the number of sweeps nonlinear GS will use in `SNESNCG`
120
121 Input Parameter:
122 . snes - the `SNES` context
123
124 Output Parameter:
125 . sweeps - the number of sweeps of nonlinear GS to perform.
126
127 Level: intermediate
128
129 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESSetNGS()`, `SNESGetNGS()`, `SNESSetNPC()`, `SNESNGSSetSweeps()`
130 @*/
SNESNGSGetSweeps(SNES snes,PetscInt * sweeps)131 PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt *sweeps)
132 {
133 SNES_NGS *gs = (SNES_NGS *)snes->data;
134
135 PetscFunctionBegin;
136 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
137 *sweeps = gs->sweeps;
138 PetscFunctionReturn(PETSC_SUCCESS);
139 }
140
SNESReset_NGS(SNES snes)141 static PetscErrorCode SNESReset_NGS(SNES snes)
142 {
143 SNES_NGS *gs = (SNES_NGS *)snes->data;
144
145 PetscFunctionBegin;
146 PetscCall(ISColoringDestroy(&gs->coloring));
147 PetscFunctionReturn(PETSC_SUCCESS);
148 }
149
SNESDestroy_NGS(SNES snes)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 }
157
SNESSetUp_NGS(SNES snes)158 static PetscErrorCode SNESSetUp_NGS(SNES snes)
159 {
160 PetscErrorCode (*f)(SNES, Vec, Vec, void *);
161
162 PetscFunctionBegin;
163 PetscCall(SNESGetNGS(snes, &f, NULL));
164 if (!f) PetscCall(SNESSetNGS(snes, SNESComputeNGSDefaultSecant, NULL));
165 PetscFunctionReturn(PETSC_SUCCESS);
166 }
167
SNESSetFromOptions_NGS(SNES snes,PetscOptionItems PetscOptionsObject)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;
174
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));
193
194 PetscOptionsHeadEnd();
195 PetscFunctionReturn(PETSC_SUCCESS);
196 }
197
SNESView_NGS(SNES snes,PetscViewer viewer)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 isascii;
203
204 PetscFunctionBegin;
205 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
206 if (isascii) {
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 }
212
SNESSolve_NGS(SNES snes)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;
221
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);
224
225 PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
226 X = snes->vec_sol;
227 F = snes->vec_func;
228 B = snes->vec_rhs;
229
230 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
231 snes->iter = 0;
232 snes->norm = 0.;
233 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
234 snes->reason = SNES_CONVERGED_ITERATING;
235
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) PetscCall(SNESComputeFunction(snes, X, F));
240 else snes->vec_func_init_set = PETSC_FALSE;
241
242 PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
243 SNESCheckFunctionDomainError(snes, fnorm);
244 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
245 snes->iter = 0;
246 snes->norm = fnorm;
247 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
248 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
249
250 /* test convergence */
251 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
252 PetscCall(SNESMonitor(snes, 0, snes->norm));
253 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
254 } else {
255 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
256 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
257 }
258
259 /* Call general purpose update function */
260 PetscTryTypeMethod(snes, update, snes->iter);
261
262 for (i = 0; i < snes->max_its; i++) {
263 PetscCall(SNESComputeNGS(snes, B, X));
264 /* only compute norms if requested or about to exit due to maximum iterations */
265 if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
266 PetscCall(SNESComputeFunction(snes, X, F));
267 PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
268 SNESCheckFunctionDomainError(snes, fnorm);
269 }
270 /* Monitor convergence */
271 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
272 snes->iter = i + 1;
273 snes->norm = fnorm;
274 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
275 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, snes->iter));
276 /* Test for convergence */
277 PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
278 PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
279 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
280 /* Call general purpose update function */
281 PetscTryTypeMethod(snes, update, snes->iter);
282 }
283 PetscFunctionReturn(PETSC_SUCCESS);
284 }
285
286 /*MC
287 SNESNGS - Either calls the user-provided Gauss-Seidel solution routine provided with `SNESSetNGS()` or does a finite difference secant approximation
288 using coloring {cite}`bruneknepleysmithtu15`.
289
290 Level: advanced
291
292 Options Database Keys:
293 + -snes_ngs_sweeps <n> - Number of sweeps of nonlinear GS to apply
294 . -snes_ngs_atol <atol> - Absolute residual tolerance for nonlinear GS iteration
295 . -snes_ngs_rtol <rtol> - Relative residual tolerance for nonlinear GS iteration
296 . -snes_ngs_stol <stol> - Absolute update tolerance for nonlinear GS iteration
297 . -snes_ngs_max_it <maxit> - Maximum number of sweeps of nonlinea GS to apply
298 . -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine,
299 this is used by default if no user provided Gauss-Seidel routine is available.
300 Requires either that a `DM` that can compute a coloring
301 is available or a Jacobian sparse matrix is provided (from which to get the coloring).
302 . -snes_ngs_secant_h <h> - Differencing parameter for secant approximation
303 . -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a `DM` is available.
304 - -snes_norm_schedule <none, always, initialonly, finalonly, initialfinalonly> - how often the residual norms are computed
305
306 Notes:
307 the Gauss-Seidel smoother is inherited through composition. If a solver has been created with `SNESGetNPC()`, it will have
308 its parent's Gauss-Seidel routine associated with it.
309
310 By default this routine computes the solution norm at each iteration, this can be time consuming, you can turn this off with `SNESSetNormSchedule()`
311 or `-snes_norm_schedule none`
312
313 .seealso: [](ch_snes), `SNESNCG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESSetNGS()`, `SNESType`, `SNESNGSSetSweeps()`, `SNESNGSSetTolerances()`,
314 `SNESSetNormSchedule()`, `SNESNGSGetTolerances()`, `SNESNGSSetSweeps()`
315 M*/
316
SNESCreate_NGS(SNES snes)317 PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
318 {
319 SNES_NGS *gs;
320
321 PetscFunctionBegin;
322 snes->ops->destroy = SNESDestroy_NGS;
323 snes->ops->setup = SNESSetUp_NGS;
324 snes->ops->setfromoptions = SNESSetFromOptions_NGS;
325 snes->ops->view = SNESView_NGS;
326 snes->ops->solve = SNESSolve_NGS;
327 snes->ops->reset = SNESReset_NGS;
328
329 snes->usesksp = PETSC_FALSE;
330 snes->usesnpc = PETSC_FALSE;
331 snes->alwayscomputesfinalresidual = PETSC_FALSE;
332
333 PetscCall(SNESParametersInitialize(snes));
334 PetscObjectParameterSetDefault(snes, max_funcs, 10000);
335 PetscObjectParameterSetDefault(snes, max_its, 10000);
336
337 PetscCall(PetscNew(&gs));
338 gs->sweeps = 1;
339 gs->rtol = 1e-5;
340 gs->abstol = PETSC_MACHINE_EPSILON;
341 gs->stol = 1000 * PETSC_MACHINE_EPSILON;
342 gs->max_its = 50;
343 gs->h = PETSC_SQRT_MACHINE_EPSILON;
344 snes->data = (void *)gs;
345 PetscFunctionReturn(PETSC_SUCCESS);
346 }
347