xref: /petsc/src/snes/impls/gs/snesgs.c (revision 2b8d69ca7ea5fe9190df62c1dce3bbd66fce84dd)
1 #include <../src/snes/impls/gs/gsimpl.h>      /*I "petscsnes.h"  I*/
2 
3 #undef __FUNCT__
4 #define __FUNCT__ "SNESNGSSetTolerances"
5 /*@
6    SNESNGSSetTolerances - Sets various parameters used in convergence tests.
7 
8    Logically Collective on SNES
9 
10    Input Parameters:
11 +  snes - the SNES context
12 .  abstol - absolute convergence tolerance
13 .  rtol - relative convergence tolerance
14 .  stol -  convergence tolerance in terms of the norm of the change in the solution between steps,  || delta x || < stol*|| x ||
15 -  maxit - maximum number of iterations
16 
17 
18    Options Database Keys:
19 +    -snes_ngs_atol <abstol> - Sets abstol
20 .    -snes_ngs_rtol <rtol> - Sets rtol
21 .    -snes_ngs_stol <stol> - Sets stol
22 -    -snes_max_it <maxit> - Sets maxit
23 
24    Level: intermediate
25 
26 .keywords: SNES, nonlinear, gauss-seidel, set, convergence, tolerances
27 
28 .seealso: SNESSetTrustRegionTolerance()
29 @*/
30 PetscErrorCode  SNESNGSSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit)
31 {
32   SNES_NGS *gs = (SNES_NGS*)snes->data;
33 
34   PetscFunctionBegin;
35   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
36 
37   if (abstol != PETSC_DEFAULT) {
38     if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
39     gs->abstol = abstol;
40   }
41   if (rtol != PETSC_DEFAULT) {
42     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);
43     gs->rtol = rtol;
44   }
45   if (stol != PETSC_DEFAULT) {
46     if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol);
47     gs->stol = stol;
48   }
49   if (maxit != PETSC_DEFAULT) {
50     if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
51     gs->max_its = maxit;
52   }
53   PetscFunctionReturn(0);
54 }
55 
56 #undef __FUNCT__
57 #define __FUNCT__ "SNESNGSGetTolerances"
58 /*@
59    SNESNGSGetTolerances - Gets various parameters used in convergence tests.
60 
61    Not Collective
62 
63    Input Parameters:
64 +  snes - the SNES context
65 .  atol - absolute convergence tolerance
66 .  rtol - relative convergence tolerance
67 .  stol -  convergence tolerance in terms of the norm
68            of the change in the solution between steps
69 -  maxit - maximum number of iterations
70 
71    Notes:
72    The user can specify NULL for any parameter that is not needed.
73 
74    Level: intermediate
75 
76 .keywords: SNES, nonlinear, get, convergence, tolerances
77 
78 .seealso: SNESSetTolerances()
79 @*/
80 PetscErrorCode  SNESNGSGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit)
81 {
82   SNES_NGS *gs = (SNES_NGS*)snes->data;
83 
84   PetscFunctionBegin;
85   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
86   if (atol) *atol = gs->abstol;
87   if (rtol) *rtol = gs->rtol;
88   if (stol) *stol = gs->stol;
89   if (maxit) *maxit = gs->max_its;
90   PetscFunctionReturn(0);
91 }
92 
93 #undef __FUNCT__
94 #define __FUNCT__ "SNESNGSSetSweeps"
95 /*@
96    SNESNGSSetSweeps - Sets the number of sweeps of GS to use.
97 
98    Input Parameters:
99 +  snes   - the SNES context
100 -  sweeps  - the number of sweeps of GS to perform.
101 
102    Level: intermediate
103 
104 .keywords: SNES, nonlinear, set, Gauss-Siedel
105 
106 .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetPC(), SNESNGSGetSweeps()
107 @*/
108 
109 PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps)
110 {
111   SNES_NGS *gs = (SNES_NGS*)snes->data;
112 
113   PetscFunctionBegin;
114   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
115   gs->sweeps = sweeps;
116   PetscFunctionReturn(0);
117 }
118 
119 #undef __FUNCT__
120 #define __FUNCT__ "SNESNGSGetSweeps"
121 /*@
122    SNESNGSGetSweeps - Gets the number of sweeps GS will use.
123 
124    Input Parameters:
125 .  snes   - the SNES context
126 
127    Output Parameters:
128 .  sweeps  - the number of sweeps of GS to perform.
129 
130    Level: intermediate
131 
132 .keywords: SNES, nonlinear, set, Gauss-Siedel
133 
134 .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetPC(), SNESNGSSetSweeps()
135 @*/
136 PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps)
137 {
138   SNES_NGS *gs = (SNES_NGS*)snes->data;
139 
140   PetscFunctionBegin;
141   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
142   *sweeps = gs->sweeps;
143   PetscFunctionReturn(0);
144 }
145 
146 
147 #undef __FUNCT__
148 #define __FUNCT__ "SNESDefaultApplyNGS"
149 PetscErrorCode SNESDefaultApplyNGS(SNES snes,Vec X,Vec F,void *ctx)
150 {
151   PetscFunctionBegin;
152   /* see if there's a coloring on the DM */
153   PetscFunctionReturn(0);
154 }
155 
156 #undef __FUNCT__
157 #define __FUNCT__ "SNESReset_NGS"
158 PetscErrorCode SNESReset_NGS(SNES snes)
159 {
160   PetscFunctionBegin;
161   PetscFunctionReturn(0);
162 }
163 
164 #undef __FUNCT__
165 #define __FUNCT__ "SNESDestroy_NGS"
166 PetscErrorCode SNESDestroy_NGS(SNES snes)
167 {
168   PetscErrorCode ierr;
169 
170   PetscFunctionBegin;
171   ierr = SNESReset_NGS(snes);CHKERRQ(ierr);
172   ierr = PetscFree(snes->data);CHKERRQ(ierr);
173   PetscFunctionReturn(0);
174 }
175 
176 #undef __FUNCT__
177 #define __FUNCT__ "SNESSetUp_NGS"
178 PetscErrorCode SNESSetUp_NGS(SNES snes)
179 {
180   PetscErrorCode ierr;
181   PetscErrorCode (*f)(SNES,Vec,Vec,void*);
182 
183   PetscFunctionBegin;
184   ierr = SNESGetNGS(snes,&f,NULL);CHKERRQ(ierr);
185   if (!f) {
186     ierr = SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);CHKERRQ(ierr);
187   }
188   PetscFunctionReturn(0);
189 }
190 
191 #undef __FUNCT__
192 #define __FUNCT__ "SNESSetFromOptions_NGS"
193 PetscErrorCode SNESSetFromOptions_NGS(PetscOptionItems *PetscOptionsObject,SNES snes)
194 {
195   SNES_NGS       *gs = (SNES_NGS*)snes->data;
196   PetscErrorCode ierr;
197   PetscInt       sweeps,max_its=PETSC_DEFAULT;
198   PetscReal      rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT;
199   PetscBool      flg,flg1,flg2,flg3;
200 
201   PetscFunctionBegin;
202   ierr = PetscOptionsHead(PetscOptionsObject,"SNES GS options");CHKERRQ(ierr);
203   /* GS Options */
204   ierr = PetscOptionsInt("-snes_ngs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",gs->sweeps,&sweeps,&flg);CHKERRQ(ierr);
205   if (flg) {
206     ierr = SNESNGSSetSweeps(snes,sweeps);CHKERRQ(ierr);
207   }
208   ierr = PetscOptionsReal("-snes_ngs_atol","Absolute residual tolerance for GS iteration","SNESComputeGS",gs->abstol,&atol,&flg);CHKERRQ(ierr);
209   ierr = PetscOptionsReal("-snes_ngs_rtol","Relative residual tolerance for GS iteration","SNESComputeGS",gs->rtol,&rtol,&flg1);CHKERRQ(ierr);
210   ierr = PetscOptionsReal("-snes_ngs_stol","Absolute update tolerance for GS iteration","SNESComputeGS",gs->stol,&stol,&flg2);CHKERRQ(ierr);
211   ierr = PetscOptionsInt("-snes_ngs_max_it","Maximum number of sweeps of GS to apply","SNESComputeGS",gs->max_its,&max_its,&flg3);CHKERRQ(ierr);
212   if (flg || flg1 || flg2 || flg3) {
213     ierr = SNESNGSSetTolerances(snes,atol,rtol,stol,max_its);CHKERRQ(ierr);
214   }
215   flg  = PETSC_FALSE;
216   ierr = PetscOptionsBool("-snes_ngs_secant","Use pointwise secant local Jacobian approximation","",flg,&flg,NULL);CHKERRQ(ierr);
217   if (flg) {
218     ierr = SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);CHKERRQ(ierr);
219     ierr = PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");CHKERRQ(ierr);
220   }
221   ierr = PetscOptionsReal("-snes_ngs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);CHKERRQ(ierr);
222   ierr = PetscOptionsBool("-snes_ngs_secant_mat_coloring","Use the Jacobian coloring for the secant GS","",gs->secant_mat,&gs->secant_mat,&flg);CHKERRQ(ierr);
223 
224   ierr = PetscOptionsTail();CHKERRQ(ierr);
225   PetscFunctionReturn(0);
226 }
227 
228 #undef __FUNCT__
229 #define __FUNCT__ "SNESView_NGS"
230 PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer)
231 {
232   PetscFunctionBegin;
233   PetscFunctionReturn(0);
234 }
235 
236 #undef __FUNCT__
237 #define __FUNCT__ "SNESSolve_NGS"
238 PetscErrorCode SNESSolve_NGS(SNES snes)
239 {
240   Vec              F;
241   Vec              X;
242   Vec              B;
243   PetscInt         i;
244   PetscReal        fnorm;
245   PetscErrorCode   ierr;
246   SNESNormSchedule normschedule;
247 
248   PetscFunctionBegin;
249 
250   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);
251 
252   ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
253   X = snes->vec_sol;
254   F = snes->vec_func;
255   B = snes->vec_rhs;
256 
257   ierr         = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
258   snes->iter   = 0;
259   snes->norm   = 0.;
260   ierr         = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
261   snes->reason = SNES_CONVERGED_ITERATING;
262 
263   ierr = SNESGetNormSchedule(snes, &normschedule);CHKERRQ(ierr);
264   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
265     /* compute the initial function and preconditioned update delX */
266     if (!snes->vec_func_init_set) {
267       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
268     } else snes->vec_func_init_set = PETSC_FALSE;
269 
270     ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
271     SNESCheckFunctionNorm(snes,fnorm);
272     ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
273     snes->iter = 0;
274     snes->norm = fnorm;
275     ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
276     ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
277     ierr       = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr);
278 
279     /* test convergence */
280     ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
281     if (snes->reason) PetscFunctionReturn(0);
282   } else {
283     ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
284     ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
285   }
286 
287   /* Call general purpose update function */
288   if (snes->ops->update) {
289     ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
290   }
291 
292   for (i = 0; i < snes->max_its; i++) {
293     ierr = SNESComputeNGS(snes, B, X);CHKERRQ(ierr);
294     /* only compute norms if requested or about to exit due to maximum iterations */
295     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
296       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
297       ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
298       SNESCheckFunctionNorm(snes,fnorm);
299       /* Monitor convergence */
300       ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
301       snes->iter = i+1;
302       snes->norm = fnorm;
303       ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
304       ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
305       ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
306     }
307     /* Test for convergence */
308     if (normschedule == SNES_NORM_ALWAYS) ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
309     if (snes->reason) PetscFunctionReturn(0);
310     /* Call general purpose update function */
311     if (snes->ops->update) {
312       ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
313     }
314   }
315   if (normschedule == SNES_NORM_ALWAYS) {
316     if (i == snes->max_its) {
317       ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);CHKERRQ(ierr);
318       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
319     }
320   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
321   PetscFunctionReturn(0);
322 }
323 
324 /*MC
325   SNESNGS - Just calls the user-provided solution routine provided with SNESSetNGS()
326 
327    Level: advanced
328 
329   Notes:
330   the Gauss-Seidel smoother is inherited through composition.  If a solver has been created with SNESGetPC(), it will have
331   its parent's Gauss-Seidel routine associated with it.
332 
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
336 
337 .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types)
338 M*/
339 
340 #undef __FUNCT__
341 #define __FUNCT__ "SNESCreate_NGS"
342 PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes)
343 {
344   SNES_NGS        *gs;
345   PetscErrorCode ierr;
346 
347   PetscFunctionBegin;
348   snes->ops->destroy        = SNESDestroy_NGS;
349   snes->ops->setup          = SNESSetUp_NGS;
350   snes->ops->setfromoptions = SNESSetFromOptions_NGS;
351   snes->ops->view           = SNESView_NGS;
352   snes->ops->solve          = SNESSolve_NGS;
353   snes->ops->reset          = SNESReset_NGS;
354 
355   snes->usesksp = PETSC_FALSE;
356   snes->usespc  = PETSC_FALSE;
357 
358   if (!snes->tolerancesset) {
359     snes->max_its   = 10000;
360     snes->max_funcs = 10000;
361   }
362 
363   ierr = PetscNewLog(snes,&gs);CHKERRQ(ierr);
364 
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;
371 
372   snes->data = (void*) gs;
373   PetscFunctionReturn(0);
374 }
375