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