xref: /petsc/src/snes/impls/ms/ms.c (revision fe998a80077c9ee0917a39496df43fc256e1b478)
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          = PetscNew(&link);CHKERRQ(ierr);
224   t             = &link->tab;
225   ierr          = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
226   t->nstages    = nstages;
227   t->nregisters = nregisters;
228   t->stability  = stability;
229 
230   ierr = PetscMalloc3(nstages*nregisters,&t->gamma,nstages,&t->delta,nstages,&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 
269     scoeff[0] = gamma[0*nstages+i]-(PetscReal)1.0;
270     scoeff[1] = gamma[1*nstages+i];
271     scoeff[2] = gamma[2*nstages+i];
272     scoeff[3] = -betasub[i]*ms->damping;
273 
274     ierr = VecAXPY(S2,delta[i],S1);CHKERRQ(ierr);
275     if (i>0) {
276       ierr = SNESComputeFunction(snes,S1,F);CHKERRQ(ierr);
277       if (snes->domainerror) {
278         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
279         PetscFunctionReturn(0);
280       }
281     }
282     ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr);
283     ierr = VecMAXPY(S1,4,scoeff,Ss);CHKERRQ(ierr);
284   }
285   PetscFunctionReturn(0);
286 }
287 
288 #undef __FUNCT__
289 #define __FUNCT__ "SNESSolve_MS"
290 static PetscErrorCode SNESSolve_MS(SNES snes)
291 {
292   SNES_MS        *ms = (SNES_MS*)snes->data;
293   Vec            X   = snes->vec_sol,F = snes->vec_func;
294   PetscReal      fnorm;
295   PetscInt       i;
296   PetscErrorCode ierr;
297 
298   PetscFunctionBegin;
299   ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
300   snes->reason = SNES_CONVERGED_ITERATING;
301   ierr         = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
302   snes->iter   = 0;
303   snes->norm   = 0.;
304   ierr         = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
305   if (!snes->vec_func_init_set) {
306     ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
307     if (snes->domainerror) {
308       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
309       PetscFunctionReturn(0);
310     }
311   } else snes->vec_func_init_set = PETSC_FALSE;
312 
313   if (snes->jacobian) {         /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
314     ierr = SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
315   }
316   if (ms->norms) {
317     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
318     if (PetscIsInfOrNanReal(fnorm)) {
319       snes->reason = SNES_DIVERGED_FNORM_NAN;
320       PetscFunctionReturn(0);
321     }
322     /* Monitor convergence */
323     ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
324     snes->iter = 0;
325     snes->norm = fnorm;
326     ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
327     ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
328     ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
329 
330     /* Test for convergence */
331     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
332     if (snes->reason) PetscFunctionReturn(0);
333   }
334 
335   /* Call general purpose update function */
336   if (snes->ops->update) {
337     ierr = (*snes->ops->update)(snes,snes->iter);CHKERRQ(ierr);
338   }
339   for (i = 0; i < snes->max_its; i++) {
340     ierr = SNESMSStep_3Sstar(snes,X,F);CHKERRQ(ierr);
341 
342     if (i+1 < snes->max_its || ms->norms) {
343       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
344       if (snes->domainerror) {
345         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
346         PetscFunctionReturn(0);
347       }
348     }
349 
350     if (ms->norms) {
351       ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
352       if (PetscIsInfOrNanReal(fnorm)) {
353         snes->reason = SNES_DIVERGED_FNORM_NAN;
354         PetscFunctionReturn(0);
355       }
356 
357       /* Monitor convergence */
358       ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
359       snes->iter = i+1;
360       snes->norm = fnorm;
361       ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
362       ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
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 = SNESSetWorkVecs(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 = PetscObjectComposeFunction((PetscObject)snes,"",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 = PetscMalloc1(count,&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,NULL);CHKERRQ(ierr);
455     ierr = PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);CHKERRQ(ierr);
456   }
457   ierr = PetscOptionsTail();CHKERRQ(ierr);
458   PetscFunctionReturn(0);
459 }
460 
461 #undef __FUNCT__
462 #define __FUNCT__ "SNESMSSetType_MS"
463 PetscErrorCode  SNESMSSetType_MS(SNES snes,SNESMSType mstype)
464 {
465   PetscErrorCode    ierr;
466   SNES_MS           *ms = (SNES_MS*)snes->data;
467   SNESMSTableauLink link;
468   PetscBool         match;
469 
470   PetscFunctionBegin;
471   if (ms->tableau) {
472     ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
473     if (match) PetscFunctionReturn(0);
474   }
475   for (link = SNESMSTableauList; link; link=link->next) {
476     ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
477     if (match) {
478       ierr        = SNESReset_MS(snes);CHKERRQ(ierr);
479       ms->tableau = &link->tab;
480       PetscFunctionReturn(0);
481     }
482   }
483   SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
484   PetscFunctionReturn(0);
485 }
486 
487 #undef __FUNCT__
488 #define __FUNCT__ "SNESMSSetType"
489 /*@C
490   SNESMSSetType - Set the type of multistage smoother
491 
492   Logically collective
493 
494   Input Parameter:
495 +  snes - nonlinear solver context
496 -  mstype - type of multistage method
497 
498   Level: beginner
499 
500 .seealso: SNESMSGetType(), SNESMS
501 @*/
502 PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
503 {
504   PetscErrorCode ierr;
505 
506   PetscFunctionBegin;
507   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
508   ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));CHKERRQ(ierr);
509   PetscFunctionReturn(0);
510 }
511 
512 /* -------------------------------------------------------------------------- */
513 /*MC
514       SNESMS - multi-stage smoothers
515 
516       Options Database:
517 
518 +     -snes_ms_type - type of multi-stage smoother
519 -     -snes_ms_damping - damping for multi-stage method
520 
521       Notes:
522       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
523       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
524 
525       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
526 
527       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
528 
529       References:
530 
531       Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
532 
533       Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
534 
535       Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
536 
537       Level: beginner
538 
539 .seealso:  SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
540 
541 M*/
542 #undef __FUNCT__
543 #define __FUNCT__ "SNESCreate_MS"
544 PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
545 {
546   PetscErrorCode ierr;
547   SNES_MS        *ms;
548 
549   PetscFunctionBegin;
550   ierr = SNESMSInitializePackage();CHKERRQ(ierr);
551 
552   snes->ops->setup          = SNESSetUp_MS;
553   snes->ops->solve          = SNESSolve_MS;
554   snes->ops->destroy        = SNESDestroy_MS;
555   snes->ops->setfromoptions = SNESSetFromOptions_MS;
556   snes->ops->view           = SNESView_MS;
557   snes->ops->reset          = SNESReset_MS;
558 
559   snes->usespc  = PETSC_FALSE;
560   snes->usesksp = PETSC_TRUE;
561 
562   ierr        = PetscNewLog(snes,&ms);CHKERRQ(ierr);
563   snes->data  = (void*)ms;
564   ms->damping = 0.9;
565   ms->norms   = PETSC_FALSE;
566 
567   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);CHKERRQ(ierr);
568   PetscFunctionReturn(0);
569 }
570 
571