Actual source code: snesgs.c
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
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: 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.
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: Notes:
64: The user can specify NULL for any parameter that is not needed.
66: Level: intermediate
68: .seealso: 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 GS to use.
85: Input Parameters:
86: + snes - the SNES context
87: - sweeps - the number of sweeps of GS to perform.
89: Level: intermediate
91: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSGetSweeps()
92: @*/
94: PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
95: {
96: SNES_NGS *gs = (SNES_NGS*)snes->data;
99: gs->sweeps = sweeps;
100: return 0;
101: }
103: /*@
104: SNESNGSGetSweeps - Gets the number of sweeps GS will use.
106: Input Parameters:
107: . snes - the SNES context
109: Output Parameters:
110: . sweeps - the number of sweeps of GS to perform.
112: Level: intermediate
114: .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSSetSweeps()
115: @*/
116: PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps)
117: {
118: SNES_NGS *gs = (SNES_NGS*)snes->data;
121: *sweeps = gs->sweeps;
122: return 0;
123: }
125: PetscErrorCode SNESReset_NGS(SNES snes)
126: {
127: SNES_NGS *gs = (SNES_NGS*)snes->data;
129: ISColoringDestroy(&gs->coloring);
130: return 0;
131: }
133: PetscErrorCode SNESDestroy_NGS(SNES snes)
134: {
135: SNESReset_NGS(snes);
136: PetscFree(snes->data);
137: return 0;
138: }
140: PetscErrorCode SNESSetUp_NGS(SNES snes)
141: {
142: PetscErrorCode (*f)(SNES,Vec,Vec,void*);
144: SNESGetNGS(snes,&f,NULL);
145: if (!f) {
146: SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
147: }
148: return 0;
149: }
151: PetscErrorCode SNESSetFromOptions_NGS(PetscOptionItems *PetscOptionsObject,SNES snes)
152: {
153: SNES_NGS *gs = (SNES_NGS*)snes->data;
154: PetscInt sweeps,max_its=PETSC_DEFAULT;
155: PetscReal rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT;
156: PetscBool flg,flg1,flg2,flg3;
158: PetscOptionsHead(PetscOptionsObject,"SNES GS options");
159: /* GS Options */
160: PetscOptionsInt("-snes_ngs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",gs->sweeps,&sweeps,&flg);
161: if (flg) {
162: SNESNGSSetSweeps(snes,sweeps);
163: }
164: PetscOptionsReal("-snes_ngs_atol","Absolute residual tolerance for GS iteration","SNESComputeGS",gs->abstol,&atol,&flg);
165: PetscOptionsReal("-snes_ngs_rtol","Relative residual tolerance for GS iteration","SNESComputeGS",gs->rtol,&rtol,&flg1);
166: PetscOptionsReal("-snes_ngs_stol","Absolute update tolerance for GS iteration","SNESComputeGS",gs->stol,&stol,&flg2);
167: PetscOptionsInt("-snes_ngs_max_it","Maximum number of sweeps of GS to apply","SNESComputeGS",gs->max_its,&max_its,&flg3);
168: if (flg || flg1 || flg2 || flg3) {
169: SNESNGSSetTolerances(snes,atol,rtol,stol,max_its);
170: }
171: flg = PETSC_FALSE;
172: PetscOptionsBool("-snes_ngs_secant","Use finite difference secant approximation with coloring","",flg,&flg,NULL);
173: if (flg) {
174: SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);
175: PetscInfo(snes,"Setting default finite difference secant approximation with coloring\n");
176: }
177: PetscOptionsReal("-snes_ngs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);
178: PetscOptionsBool("-snes_ngs_secant_mat_coloring","Use the graph coloring of the Jacobian for the secant GS","",gs->secant_mat,&gs->secant_mat,&flg);
180: PetscOptionsTail();
181: return 0;
182: }
184: PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
185: {
186: PetscErrorCode (*f)(SNES,Vec,Vec,void*);
187: SNES_NGS *gs = (SNES_NGS*)snes->data;
188: PetscBool iascii;
190: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
191: if (iascii) {
192: DMSNESGetNGS(snes->dm,&f,NULL);
193: if (f == SNESComputeNGSDefaultSecant) {
194: PetscViewerASCIIPrintf(viewer," Use finite difference secant approximation with coloring with h = %g \n",(double)gs->h);
195: }
196: }
197: return 0;
198: }
200: PetscErrorCode SNESSolve_NGS(SNES snes)
201: {
202: Vec F;
203: Vec X;
204: Vec B;
205: PetscInt i;
206: PetscReal fnorm;
207: SNESNormSchedule normschedule;
212: PetscCitationsRegister(SNESCitation,&SNEScite);
213: X = snes->vec_sol;
214: F = snes->vec_func;
215: B = snes->vec_rhs;
217: PetscObjectSAWsTakeAccess((PetscObject)snes);
218: snes->iter = 0;
219: snes->norm = 0.;
220: PetscObjectSAWsGrantAccess((PetscObject)snes);
221: snes->reason = SNES_CONVERGED_ITERATING;
223: SNESGetNormSchedule(snes, &normschedule);
224: if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
225: /* compute the initial function and preconditioned update delX */
226: if (!snes->vec_func_init_set) {
227: SNESComputeFunction(snes,X,F);
228: } else snes->vec_func_init_set = PETSC_FALSE;
230: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
231: SNESCheckFunctionNorm(snes,fnorm);
232: PetscObjectSAWsTakeAccess((PetscObject)snes);
233: snes->iter = 0;
234: snes->norm = fnorm;
235: PetscObjectSAWsGrantAccess((PetscObject)snes);
236: SNESLogConvergenceHistory(snes,snes->norm,0);
237: SNESMonitor(snes,0,snes->norm);
239: /* test convergence */
240: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
241: if (snes->reason) return 0;
242: } else {
243: PetscObjectSAWsGrantAccess((PetscObject)snes);
244: SNESLogConvergenceHistory(snes,snes->norm,0);
245: }
247: /* Call general purpose update function */
248: if (snes->ops->update) {
249: (*snes->ops->update)(snes, snes->iter);
250: }
252: for (i = 0; i < snes->max_its; i++) {
253: SNESComputeNGS(snes, B, X);
254: /* only compute norms if requested or about to exit due to maximum iterations */
255: if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
256: SNESComputeFunction(snes,X,F);
257: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
258: SNESCheckFunctionNorm(snes,fnorm);
259: /* Monitor convergence */
260: PetscObjectSAWsTakeAccess((PetscObject)snes);
261: snes->iter = i+1;
262: snes->norm = fnorm;
263: PetscObjectSAWsGrantAccess((PetscObject)snes);
264: SNESLogConvergenceHistory(snes,snes->norm,0);
265: SNESMonitor(snes,snes->iter,snes->norm);
266: }
267: /* Test for convergence */
268: if (normschedule == SNES_NORM_ALWAYS) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
269: if (snes->reason) return 0;
270: /* Call general purpose update function */
271: if (snes->ops->update) {
272: (*snes->ops->update)(snes, snes->iter);
273: }
274: }
275: if (normschedule == SNES_NORM_ALWAYS) {
276: if (i == snes->max_its) {
277: PetscInfo(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
278: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
279: }
280: } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
281: return 0;
282: }
284: /*MC
285: SNESNGS - Either calls the user-provided solution routine provided with SNESSetNGS() or does a finite difference secant approximation
286: using coloring.
288: Level: advanced
290: Options Database:
291: + -snes_ngs_sweeps <n> - Number of sweeps of GS to apply
292: . -snes_ngs_atol <atol> - Absolute residual tolerance for GS iteration
293: . -snes_ngs_rtol <rtol> - Relative residual tolerance for GS iteration
294: . -snes_ngs_stol <stol> - Absolute update tolerance for GS iteration
295: . -snes_ngs_max_it <maxit> - Maximum number of sweeps of GS to apply
296: . -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine, this is
297: used by default if no user provided Gauss-Seidel routine is available. Requires either that a DM that can compute a coloring
298: is available or a Jacobian sparse matrix is provided (from which to get the coloring).
299: . -snes_ngs_secant_h <h> - Differencing parameter for secant approximation
300: . -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a DM is available.
301: - -snes_norm_schedule <none, always, initialonly, finalonly, initialfinalonly> - how often the residual norms are computed
303: Notes:
304: the Gauss-Seidel smoother is inherited through composition. If a solver has been created with SNESGetNPC(), it will have
305: its parent's Gauss-Seidel routine associated with it.
307: 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
308: References:
309: . * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
310: SIAM Review, 57(4), 2015
312: .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types), SNESNGSSetSweeps(), SNESNGSSetTolerances(),
313: SNESSetNormSchedule()
314: M*/
316: PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
317: {
318: SNES_NGS *gs;
320: snes->ops->destroy = SNESDestroy_NGS;
321: snes->ops->setup = SNESSetUp_NGS;
322: snes->ops->setfromoptions = SNESSetFromOptions_NGS;
323: snes->ops->view = SNESView_NGS;
324: snes->ops->solve = SNESSolve_NGS;
325: snes->ops->reset = SNESReset_NGS;
327: snes->usesksp = PETSC_FALSE;
328: snes->usesnpc = PETSC_FALSE;
330: snes->alwayscomputesfinalresidual = PETSC_FALSE;
332: if (!snes->tolerancesset) {
333: snes->max_its = 10000;
334: snes->max_funcs = 10000;
335: }
337: PetscNewLog(snes,&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;
346: snes->data = (void*) gs;
347: return 0;
348: }