xref: /petsc/src/snes/impls/ms/ms.c (revision 009bbdc485cd9ad46be9940d3549e2dde9cdc322)
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     ierr = PetscFree3(t->gamma,t->delta,t->betasub);CHKERRQ(ierr);
130     ierr = PetscFree(t->name);CHKERRQ(ierr);
131     ierr = PetscFree(link);CHKERRQ(ierr);
132   }
133   SNESMSRegisterAllCalled = PETSC_FALSE;
134   PetscFunctionReturn(0);
135 }
136 
137 #undef __FUNCT__
138 #define __FUNCT__ "SNESMSInitializePackage"
139 /*@C
140   SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
141   from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to SNESCreate_MS()
142   when using static libraries.
143 
144   Input Parameter:
145   path - The dynamic library path, or PETSC_NULL
146 
147   Level: developer
148 
149 .keywords: SNES, SNESMS, initialize, package
150 .seealso: PetscInitialize()
151 @*/
152 PetscErrorCode SNESMSInitializePackage(const char path[])
153 {
154   PetscErrorCode ierr;
155 
156   PetscFunctionBegin;
157   if (SNESMSPackageInitialized) PetscFunctionReturn(0);
158   SNESMSPackageInitialized = PETSC_TRUE;
159   ierr = SNESMSRegisterAll();CHKERRQ(ierr);
160   ierr = PetscRegisterFinalize(SNESMSFinalizePackage);CHKERRQ(ierr);
161   PetscFunctionReturn(0);
162 }
163 
164 #undef __FUNCT__
165 #define __FUNCT__ "SNESMSFinalizePackage"
166 /*@C
167   SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
168   called from PetscFinalize().
169 
170   Level: developer
171 
172 .keywords: Petsc, destroy, package
173 .seealso: PetscFinalize()
174 @*/
175 PetscErrorCode SNESMSFinalizePackage(void)
176 {
177   PetscErrorCode ierr;
178 
179   PetscFunctionBegin;
180   SNESMSPackageInitialized = PETSC_FALSE;
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   ierr = PetscMalloc3(nstages*nregisters,PetscReal,&t->gamma,nstages,PetscReal,&t->delta,nstages,PetscReal,&t->betasub);CHKERRQ(ierr);
231   ierr = PetscMemcpy(t->gamma,gamma,nstages*nregisters*sizeof(PetscReal));CHKERRQ(ierr);
232   ierr = PetscMemcpy(t->delta,delta,nstages*sizeof(PetscReal));CHKERRQ(ierr);
233   ierr = PetscMemcpy(t->betasub,betasub,nstages*sizeof(PetscReal));CHKERRQ(ierr);
234 
235   link->next = SNESMSTableauList;
236   SNESMSTableauList = link;
237   PetscFunctionReturn(0);
238 }
239 
240 #undef __FUNCT__
241 #define __FUNCT__ "SNESMSStep_3Sstar"
242 /*
243   X - initial state, updated in-place.
244   F - residual, computed at the initial X on input
245 */
246 static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
247 {
248   PetscErrorCode ierr;
249   SNES_MS        *ms = (SNES_MS*)snes->data;
250   SNESMSTableau  t = ms->tableau;
251   const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
252   Vec            S1,S2,S3,Y;
253   PetscInt       i,nstages = t->nstages;;
254 
255 
256   PetscFunctionBegin;
257   Y = snes->work[0];
258   S1 = X;
259   S2 = snes->work[1];
260   S3 = snes->work[2];
261   ierr = VecZeroEntries(S2);CHKERRQ(ierr);
262   ierr = VecCopy(X,S3);CHKERRQ(ierr);
263   for (i=0; i<nstages; i++) {
264     Vec Ss[4];
265     PetscScalar scoeff[4];
266 
267     Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;
268     scoeff[0] = gamma[0*nstages+i]-(PetscReal)1.0;
269     scoeff[1] = gamma[1*nstages+i];
270     scoeff[2] = gamma[2*nstages+i];
271     scoeff[3] = -betasub[i]*ms->damping;
272     ierr = VecAXPY(S2,delta[i],S1);CHKERRQ(ierr);
273     if (i>0) {
274       ierr = SNESComputeFunction(snes,S1,F);CHKERRQ(ierr);
275       if (snes->domainerror) {
276         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
277         PetscFunctionReturn(0);
278       }
279     }
280     ierr = SNES_KSPSolve(snes,snes->ksp,F,Y);CHKERRQ(ierr);
281     ierr = VecMAXPY(S1,4,scoeff,Ss);CHKERRQ(ierr);
282   }
283   PetscFunctionReturn(0);
284 }
285 
286 #undef __FUNCT__
287 #define __FUNCT__ "SNESSolve_MS"
288 static PetscErrorCode SNESSolve_MS(SNES snes)
289 {
290   SNES_MS        *ms = (SNES_MS*)snes->data;
291   Vec            X = snes->vec_sol,F = snes->vec_func;
292   PetscReal      fnorm;
293   MatStructure   mstruct;
294   PetscInt       i;
295   PetscErrorCode ierr;
296 
297   PetscFunctionBegin;
298   snes->reason = SNES_CONVERGED_ITERATING;
299   ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
300   snes->iter = 0;
301   snes->norm = 0.;
302   ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
303   if (!snes->vec_func_init_set) {
304     ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
305     if (snes->domainerror) {
306       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
307       PetscFunctionReturn(0);
308     }
309   } else {
310     snes->vec_func_init_set = PETSC_FALSE;
311   }
312   if (snes->jacobian) {         /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
313     ierr = SNESComputeJacobian(snes,snes->vec_sol,&snes->jacobian,&snes->jacobian_pre,&mstruct);CHKERRQ(ierr);
314     ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,mstruct);CHKERRQ(ierr);
315   }
316   if (ms->norms) {
317     if (!snes->norm_init_set) {
318       ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
319       if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
320     } else {
321       fnorm = snes->norm_init;
322       snes->norm_init_set = PETSC_FALSE;
323     }
324     /* Monitor convergence */
325     ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
326     snes->iter = 0;
327     snes->norm = fnorm;
328     ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
329     SNESLogConvHistory(snes,snes->norm,0);
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)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
357       /* Monitor convergence */
358       ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
359       snes->iter = i+1;
360       snes->norm = fnorm;
361       ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
362       SNESLogConvHistory(snes,snes->norm,0);
363       ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
364 
365       /* Test for convergence */
366       ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
367       if (snes->reason) PetscFunctionReturn(0);
368     }
369 
370     /* Call general purpose update function */
371     if (snes->ops->update) {
372       ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
373     }
374   }
375   if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
376   PetscFunctionReturn(0);
377 }
378 
379 #undef __FUNCT__
380 #define __FUNCT__ "SNESSetUp_MS"
381 static PetscErrorCode SNESSetUp_MS(SNES snes)
382 {
383   SNES_MS * ms = (SNES_MS *)snes->data;
384   PetscErrorCode ierr;
385 
386   PetscFunctionBegin;
387   if (!ms->tableau) {ierr = SNESMSSetType(snes,SNESMSDefault);CHKERRQ(ierr);}
388   ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr);
389   ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
390   PetscFunctionReturn(0);
391 }
392 
393 #undef __FUNCT__
394 #define __FUNCT__ "SNESReset_MS"
395 static PetscErrorCode SNESReset_MS(SNES snes)
396 {
397 
398   PetscFunctionBegin;
399   PetscFunctionReturn(0);
400 }
401 
402 #undef __FUNCT__
403 #define __FUNCT__ "SNESDestroy_MS"
404 static PetscErrorCode SNESDestroy_MS(SNES snes)
405 {
406   PetscErrorCode ierr;
407 
408   PetscFunctionBegin;
409   ierr = PetscFree(snes->data);CHKERRQ(ierr);
410   ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"","",PETSC_NULL);CHKERRQ(ierr);
411   PetscFunctionReturn(0);
412 }
413 
414 #undef __FUNCT__
415 #define __FUNCT__ "SNESView_MS"
416 static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
417 {
418   PetscBool        iascii;
419   PetscErrorCode   ierr;
420   SNES_MS          *ms = (SNES_MS*)snes->data;
421 
422   PetscFunctionBegin;
423   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
424   if (iascii) {
425     SNESMSTableau tab = ms->tableau;
426     ierr = PetscViewerASCIIPrintf(viewer,"  multi-stage method type: %s\n",tab?tab->name:"not yet set");CHKERRQ(ierr);
427   }
428   PetscFunctionReturn(0);
429 }
430 
431 #undef __FUNCT__
432 #define __FUNCT__ "SNESSetFromOptions_MS"
433 static PetscErrorCode SNESSetFromOptions_MS(SNES snes)
434 {
435   SNES_MS        *ms = (SNES_MS*)snes->data;
436   PetscErrorCode ierr;
437 
438   PetscFunctionBegin;
439   ierr = PetscOptionsHead("SNES MS options");CHKERRQ(ierr);
440   {
441     SNESMSTableauLink link;
442     PetscInt count,choice;
443     PetscBool flg;
444     const char **namelist;
445     char mstype[256];
446 
447     ierr = PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));CHKERRQ(ierr);
448     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
449     ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr);
450     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
451     ierr = PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);CHKERRQ(ierr);
452     ierr = SNESMSSetType(snes,flg ? namelist[choice] : mstype);CHKERRQ(ierr);
453     ierr = PetscFree(namelist);CHKERRQ(ierr);
454     ierr = PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,PETSC_NULL);CHKERRQ(ierr);
455     ierr = PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,PETSC_NULL);CHKERRQ(ierr);
456   }
457   ierr = PetscOptionsTail();CHKERRQ(ierr);
458   PetscFunctionReturn(0);
459 }
460 
461 EXTERN_C_BEGIN
462 #undef __FUNCT__
463 #define __FUNCT__ "SNESMSSetType_MS"
464 PetscErrorCode  SNESMSSetType_MS(SNES snes,SNESMSType mstype)
465 {
466   PetscErrorCode    ierr;
467   SNES_MS           *ms = (SNES_MS*)snes->data;
468   SNESMSTableauLink link;
469   PetscBool         match;
470 
471   PetscFunctionBegin;
472   if (ms->tableau) {
473     ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
474     if (match) PetscFunctionReturn(0);
475   }
476   for (link = SNESMSTableauList; link; link=link->next) {
477     ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
478     if (match) {
479       ierr = SNESReset_MS(snes);CHKERRQ(ierr);
480       ms->tableau = &link->tab;
481       PetscFunctionReturn(0);
482     }
483   }
484   SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
485   PetscFunctionReturn(0);
486 }
487 EXTERN_C_END
488 
489 
490 #undef __FUNCT__
491 #define __FUNCT__ "SNESMSSetType"
492 /*@C
493   SNESMSSetType - Set the type of multistage smoother
494 
495   Logically collective
496 
497   Input Parameter:
498 +  snes - nonlinear solver context
499 -  mstype - type of multistage method
500 
501   Level: beginner
502 
503 .seealso: SNESMSGetType(), SNESMS
504 @*/
505 PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
506 {
507   PetscErrorCode ierr;
508 
509   PetscFunctionBegin;
510   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
511   ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));CHKERRQ(ierr);
512   PetscFunctionReturn(0);
513 }
514 
515 /* -------------------------------------------------------------------------- */
516 /*MC
517       SNESMS - multi-stage smoothers
518 
519       Options Database:
520 
521 +     -snes_ms_type - type of multi-stage smoother
522 -     -snes_ms_damping - damping for multi-stage method
523 
524       Notes:
525       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
526       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
527 
528       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
529 
530       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
531 
532       References:
533 
534       Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
535 
536       Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
537 
538       Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
539 
540       Level: beginner
541 
542 .seealso:  SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
543 
544 M*/
545 EXTERN_C_BEGIN
546 #undef __FUNCT__
547 #define __FUNCT__ "SNESCreate_MS"
548 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(PETSC_NULL);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 = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESMSSetType_C","SNESMSSetType_MS",SNESMSSetType_MS);CHKERRQ(ierr);
574   PetscFunctionReturn(0);
575 }
576 EXTERN_C_END
577