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