xref: /petsc/src/snes/impls/ms/ms.c (revision 3847c7255eb26e2c1d88b23c026eee0391de8e56)
1 #include <petsc/private/snesimpl.h>   /*I "petscsnes.h" I*/
2 
3 static SNESMSType SNESMSDefault = SNESMSM62;
4 static PetscBool  SNESMSRegisterAllCalled;
5 static PetscBool  SNESMSPackageInitialized;
6 
7 typedef struct _SNESMSTableau *SNESMSTableau;
8 struct _SNESMSTableau {
9   char      *name;
10   PetscInt  nstages;            /* Number of stages */
11   PetscInt  nregisters;         /* Number of registers */
12   PetscReal stability;          /* Scaled stability region */
13   PetscReal *gamma;             /* Coefficients of 3S* method */
14   PetscReal *delta;             /* Coefficients of 3S* method */
15   PetscReal *betasub;           /* Subdiagonal of beta in Shu-Osher form */
16 };
17 
18 typedef struct _SNESMSTableauLink *SNESMSTableauLink;
19 struct _SNESMSTableauLink {
20   struct _SNESMSTableau tab;
21   SNESMSTableauLink     next;
22 };
23 static SNESMSTableauLink SNESMSTableauList;
24 
25 typedef struct {
26   SNESMSTableau tableau;        /* Tableau in low-storage form */
27   PetscReal     damping;        /* Damping parameter, like length of (pseudo) time step */
28   PetscBool     norms;          /* Compute norms, usually only for monitoring purposes */
29 } SNES_MS;
30 
31 /*@C
32   SNESMSRegisterAll - Registers all of the multi-stage methods in SNESMS
33 
34   Not Collective, but should be called by all processes which will need the schemes to be registered
35 
36   Level: advanced
37 
38 .seealso:  SNESMSRegisterDestroy()
39 @*/
40 PetscErrorCode SNESMSRegisterAll(void)
41 {
42   PetscErrorCode ierr;
43 
44   PetscFunctionBegin;
45   if (SNESMSRegisterAllCalled) PetscFunctionReturn(0);
46   SNESMSRegisterAllCalled = PETSC_TRUE;
47 
48   {
49     PetscReal gamma[3][1] = {{1.0},{0.0},{0.0}};
50     PetscReal delta[1]    = {0.0};
51     PetscReal betasub[1]  = {1.0};
52     ierr = SNESMSRegister(SNESMSEULER,1,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
53   }
54 
55   {
56     PetscReal gamma[3][6] = {
57       {0.0000000000000000E+00,  -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02,  -1.4204296130641869E-01},
58       {1.0000000000000000E+00,  1.1111025767083920E+00,  5.6150921583923230E-01,  7.4151723494934041E-01,  3.1714538168600587E-01,  4.6479276238548706E-01},
59       {0.0000000000000000E+00,  0.0000000000000000E+00,  0.0000000000000000E+00,  6.7968174970583317E-01,  -4.1755042846051737E-03, -1.9115668129923846E-01}
60     };
61     PetscReal delta[6]   = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
62     PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
63     ierr = SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
64   }
65   {
66     PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
67     PetscReal delta[4]    = {0,0,0,0};
68     PetscReal betasub[4]  = {0.25, 0.5, 0.55, 1.0};
69     ierr = SNESMSRegister(SNESMSJAMESON83,4,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
70   }
71 
72   { /* Van Leer, Tai, and Powell (1989) 1 stage, order 1 */
73     PetscReal gamma[3][1] = {{0},{0},{1}};
74     PetscReal delta[1]    = {0};
75     PetscReal betasub[1]  = {1.0};
76     ierr = SNESMSRegister(SNESMSVLTP11,1,3,0.5,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
77   }
78   { /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
79     PetscReal gamma[3][2] = {{0,0},{0,0},{1,1}};
80     PetscReal delta[2]    = {0,0};
81     PetscReal betasub[2]  = {0.3333,1.0};
82     ierr = SNESMSRegister(SNESMSVLTP21,2,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
83   }
84   { /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
85     PetscReal gamma[3][3] = {{0,0,0},{0,0,0},{1,1,1}};
86     PetscReal delta[3]    = {0,0,0};
87     PetscReal betasub[3]  = {0.1481,0.4000,1.0};
88     ierr = SNESMSRegister(SNESMSVLTP31,3,3,1.5,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
89   }
90   { /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
91     PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
92     PetscReal delta[4]    = {0,0,0,0};
93     PetscReal betasub[4]  = {0.0833,0.2069,0.4265,1.0};
94     ierr = SNESMSRegister(SNESMSVLTP41,4,3,2.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
95   }
96   { /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
97     PetscReal gamma[3][5] = {{0,0,0,0,0},{0,0,0,0,0},{1,1,1,1,1}};
98     PetscReal delta[5]    = {0,0,0,0,0};
99     PetscReal betasub[5]  = {0.0533,0.1263,0.2375,0.4414,1.0};
100     ierr = SNESMSRegister(SNESMSVLTP51,5,3,2.5,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
101   }
102   { /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
103     PetscReal gamma[3][6] = {{0,0,0,0,0,0},{0,0,0,0,0,0},{1,1,1,1,1,1}};
104     PetscReal delta[6]    = {0,0,0,0,0,0};
105     PetscReal betasub[6]  = {0.0370,0.0851,0.1521,0.2562,0.4512,1.0};
106     ierr = SNESMSRegister(SNESMSVLTP61,6,3,3.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
107   }
108   PetscFunctionReturn(0);
109 }
110 
111 /*@C
112    SNESMSRegisterDestroy - Frees the list of schemes that were registered by SNESMSRegister().
113 
114    Not Collective
115 
116    Level: advanced
117 
118 .seealso: SNESMSRegister(), SNESMSRegisterAll()
119 @*/
120 PetscErrorCode SNESMSRegisterDestroy(void)
121 {
122   PetscErrorCode    ierr;
123   SNESMSTableauLink link;
124 
125   PetscFunctionBegin;
126   while ((link = SNESMSTableauList)) {
127     SNESMSTableau t = &link->tab;
128     SNESMSTableauList = link->next;
129 
130     ierr = PetscFree3(t->gamma,t->delta,t->betasub);CHKERRQ(ierr);
131     ierr = PetscFree(t->name);CHKERRQ(ierr);
132     ierr = PetscFree(link);CHKERRQ(ierr);
133   }
134   SNESMSRegisterAllCalled = PETSC_FALSE;
135   PetscFunctionReturn(0);
136 }
137 
138 /*@C
139   SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
140   from SNESInitializePackage().
141 
142   Level: developer
143 
144 .seealso: PetscInitialize()
145 @*/
146 PetscErrorCode SNESMSInitializePackage(void)
147 {
148   PetscErrorCode ierr;
149 
150   PetscFunctionBegin;
151   if (SNESMSPackageInitialized) PetscFunctionReturn(0);
152   SNESMSPackageInitialized = PETSC_TRUE;
153 
154   ierr = SNESMSRegisterAll();CHKERRQ(ierr);
155   ierr = PetscRegisterFinalize(SNESMSFinalizePackage);CHKERRQ(ierr);
156   PetscFunctionReturn(0);
157 }
158 
159 /*@C
160   SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
161   called from PetscFinalize().
162 
163   Level: developer
164 
165 .seealso: PetscFinalize()
166 @*/
167 PetscErrorCode SNESMSFinalizePackage(void)
168 {
169   PetscErrorCode ierr;
170 
171   PetscFunctionBegin;
172   SNESMSPackageInitialized = PETSC_FALSE;
173 
174   ierr = SNESMSRegisterDestroy();CHKERRQ(ierr);
175   PetscFunctionReturn(0);
176 }
177 
178 /*@C
179    SNESMSRegister - register a multistage scheme
180 
181    Not Collective, but the same schemes should be registered on all processes on which they will be used
182 
183    Input Parameters:
184 +  name - identifier for method
185 .  nstages - number of stages
186 .  nregisters - number of registers used by low-storage implementation
187 .  gamma - coefficients, see Ketcheson's paper
188 .  delta - coefficients, see Ketcheson's paper
189 -  betasub - subdiagonal of Shu-Osher form
190 
191    Notes:
192    The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
193 
194    Level: advanced
195 
196 .seealso: SNESMS
197 @*/
198 PetscErrorCode SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
199 {
200   PetscErrorCode    ierr;
201   SNESMSTableauLink link;
202   SNESMSTableau     t;
203 
204   PetscFunctionBegin;
205   PetscValidCharPointer(name,1);
206   if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage");
207   if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form");
208   PetscValidRealPointer(gamma,4);
209   PetscValidRealPointer(delta,5);
210   PetscValidRealPointer(betasub,6);
211 
212   ierr          = SNESMSInitializePackage();CHKERRQ(ierr);
213   ierr          = PetscNew(&link);CHKERRQ(ierr);
214   t             = &link->tab;
215   ierr          = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
216   t->nstages    = nstages;
217   t->nregisters = nregisters;
218   t->stability  = stability;
219 
220   ierr = PetscMalloc3(nstages*nregisters,&t->gamma,nstages,&t->delta,nstages,&t->betasub);CHKERRQ(ierr);
221   ierr = PetscArraycpy(t->gamma,gamma,nstages*nregisters);CHKERRQ(ierr);
222   ierr = PetscArraycpy(t->delta,delta,nstages);CHKERRQ(ierr);
223   ierr = PetscArraycpy(t->betasub,betasub,nstages);CHKERRQ(ierr);
224 
225   link->next        = SNESMSTableauList;
226   SNESMSTableauList = link;
227   PetscFunctionReturn(0);
228 }
229 
230 /*
231   X - initial state, updated in-place.
232   F - residual, computed at the initial X on input
233 */
234 static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
235 {
236   PetscErrorCode  ierr;
237   SNES_MS         *ms    = (SNES_MS*)snes->data;
238   SNESMSTableau   t      = ms->tableau;
239   const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
240   Vec             S1,S2,S3,Y;
241   PetscInt        i,nstages = t->nstages;
242 
243 
244   PetscFunctionBegin;
245   Y    = snes->work[0];
246   S1   = X;
247   S2   = snes->work[1];
248   S3   = snes->work[2];
249   ierr = VecZeroEntries(S2);CHKERRQ(ierr);
250   ierr = VecCopy(X,S3);CHKERRQ(ierr);
251   for (i = 0; i < nstages; i++) {
252     Vec         Ss[4];
253     PetscScalar scoeff[4];
254 
255     Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;
256 
257     scoeff[0] = gamma[0*nstages+i] - 1;
258     scoeff[1] = gamma[1*nstages+i];
259     scoeff[2] = gamma[2*nstages+i];
260     scoeff[3] = -betasub[i]*ms->damping;
261 
262     ierr = VecAXPY(S2,delta[i],S1);CHKERRQ(ierr);
263     if (i > 0) {
264       ierr = SNESComputeFunction(snes,S1,F);CHKERRQ(ierr);
265     }
266     ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr);
267     ierr = VecMAXPY(S1,4,scoeff,Ss);CHKERRQ(ierr);
268   }
269   PetscFunctionReturn(0);
270 }
271 
272 static PetscErrorCode SNESMSStep_Norms(SNES snes,PetscInt iter,Vec F)
273 {
274   SNES_MS        *ms = (SNES_MS*)snes->data;
275   PetscReal      fnorm;
276   PetscErrorCode ierr;
277 
278   PetscFunctionBegin;
279   if (ms->norms) {
280     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
281     SNESCheckFunctionNorm(snes,fnorm);
282     /* Monitor convergence */
283     ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
284     snes->iter = iter;
285     snes->norm = fnorm;
286     ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
287     ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
288     ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
289     /* Test for convergence */
290     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
291   } else if (iter > 0) {
292     ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
293     snes->iter = iter;
294     ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
295   }
296   PetscFunctionReturn(0);
297 }
298 
299 
300 static PetscErrorCode SNESSolve_MS(SNES snes)
301 {
302   SNES_MS        *ms = (SNES_MS*)snes->data;
303   Vec            X   = snes->vec_sol,F = snes->vec_func;
304   PetscInt       i;
305   PetscErrorCode ierr;
306 
307   PetscFunctionBegin;
308   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);
309   ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
310 
311   snes->reason = SNES_CONVERGED_ITERATING;
312   ierr         = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
313   snes->iter   = 0;
314   snes->norm   = 0;
315   ierr         = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
316 
317   if (!snes->vec_func_init_set) {
318     ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
319   } else snes->vec_func_init_set = PETSC_FALSE;
320 
321   ierr = SNESMSStep_Norms(snes,0,F);CHKERRQ(ierr);
322   if (snes->reason) PetscFunctionReturn(0);
323 
324   for (i = 0; i < snes->max_its; i++) {
325 
326     /* Call general purpose update function */
327     if (snes->ops->update) {
328       ierr = (*snes->ops->update)(snes,snes->iter);CHKERRQ(ierr);
329     }
330 
331     if (i == 0 && snes->jacobian) {
332       /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
333       ierr = SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
334       SNESCheckJacobianDomainerror(snes);
335       ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
336     }
337 
338     ierr = SNESMSStep_3Sstar(snes,X,F);CHKERRQ(ierr);
339 
340     if (i < snes->max_its-1 || ms->norms) {
341       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
342     }
343 
344     ierr = SNESMSStep_Norms(snes,i+1,F);CHKERRQ(ierr);
345     if (snes->reason) PetscFunctionReturn(0);
346   }
347   if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
348   PetscFunctionReturn(0);
349 }
350 
351 static PetscErrorCode SNESSetUp_MS(SNES snes)
352 {
353   PetscErrorCode ierr;
354 
355   PetscFunctionBegin;
356   ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr);
357   ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
358   PetscFunctionReturn(0);
359 }
360 
361 static PetscErrorCode SNESReset_MS(SNES snes)
362 {
363   PetscFunctionBegin;
364   PetscFunctionReturn(0);
365 }
366 
367 static PetscErrorCode SNESDestroy_MS(SNES snes)
368 {
369   PetscErrorCode ierr;
370 
371   PetscFunctionBegin;
372   ierr = SNESReset_MS(snes);CHKERRQ(ierr);
373   ierr = PetscFree(snes->data);CHKERRQ(ierr);
374   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetType_C",NULL);CHKERRQ(ierr);
375   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",NULL);CHKERRQ(ierr);
376   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetDamping_C",NULL);CHKERRQ(ierr);
377   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetDamping_C",NULL);CHKERRQ(ierr);
378   PetscFunctionReturn(0);
379 }
380 
381 static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
382 {
383   PetscBool      iascii;
384   PetscErrorCode ierr;
385   SNES_MS        *ms = (SNES_MS*)snes->data;
386   SNESMSTableau  tab = ms->tableau;
387 
388   PetscFunctionBegin;
389   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
390   if (iascii) {
391     ierr = PetscViewerASCIIPrintf(viewer,"  multi-stage method type: %s\n",tab->name);CHKERRQ(ierr);
392   }
393   PetscFunctionReturn(0);
394 }
395 
396 static PetscErrorCode SNESSetFromOptions_MS(PetscOptionItems *PetscOptionsObject,SNES snes)
397 {
398   SNES_MS        *ms = (SNES_MS*)snes->data;
399   PetscErrorCode ierr;
400 
401   PetscFunctionBegin;
402   ierr = PetscOptionsHead(PetscOptionsObject,"SNES MS options");CHKERRQ(ierr);
403   {
404     SNESMSTableauLink link;
405     PetscInt          count,choice;
406     PetscBool         flg;
407     const char        **namelist;
408     SNESMSType        mstype;
409     PetscReal         damping;
410 
411     ierr = SNESMSGetType(snes,&mstype);CHKERRQ(ierr);
412     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
413     ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr);
414     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
415     ierr = PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);CHKERRQ(ierr);
416     if (flg) {ierr = SNESMSSetType(snes,namelist[choice]);CHKERRQ(ierr);}
417     ierr = PetscFree(namelist);CHKERRQ(ierr);
418     ierr = SNESMSGetDamping(snes,&damping);CHKERRQ(ierr);
419     ierr = PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",damping,&damping,&flg);CHKERRQ(ierr);
420     if (flg) {ierr = SNESMSSetDamping(snes,damping);CHKERRQ(ierr);}
421     ierr = PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);CHKERRQ(ierr);
422   }
423   ierr = PetscOptionsTail();CHKERRQ(ierr);
424   PetscFunctionReturn(0);
425 }
426 
427 static PetscErrorCode SNESMSGetType_MS(SNES snes,SNESMSType *mstype)
428 {
429   SNES_MS        *ms = (SNES_MS*)snes->data;
430   SNESMSTableau  tab = ms->tableau;
431 
432   PetscFunctionBegin;
433   *mstype = tab->name;
434   PetscFunctionReturn(0);
435 }
436 
437 static PetscErrorCode SNESMSSetType_MS(SNES snes,SNESMSType mstype)
438 {
439   SNES_MS           *ms = (SNES_MS*)snes->data;
440   SNESMSTableauLink link;
441   PetscBool         match;
442   PetscErrorCode    ierr;
443 
444   PetscFunctionBegin;
445   if (ms->tableau) {
446     ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
447     if (match) PetscFunctionReturn(0);
448   }
449   for (link = SNESMSTableauList; link; link=link->next) {
450     ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
451     if (match) {
452       if (snes->setupcalled)  {ierr = SNESReset_MS(snes);CHKERRQ(ierr);}
453       ms->tableau = &link->tab;
454       if (snes->setupcalled)  {ierr = SNESSetUp_MS(snes);CHKERRQ(ierr);}
455       PetscFunctionReturn(0);
456     }
457   }
458   SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
459   PetscFunctionReturn(0);
460 }
461 
462 /*@C
463   SNESMSGetType - Get the type of multistage smoother
464 
465   Not collective
466 
467   Input Parameter:
468 .  snes - nonlinear solver context
469 
470   Output Parameter:
471 .  mstype - type of multistage method
472 
473   Level: beginner
474 
475 .seealso: SNESMSSetType(), SNESMSType, SNESMS
476 @*/
477 PetscErrorCode SNESMSGetType(SNES snes,SNESMSType *mstype)
478 {
479   PetscErrorCode ierr;
480 
481   PetscFunctionBegin;
482   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
483   PetscValidPointer(mstype,2);
484   ierr = PetscUseMethod(snes,"SNESMSGetType_C",(SNES,SNESMSType*),(snes,mstype));CHKERRQ(ierr);
485   PetscFunctionReturn(0);
486 }
487 
488 /*@C
489   SNESMSSetType - Set the type of multistage smoother
490 
491   Logically collective
492 
493   Input Parameter:
494 +  snes - nonlinear solver context
495 -  mstype - type of multistage method
496 
497   Level: beginner
498 
499 .seealso: SNESMSGetType(), SNESMSType, SNESMS
500 @*/
501 PetscErrorCode SNESMSSetType(SNES snes,SNESMSType mstype)
502 {
503   PetscErrorCode ierr;
504 
505   PetscFunctionBegin;
506   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
507   PetscValidPointer(mstype,2);
508   ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,mstype));CHKERRQ(ierr);
509   PetscFunctionReturn(0);
510 }
511 
512 static PetscErrorCode SNESMSGetDamping_MS(SNES snes,PetscReal *damping)
513 {
514   SNES_MS        *ms = (SNES_MS*)snes->data;
515 
516   PetscFunctionBegin;
517   *damping = ms->damping;
518   PetscFunctionReturn(0);
519 }
520 
521 static PetscErrorCode SNESMSSetDamping_MS(SNES snes,PetscReal damping)
522 {
523   SNES_MS           *ms = (SNES_MS*)snes->data;
524 
525   PetscFunctionBegin;
526   ms->damping = damping;
527   PetscFunctionReturn(0);
528 }
529 
530 /*@C
531   SNESMSGetDamping - Get the damping parameter
532 
533   Not collective
534 
535   Input Parameter:
536 .  snes - nonlinear solver context
537 
538   Output Parameter:
539 .  damping - damping parameter
540 
541   Level: advanced
542 
543 .seealso: SNESMSSetDamping(), SNESMS
544 @*/
545 PetscErrorCode SNESMSGetDamping(SNES snes,PetscReal *damping)
546 {
547   PetscErrorCode ierr;
548 
549   PetscFunctionBegin;
550   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
551   PetscValidPointer(damping,2);
552   ierr = PetscUseMethod(snes,"SNESMSGetDamping_C",(SNES,PetscReal*),(snes,damping));CHKERRQ(ierr);
553   PetscFunctionReturn(0);
554 }
555 
556 /*@C
557   SNESMSSetDamping - Set the damping parameter
558 
559   Logically collective
560 
561   Input Parameter:
562 +  snes - nonlinear solver context
563 -  damping - damping parameter
564 
565   Level: advanced
566 
567 .seealso: SNESMSGetDamping(), SNESMS
568 @*/
569 PetscErrorCode SNESMSSetDamping(SNES snes,PetscReal damping)
570 {
571   PetscErrorCode ierr;
572 
573   PetscFunctionBegin;
574   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
575   PetscValidLogicalCollectiveReal(snes,damping,2);
576   ierr = PetscTryMethod(snes,"SNESMSSetDamping_C",(SNES,PetscReal),(snes,damping));CHKERRQ(ierr);
577   PetscFunctionReturn(0);
578 }
579 
580 /* -------------------------------------------------------------------------- */
581 /*MC
582       SNESMS - multi-stage smoothers
583 
584       Options Database:
585 
586 +     -snes_ms_type - type of multi-stage smoother
587 -     -snes_ms_damping - damping for multi-stage method
588 
589       Notes:
590       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
591       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
592 
593       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
594 
595       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
596 
597       References:
598 +   1. -   Ketcheson (2010) Runge Kutta methods with minimum storage implementations (https://doi.org/10.1016/j.jcp.2009.11.006).
599 .   2. -   Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method (https://doi.org/10.1016/0096-3003(83)90019-X).
600 .   3. -   Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes (https://doi.org/10.1006/jcph.1997.5772).
601 -   4. -   Van Leer, Tai, and Powell (1989) Design of optimally smoothing multi-stage schemes for the Euler equations (https://doi.org/10.2514/6.1989-1933).
602 
603       Level: beginner
604 
605 .seealso:  SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
606 
607 M*/
608 PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
609 {
610   PetscErrorCode ierr;
611   SNES_MS        *ms;
612 
613   PetscFunctionBegin;
614   ierr = SNESMSInitializePackage();CHKERRQ(ierr);
615 
616   snes->ops->setup          = SNESSetUp_MS;
617   snes->ops->solve          = SNESSolve_MS;
618   snes->ops->destroy        = SNESDestroy_MS;
619   snes->ops->setfromoptions = SNESSetFromOptions_MS;
620   snes->ops->view           = SNESView_MS;
621   snes->ops->reset          = SNESReset_MS;
622 
623   snes->usesnpc = PETSC_FALSE;
624   snes->usesksp = PETSC_TRUE;
625 
626   snes->alwayscomputesfinalresidual = PETSC_FALSE;
627 
628   ierr        = PetscNewLog(snes,&ms);CHKERRQ(ierr);
629   snes->data  = (void*)ms;
630   ms->damping = 0.9;
631   ms->norms   = PETSC_FALSE;
632 
633   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetType_C",SNESMSGetType_MS);CHKERRQ(ierr);
634   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);CHKERRQ(ierr);
635   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetDamping_C",SNESMSGetDamping_MS);CHKERRQ(ierr);
636   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetDamping_C",SNESMSSetDamping_MS);CHKERRQ(ierr);
637 
638   ierr = SNESMSSetType(snes,SNESMSDefault);CHKERRQ(ierr);
639   PetscFunctionReturn(0);
640 }
641 
642