xref: /petsc/src/snes/impls/ms/ms.c (revision e07f7f948f87ca2cb8721c32a91a5a8f40991471)
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         = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr);
303   snes->iter   = 0;
304   snes->norm   = 0.;
305   ierr         = PetscObjectAMSGrantAccess((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     if (!snes->norm_init_set) {
320       ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
321       if (PetscIsInfOrNanReal(fnorm)) {
322         snes->reason = SNES_DIVERGED_FNORM_NAN;
323         PetscFunctionReturn(0);
324       }
325     } else {
326       fnorm               = snes->norm_init;
327       snes->norm_init_set = PETSC_FALSE;
328     }
329     /* Monitor convergence */
330     ierr       = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr);
331     snes->iter = 0;
332     snes->norm = fnorm;
333     ierr       = PetscObjectAMSGrantAccess((PetscObject)snes);CHKERRQ(ierr);
334     ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
335     ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
336 
337     /* Test for convergence */
338     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
339     if (snes->reason) PetscFunctionReturn(0);
340   }
341 
342   /* Call general purpose update function */
343   if (snes->ops->update) {
344     ierr = (*snes->ops->update)(snes,snes->iter);CHKERRQ(ierr);
345   }
346   for (i = 0; i < snes->max_its; i++) {
347     ierr = SNESMSStep_3Sstar(snes,X,F);CHKERRQ(ierr);
348 
349     if (i+1 < snes->max_its || ms->norms) {
350       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
351       if (snes->domainerror) {
352         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
353         PetscFunctionReturn(0);
354       }
355     }
356 
357     if (ms->norms) {
358       ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
359       if (PetscIsInfOrNanReal(fnorm)) {
360         snes->reason = SNES_DIVERGED_FNORM_NAN;
361         PetscFunctionReturn(0);
362       }
363 
364       /* Monitor convergence */
365       ierr       = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr);
366       snes->iter = i+1;
367       snes->norm = fnorm;
368       ierr       = PetscObjectAMSGrantAccess((PetscObject)snes);CHKERRQ(ierr);
369       ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
370       ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
371 
372       /* Test for convergence */
373       ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
374       if (snes->reason) PetscFunctionReturn(0);
375     }
376 
377     /* Call general purpose update function */
378     if (snes->ops->update) {
379       ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
380     }
381   }
382   if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
383   PetscFunctionReturn(0);
384 }
385 
386 #undef __FUNCT__
387 #define __FUNCT__ "SNESSetUp_MS"
388 static PetscErrorCode SNESSetUp_MS(SNES snes)
389 {
390   SNES_MS        *ms = (SNES_MS*)snes->data;
391   PetscErrorCode ierr;
392 
393   PetscFunctionBegin;
394   if (!ms->tableau) {ierr = SNESMSSetType(snes,SNESMSDefault);CHKERRQ(ierr);}
395   ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr);
396   ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
397   PetscFunctionReturn(0);
398 }
399 
400 #undef __FUNCT__
401 #define __FUNCT__ "SNESReset_MS"
402 static PetscErrorCode SNESReset_MS(SNES snes)
403 {
404 
405   PetscFunctionBegin;
406   PetscFunctionReturn(0);
407 }
408 
409 #undef __FUNCT__
410 #define __FUNCT__ "SNESDestroy_MS"
411 static PetscErrorCode SNESDestroy_MS(SNES snes)
412 {
413   PetscErrorCode ierr;
414 
415   PetscFunctionBegin;
416   ierr = PetscFree(snes->data);CHKERRQ(ierr);
417   ierr = PetscObjectComposeFunction((PetscObject)snes,"",NULL);CHKERRQ(ierr);
418   PetscFunctionReturn(0);
419 }
420 
421 #undef __FUNCT__
422 #define __FUNCT__ "SNESView_MS"
423 static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
424 {
425   PetscBool      iascii;
426   PetscErrorCode ierr;
427   SNES_MS        *ms = (SNES_MS*)snes->data;
428 
429   PetscFunctionBegin;
430   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
431   if (iascii) {
432     SNESMSTableau tab = ms->tableau;
433     ierr = PetscViewerASCIIPrintf(viewer,"  multi-stage method type: %s\n",tab ? tab->name : "not yet set");CHKERRQ(ierr);
434   }
435   PetscFunctionReturn(0);
436 }
437 
438 #undef __FUNCT__
439 #define __FUNCT__ "SNESSetFromOptions_MS"
440 static PetscErrorCode SNESSetFromOptions_MS(SNES snes)
441 {
442   SNES_MS        *ms = (SNES_MS*)snes->data;
443   PetscErrorCode ierr;
444 
445   PetscFunctionBegin;
446   ierr = PetscOptionsHead("SNES MS options");CHKERRQ(ierr);
447   {
448     SNESMSTableauLink link;
449     PetscInt          count,choice;
450     PetscBool         flg;
451     const char        **namelist;
452     char              mstype[256];
453 
454     ierr = PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));CHKERRQ(ierr);
455     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
456     ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr);
457     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
458     ierr = PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);CHKERRQ(ierr);
459     ierr = SNESMSSetType(snes,flg ? namelist[choice] : mstype);CHKERRQ(ierr);
460     ierr = PetscFree(namelist);CHKERRQ(ierr);
461     ierr = PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,NULL);CHKERRQ(ierr);
462     ierr = PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);CHKERRQ(ierr);
463   }
464   ierr = PetscOptionsTail();CHKERRQ(ierr);
465   PetscFunctionReturn(0);
466 }
467 
468 #undef __FUNCT__
469 #define __FUNCT__ "SNESMSSetType_MS"
470 PetscErrorCode  SNESMSSetType_MS(SNES snes,SNESMSType mstype)
471 {
472   PetscErrorCode    ierr;
473   SNES_MS           *ms = (SNES_MS*)snes->data;
474   SNESMSTableauLink link;
475   PetscBool         match;
476 
477   PetscFunctionBegin;
478   if (ms->tableau) {
479     ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
480     if (match) PetscFunctionReturn(0);
481   }
482   for (link = SNESMSTableauList; link; link=link->next) {
483     ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
484     if (match) {
485       ierr        = SNESReset_MS(snes);CHKERRQ(ierr);
486       ms->tableau = &link->tab;
487       PetscFunctionReturn(0);
488     }
489   }
490   SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
491   PetscFunctionReturn(0);
492 }
493 
494 #undef __FUNCT__
495 #define __FUNCT__ "SNESMSSetType"
496 /*@C
497   SNESMSSetType - Set the type of multistage smoother
498 
499   Logically collective
500 
501   Input Parameter:
502 +  snes - nonlinear solver context
503 -  mstype - type of multistage method
504 
505   Level: beginner
506 
507 .seealso: SNESMSGetType(), SNESMS
508 @*/
509 PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
510 {
511   PetscErrorCode ierr;
512 
513   PetscFunctionBegin;
514   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
515   ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));CHKERRQ(ierr);
516   PetscFunctionReturn(0);
517 }
518 
519 /* -------------------------------------------------------------------------- */
520 /*MC
521       SNESMS - multi-stage smoothers
522 
523       Options Database:
524 
525 +     -snes_ms_type - type of multi-stage smoother
526 -     -snes_ms_damping - damping for multi-stage method
527 
528       Notes:
529       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
530       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
531 
532       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
533 
534       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
535 
536       References:
537 
538       Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
539 
540       Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
541 
542       Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
543 
544       Level: beginner
545 
546 .seealso:  SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
547 
548 M*/
549 #undef __FUNCT__
550 #define __FUNCT__ "SNESCreate_MS"
551 PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
552 {
553   PetscErrorCode ierr;
554   SNES_MS        *ms;
555 
556   PetscFunctionBegin;
557 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
558   ierr = SNESMSInitializePackage();CHKERRQ(ierr);
559 #endif
560 
561   snes->ops->setup          = SNESSetUp_MS;
562   snes->ops->solve          = SNESSolve_MS;
563   snes->ops->destroy        = SNESDestroy_MS;
564   snes->ops->setfromoptions = SNESSetFromOptions_MS;
565   snes->ops->view           = SNESView_MS;
566   snes->ops->reset          = SNESReset_MS;
567 
568   snes->usespc  = PETSC_FALSE;
569   snes->usesksp = PETSC_TRUE;
570 
571   ierr        = PetscNewLog(snes,SNES_MS,&ms);CHKERRQ(ierr);
572   snes->data  = (void*)ms;
573   ms->damping = 0.9;
574   ms->norms   = PETSC_FALSE;
575 
576   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);CHKERRQ(ierr);
577   PetscFunctionReturn(0);
578 }
579 
580