Actual source code: snesgs.c
petsc-3.12.5 2020-03-29
1: #include <../src/snes/impls/gs/gsimpl.h>
3: /*@
4: SNESNGSSetTolerances - Sets various parameters used in convergence tests.
6: Logically Collective on SNES
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
16: Options Database Keys:
17: + -snes_ngs_atol <abstol> - Sets abstol
18: . -snes_ngs_rtol <rtol> - Sets rtol
19: . -snes_ngs_stol <stol> - Sets stol
20: - -snes_max_it <maxit> - Sets maxit
22: Level: intermediate
24: .seealso: SNESSetTrustRegionTolerance()
25: @*/
26: PetscErrorCode SNESNGSSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit)
27: {
28: SNES_NGS *gs = (SNES_NGS*)snes->data;
33: if (abstol != PETSC_DEFAULT) {
34: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
35: gs->abstol = abstol;
36: }
37: if (rtol != PETSC_DEFAULT) {
38: if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol);
39: gs->rtol = rtol;
40: }
41: if (stol != PETSC_DEFAULT) {
42: if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol);
43: gs->stol = stol;
44: }
45: if (maxit != PETSC_DEFAULT) {
46: if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
47: gs->max_its = maxit;
48: }
49: return(0);
50: }
52: /*@
53: SNESNGSGetTolerances - Gets various parameters used in convergence tests.
55: Not Collective
57: Input Parameters:
58: + snes - the SNES context
59: . atol - absolute convergence tolerance
60: . rtol - relative convergence tolerance
61: . stol - convergence tolerance in terms of the norm
62: of the change in the solution between steps
63: - maxit - maximum number of iterations
65: Notes:
66: The user can specify NULL for any parameter that is not needed.
68: Level: intermediate
70: .seealso: SNESSetTolerances()
71: @*/
72: PetscErrorCode SNESNGSGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit)
73: {
74: SNES_NGS *gs = (SNES_NGS*)snes->data;
78: if (atol) *atol = gs->abstol;
79: if (rtol) *rtol = gs->rtol;
80: if (stol) *stol = gs->stol;
81: if (maxit) *maxit = gs->max_its;
82: return(0);
83: }
85: /*@
86: SNESNGSSetSweeps - Sets the number of sweeps of GS to use.
88: Input Parameters:
89: + snes - the SNES context
90: - sweeps - the number of sweeps of GS to perform.
92: Level: intermediate
94: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSGetSweeps()
95: @*/
97: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
98: {
99: SNES_NGS *gs = (SNES_NGS*)snes->data;
103: gs->sweeps = sweeps;
104: return(0);
105: }
107: /*@
108: SNESNGSGetSweeps - Gets the number of sweeps GS will use.
110: Input Parameters:
111: . snes - the SNES context
113: Output Parameters:
114: . sweeps - the number of sweeps of GS to perform.
116: Level: intermediate
118: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSSetSweeps()
119: @*/
120: PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps)
121: {
122: SNES_NGS *gs = (SNES_NGS*)snes->data;
126: *sweeps = gs->sweeps;
127: return(0);
128: }
130: PetscErrorCode SNESReset_NGS(SNES snes)
131: {
132: SNES_NGS *gs = (SNES_NGS*)snes->data;
136: ISColoringDestroy(&gs->coloring);
137: return(0);
138: }
140: PetscErrorCode SNESDestroy_NGS(SNES snes)
141: {
145: SNESReset_NGS(snes);
146: PetscFree(snes->data);
147: return(0);
148: }
150: PetscErrorCode SNESSetUp_NGS(SNES snes)
151: {
153: PetscErrorCode (*f)(SNES,Vec,Vec,void*);
156: SNESGetNGS(snes,&f,NULL);
157: if (!f) {
158: SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
159: }
160: return(0);
161: }
163: PetscErrorCode SNESSetFromOptions_NGS(PetscOptionItems *PetscOptionsObject,SNES snes)
164: {
165: SNES_NGS *gs = (SNES_NGS*)snes->data;
167: PetscInt sweeps,max_its=PETSC_DEFAULT;
168: PetscReal rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT;
169: PetscBool flg,flg1,flg2,flg3;
172: PetscOptionsHead(PetscOptionsObject,"SNES GS options");
173: /* GS Options */
174: PetscOptionsInt("-snes_ngs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",gs->sweeps,&sweeps,&flg);
175: if (flg) {
176: SNESNGSSetSweeps(snes,sweeps);
177: }
178: PetscOptionsReal("-snes_ngs_atol","Absolute residual tolerance for GS iteration","SNESComputeGS",gs->abstol,&atol,&flg);
179: PetscOptionsReal("-snes_ngs_rtol","Relative residual tolerance for GS iteration","SNESComputeGS",gs->rtol,&rtol,&flg1);
180: PetscOptionsReal("-snes_ngs_stol","Absolute update tolerance for GS iteration","SNESComputeGS",gs->stol,&stol,&flg2);
181: PetscOptionsInt("-snes_ngs_max_it","Maximum number of sweeps of GS to apply","SNESComputeGS",gs->max_its,&max_its,&flg3);
182: if (flg || flg1 || flg2 || flg3) {
183: SNESNGSSetTolerances(snes,atol,rtol,stol,max_its);
184: }
185: flg = PETSC_FALSE;
186: PetscOptionsBool("-snes_ngs_secant","Use finite difference secant approximation with coloring","",flg,&flg,NULL);
187: if (flg) {
188: SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
189: PetscInfo(snes,"Setting default finite difference secant approximation with coloring\n");
190: }
191: PetscOptionsReal("-snes_ngs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);
192: 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: PetscOptionsTail();
195: return(0);
196: }
198: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
199: {
201: PetscErrorCode (*f)(SNES,Vec,Vec,void*);
202: SNES_NGS *gs = (SNES_NGS*)snes->data;
203: PetscBool iascii;
206: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
207: if (iascii) {
208: DMSNESGetNGS(snes->dm,&f,NULL);
209: if (f == SNESComputeNGSDefaultSecant) {
210: PetscViewerASCIIPrintf(viewer," Use finite difference secant approximation with coloring with h = %g \n",(double)gs->h);
211: }
212: }
213: return(0);
214: }
216: PetscErrorCode SNESSolve_NGS(SNES snes)
217: {
218: Vec F;
219: Vec X;
220: Vec B;
221: PetscInt i;
222: PetscReal fnorm;
223: PetscErrorCode ierr;
224: SNESNormSchedule normschedule;
228: if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
230: PetscCitationsRegister(SNESCitation,&SNEScite);
231: X = snes->vec_sol;
232: F = snes->vec_func;
233: B = snes->vec_rhs;
235: PetscObjectSAWsTakeAccess((PetscObject)snes);
236: snes->iter = 0;
237: snes->norm = 0.;
238: PetscObjectSAWsGrantAccess((PetscObject)snes);
239: snes->reason = SNES_CONVERGED_ITERATING;
241: SNESGetNormSchedule(snes, &normschedule);
242: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
243: /* compute the initial function and preconditioned update delX */
244: if (!snes->vec_func_init_set) {
245: SNESComputeFunction(snes,X,F);
246: } else snes->vec_func_init_set = PETSC_FALSE;
248: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
249: SNESCheckFunctionNorm(snes,fnorm);
250: PetscObjectSAWsTakeAccess((PetscObject)snes);
251: snes->iter = 0;
252: snes->norm = fnorm;
253: PetscObjectSAWsGrantAccess((PetscObject)snes);
254: SNESLogConvergenceHistory(snes,snes->norm,0);
255: SNESMonitor(snes,0,snes->norm);
257: /* test convergence */
258: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
259: if (snes->reason) return(0);
260: } else {
261: PetscObjectSAWsGrantAccess((PetscObject)snes);
262: SNESLogConvergenceHistory(snes,snes->norm,0);
263: }
265: /* Call general purpose update function */
266: if (snes->ops->update) {
267: (*snes->ops->update)(snes, snes->iter);
268: }
270: for (i = 0; i < snes->max_its; i++) {
271: SNESComputeNGS(snes, B, X);
272: /* only compute norms if requested or about to exit due to maximum iterations */
273: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
274: SNESComputeFunction(snes,X,F);
275: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
276: SNESCheckFunctionNorm(snes,fnorm);
277: /* Monitor convergence */
278: PetscObjectSAWsTakeAccess((PetscObject)snes);
279: snes->iter = i+1;
280: snes->norm = fnorm;
281: PetscObjectSAWsGrantAccess((PetscObject)snes);
282: SNESLogConvergenceHistory(snes,snes->norm,0);
283: SNESMonitor(snes,snes->iter,snes->norm);
284: }
285: /* Test for convergence */
286: if (normschedule == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
287: if (snes->reason) return(0);
288: /* Call general purpose update function */
289: if (snes->ops->update) {
290: (*snes->ops->update)(snes, snes->iter);
291: }
292: }
293: if (normschedule == SNES_NORM_ALWAYS) {
294: if (i == snes->max_its) {
295: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
296: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
297: }
298: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
299: return(0);
300: }
302: /*MC
303: SNESNGS - Either calls the user-provided solution routine provided with SNESSetNGS() or does a finite difference secant approximation
304: using coloring.
306: Level: advanced
308: Options Database:
309: + -snes_ngs_sweeps <n> - Number of sweeps of GS to apply
310: . -snes_ngs_atol <atol> - Absolute residual tolerance for GS iteration
311: . -snes_ngs_rtol <rtol> - Relative residual tolerance for GS iteration
312: . -snes_ngs_stol <stol> - Absolute update tolerance for GS iteration
313: . -snes_ngs_max_it <maxit> - Maximum number of sweeps of GS to apply
314: . -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine, this is
315: used by default if no user provided Gauss-Seidel routine is available. Requires either that a DM that can compute a coloring
316: is available or a Jacobian sparse matrix is provided (from which to get the coloring).
317: . -snes_ngs_secant_h <h> - Differencing parameter for secant approximation
318: . -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a DM is available.
319: - -snes_norm_schedule <none, always, initialonly, finalonly, initalfinalonly> - how often the residual norms are computed
321: Notes:
322: the Gauss-Seidel smoother is inherited through composition. If a solver has been created with SNESGetNPC(), it will have
323: its parent's Gauss-Seidel routine associated with it.
325: By default this routine computes the solution norm at each iteration, this can be time consuming, you can turn this off with SNESSetNormSchedule() or -snes_norm_schedule
326: References:
327: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
328: SIAM Review, 57(4), 2015
330: .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types), SNESNGSSetSweeps(), SNESNGSSetTolerances(),
331: SNESSetNormSchedule()
332: M*/
334: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
335: {
336: SNES_NGS *gs;
340: snes->ops->destroy = SNESDestroy_NGS;
341: snes->ops->setup = SNESSetUp_NGS;
342: snes->ops->setfromoptions = SNESSetFromOptions_NGS;
343: snes->ops->view = SNESView_NGS;
344: snes->ops->solve = SNESSolve_NGS;
345: snes->ops->reset = SNESReset_NGS;
347: snes->usesksp = PETSC_FALSE;
348: snes->usesnpc = PETSC_FALSE;
350: snes->alwayscomputesfinalresidual = PETSC_FALSE;
352: if (!snes->tolerancesset) {
353: snes->max_its = 10000;
354: snes->max_funcs = 10000;
355: }
357: PetscNewLog(snes,&gs);
359: gs->sweeps = 1;
360: gs->rtol = 1e-5;
361: gs->abstol = PETSC_MACHINE_EPSILON;
362: gs->stol = 1000*PETSC_MACHINE_EPSILON;
363: gs->max_its = 50;
364: gs->h = PETSC_SQRT_MACHINE_EPSILON;
366: snes->data = (void*) gs;
367: return(0);
368: }