Actual source code: snesgs.c
petsc-3.11.4 2019-09-28
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: .keywords: SNES, nonlinear, gauss-seidel, set, convergence, tolerances
26: .seealso: SNESSetTrustRegionTolerance()
27: @*/
28: PetscErrorCode SNESNGSSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit)
29: {
30: SNES_NGS *gs = (SNES_NGS*)snes->data;
35: if (abstol != PETSC_DEFAULT) {
36: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
37: gs->abstol = abstol;
38: }
39: if (rtol != PETSC_DEFAULT) {
40: 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);
41: gs->rtol = rtol;
42: }
43: if (stol != PETSC_DEFAULT) {
44: if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol);
45: gs->stol = stol;
46: }
47: if (maxit != PETSC_DEFAULT) {
48: if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
49: gs->max_its = maxit;
50: }
51: return(0);
52: }
54: /*@
55: SNESNGSGetTolerances - Gets various parameters used in convergence tests.
57: Not Collective
59: Input Parameters:
60: + snes - the SNES context
61: . atol - absolute convergence tolerance
62: . rtol - relative convergence tolerance
63: . stol - convergence tolerance in terms of the norm
64: of the change in the solution between steps
65: - maxit - maximum number of iterations
67: Notes:
68: The user can specify NULL for any parameter that is not needed.
70: Level: intermediate
72: .keywords: SNES, nonlinear, get, convergence, tolerances
74: .seealso: SNESSetTolerances()
75: @*/
76: PetscErrorCode SNESNGSGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit)
77: {
78: SNES_NGS *gs = (SNES_NGS*)snes->data;
82: if (atol) *atol = gs->abstol;
83: if (rtol) *rtol = gs->rtol;
84: if (stol) *stol = gs->stol;
85: if (maxit) *maxit = gs->max_its;
86: return(0);
87: }
89: /*@
90: SNESNGSSetSweeps - Sets the number of sweeps of GS to use.
92: Input Parameters:
93: + snes - the SNES context
94: - sweeps - the number of sweeps of GS to perform.
96: Level: intermediate
98: .keywords: SNES, nonlinear, set, Gauss-Siedel
100: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSGetSweeps()
101: @*/
103: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
104: {
105: SNES_NGS *gs = (SNES_NGS*)snes->data;
109: gs->sweeps = sweeps;
110: return(0);
111: }
113: /*@
114: SNESNGSGetSweeps - Gets the number of sweeps GS will use.
116: Input Parameters:
117: . snes - the SNES context
119: Output Parameters:
120: . sweeps - the number of sweeps of GS to perform.
122: Level: intermediate
124: .keywords: SNES, nonlinear, set, Gauss-Siedel
126: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSSetSweeps()
127: @*/
128: PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps)
129: {
130: SNES_NGS *gs = (SNES_NGS*)snes->data;
134: *sweeps = gs->sweeps;
135: return(0);
136: }
138: PetscErrorCode SNESReset_NGS(SNES snes)
139: {
140: SNES_NGS *gs = (SNES_NGS*)snes->data;
144: ISColoringDestroy(&gs->coloring);
145: return(0);
146: }
148: PetscErrorCode SNESDestroy_NGS(SNES snes)
149: {
153: SNESReset_NGS(snes);
154: PetscFree(snes->data);
155: return(0);
156: }
158: PetscErrorCode SNESSetUp_NGS(SNES snes)
159: {
161: PetscErrorCode (*f)(SNES,Vec,Vec,void*);
164: SNESGetNGS(snes,&f,NULL);
165: if (!f) {
166: SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
167: }
168: return(0);
169: }
171: PetscErrorCode SNESSetFromOptions_NGS(PetscOptionItems *PetscOptionsObject,SNES snes)
172: {
173: SNES_NGS *gs = (SNES_NGS*)snes->data;
175: PetscInt sweeps,max_its=PETSC_DEFAULT;
176: PetscReal rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT;
177: PetscBool flg,flg1,flg2,flg3;
180: PetscOptionsHead(PetscOptionsObject,"SNES GS options");
181: /* GS Options */
182: PetscOptionsInt("-snes_ngs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",gs->sweeps,&sweeps,&flg);
183: if (flg) {
184: SNESNGSSetSweeps(snes,sweeps);
185: }
186: PetscOptionsReal("-snes_ngs_atol","Absolute residual tolerance for GS iteration","SNESComputeGS",gs->abstol,&atol,&flg);
187: PetscOptionsReal("-snes_ngs_rtol","Relative residual tolerance for GS iteration","SNESComputeGS",gs->rtol,&rtol,&flg1);
188: PetscOptionsReal("-snes_ngs_stol","Absolute update tolerance for GS iteration","SNESComputeGS",gs->stol,&stol,&flg2);
189: PetscOptionsInt("-snes_ngs_max_it","Maximum number of sweeps of GS to apply","SNESComputeGS",gs->max_its,&max_its,&flg3);
190: if (flg || flg1 || flg2 || flg3) {
191: SNESNGSSetTolerances(snes,atol,rtol,stol,max_its);
192: }
193: flg = PETSC_FALSE;
194: PetscOptionsBool("-snes_ngs_secant","Use finite difference secant approximation with coloring","",flg,&flg,NULL);
195: if (flg) {
196: SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
197: PetscInfo(snes,"Setting default finite difference secant approximation with coloring\n");
198: }
199: PetscOptionsReal("-snes_ngs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);
200: PetscOptionsBool("-snes_ngs_secant_mat_coloring","Use the graph coloring of the Jacobian for the secant GS","",gs->secant_mat,&gs->secant_mat,&flg);
202: PetscOptionsTail();
203: return(0);
204: }
206: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
207: {
209: PetscErrorCode (*f)(SNES,Vec,Vec,void*);
210: SNES_NGS *gs = (SNES_NGS*)snes->data;
211: PetscBool iascii;
214: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
215: if (iascii) {
216: DMSNESGetNGS(snes->dm,&f,NULL);
217: if (f == SNESComputeNGSDefaultSecant) {
218: PetscViewerASCIIPrintf(viewer," Use finite difference secant approximation with coloring with h = %g \n",(double)gs->h);
219: }
220: }
221: return(0);
222: }
224: PetscErrorCode SNESSolve_NGS(SNES snes)
225: {
226: Vec F;
227: Vec X;
228: Vec B;
229: PetscInt i;
230: PetscReal fnorm;
231: PetscErrorCode ierr;
232: SNESNormSchedule normschedule;
236: 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);
238: PetscCitationsRegister(SNESCitation,&SNEScite);
239: X = snes->vec_sol;
240: F = snes->vec_func;
241: B = snes->vec_rhs;
243: PetscObjectSAWsTakeAccess((PetscObject)snes);
244: snes->iter = 0;
245: snes->norm = 0.;
246: PetscObjectSAWsGrantAccess((PetscObject)snes);
247: snes->reason = SNES_CONVERGED_ITERATING;
249: SNESGetNormSchedule(snes, &normschedule);
250: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
251: /* compute the initial function and preconditioned update delX */
252: if (!snes->vec_func_init_set) {
253: SNESComputeFunction(snes,X,F);
254: } else snes->vec_func_init_set = PETSC_FALSE;
256: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
257: SNESCheckFunctionNorm(snes,fnorm);
258: PetscObjectSAWsTakeAccess((PetscObject)snes);
259: snes->iter = 0;
260: snes->norm = fnorm;
261: PetscObjectSAWsGrantAccess((PetscObject)snes);
262: SNESLogConvergenceHistory(snes,snes->norm,0);
263: SNESMonitor(snes,0,snes->norm);
265: /* test convergence */
266: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
267: if (snes->reason) return(0);
268: } else {
269: PetscObjectSAWsGrantAccess((PetscObject)snes);
270: SNESLogConvergenceHistory(snes,snes->norm,0);
271: }
273: /* Call general purpose update function */
274: if (snes->ops->update) {
275: (*snes->ops->update)(snes, snes->iter);
276: }
278: for (i = 0; i < snes->max_its; i++) {
279: SNESComputeNGS(snes, B, X);
280: /* only compute norms if requested or about to exit due to maximum iterations */
281: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
282: SNESComputeFunction(snes,X,F);
283: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
284: SNESCheckFunctionNorm(snes,fnorm);
285: /* Monitor convergence */
286: PetscObjectSAWsTakeAccess((PetscObject)snes);
287: snes->iter = i+1;
288: snes->norm = fnorm;
289: PetscObjectSAWsGrantAccess((PetscObject)snes);
290: SNESLogConvergenceHistory(snes,snes->norm,0);
291: SNESMonitor(snes,snes->iter,snes->norm);
292: }
293: /* Test for convergence */
294: if (normschedule == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
295: if (snes->reason) return(0);
296: /* Call general purpose update function */
297: if (snes->ops->update) {
298: (*snes->ops->update)(snes, snes->iter);
299: }
300: }
301: if (normschedule == SNES_NORM_ALWAYS) {
302: if (i == snes->max_its) {
303: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
304: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
305: }
306: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
307: return(0);
308: }
310: /*MC
311: SNESNGS - Either calls the user-provided solution routine provided with SNESSetNGS() or does a finite difference secant approximation
312: using coloring.
314: Level: advanced
316: Options Database:
317: + -snes_ngs_sweeps - Number of sweeps of GS to apply
318: . -snes_ngs_atol - Absolute residual tolerance for GS iteration
319: . -snes_ngs_rtol - Relative residual tolerance for GS iteration
320: . -snes_ngs_stol - Absolute update tolerance for GS iteration
321: . -snes_ngs_max_it - Maximum number of sweeps of GS to apply
322: . -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine, this is
323: used by default if no user provided Gauss-Seidel routine is available. Requires either that a DM that can compute a coloring
324: is available or a Jacobian sparse matrix is provided (from which to get the coloring).
325: . -snes_ngs_secant_h - Differencing parameter for secant approximation
326: - -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a DM is available.
329: Notes:
330: the Gauss-Seidel smoother is inherited through composition. If a solver has been created with SNESGetNPC(), it will have
331: its parent's Gauss-Seidel routine associated with it.
333: References:
334: . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
335: SIAM Review, 57(4), 2015
337: .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types), SNESNGSSetSweeps(), SNESNGSSetTolerances()
338: M*/
340: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
341: {
342: SNES_NGS *gs;
346: snes->ops->destroy = SNESDestroy_NGS;
347: snes->ops->setup = SNESSetUp_NGS;
348: snes->ops->setfromoptions = SNESSetFromOptions_NGS;
349: snes->ops->view = SNESView_NGS;
350: snes->ops->solve = SNESSolve_NGS;
351: snes->ops->reset = SNESReset_NGS;
353: snes->usesksp = PETSC_FALSE;
354: snes->usesnpc = PETSC_FALSE;
356: snes->alwayscomputesfinalresidual = PETSC_FALSE;
358: if (!snes->tolerancesset) {
359: snes->max_its = 10000;
360: snes->max_funcs = 10000;
361: }
363: PetscNewLog(snes,&gs);
365: gs->sweeps = 1;
366: gs->rtol = 1e-5;
367: gs->abstol = 1e-15;
368: gs->stol = 1e-12;
369: gs->max_its = 50;
370: gs->h = 1e-8;
372: snes->data = (void*) gs;
373: return(0);
374: }