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