xref: /petsc/src/snes/impls/ms/ms.c (revision 5a856986583887c326abe5dfd149e8184a29cd80)
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 /*@C
32   SNESMSRegisterAll - Registers all of the multi-stage methods in SNESMS
33 
34   Not Collective, but should be called by all processes which will need the schemes to be registered
35 
36   Level: advanced
37 
38 .seealso:  SNESMSRegisterDestroy()
39 @*/
40 PetscErrorCode SNESMSRegisterAll(void)
41 {
42   PetscErrorCode ierr;
43 
44   PetscFunctionBegin;
45   if (SNESMSRegisterAllCalled) PetscFunctionReturn(0);
46   SNESMSRegisterAllCalled = PETSC_TRUE;
47 
48   {
49     PetscReal gamma[3][1] = {{1.0},{0.0},{0.0}};
50     PetscReal delta[1]    = {0.0};
51     PetscReal betasub[1]  = {1.0};
52     ierr = SNESMSRegister(SNESMSEULER,1,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
53   }
54 
55   {
56     PetscReal gamma[3][6] = {
57       {0.0000000000000000E+00,  -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02,  -1.4204296130641869E-01},
58       {1.0000000000000000E+00,  1.1111025767083920E+00,  5.6150921583923230E-01,  7.4151723494934041E-01,  3.1714538168600587E-01,  4.6479276238548706E-01},
59       {0.0000000000000000E+00,  0.0000000000000000E+00,  0.0000000000000000E+00,  6.7968174970583317E-01,  -4.1755042846051737E-03, -1.9115668129923846E-01}
60     };
61     PetscReal delta[6]   = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
62     PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
63     ierr = SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
64   }
65   {
66     PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
67     PetscReal delta[4]    = {0,0,0,0};
68     PetscReal betasub[4]  = {0.25, 0.5, 0.55, 1.0};
69     ierr = SNESMSRegister(SNESMSJAMESON83,4,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
70   }
71   {                             /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
72     PetscReal gamma[3][2] = {{0,0},{0,0},{1,1}};
73     PetscReal delta[2]    = {0,0};
74     PetscReal betasub[2]  = {0.3333,1.0};
75     ierr = SNESMSRegister(SNESMSVLTP21,2,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
76   }
77   {                             /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
78     PetscReal gamma[3][3] = {{0,0,0},{0,0,0},{1,1,1}};
79     PetscReal delta[3]    = {0,0,0};
80     PetscReal betasub[3]  = {0.1481,0.4000,1.0};
81     ierr = SNESMSRegister(SNESMSVLTP31,3,3,1.5,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
82   }
83   {                             /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
84     PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
85     PetscReal delta[4]    = {0,0,0,0};
86     PetscReal betasub[4]  = {0.0833,0.2069,0.4265,1.0};
87     ierr = SNESMSRegister(SNESMSVLTP41,4,3,2.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
88   }
89   {                             /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
90     PetscReal gamma[3][5] = {{0,0,0,0,0},{0,0,0,0,0},{1,1,1,1,1}};
91     PetscReal delta[5]    = {0,0,0,0,0};
92     PetscReal betasub[5]  = {0.0533,0.1263,0.2375,0.4414,1.0};
93     ierr = SNESMSRegister(SNESMSVLTP51,5,3,2.5,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
94   }
95   {                             /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
96     PetscReal gamma[3][6] = {{0,0,0,0,0,0},{0,0,0,0,0,0},{1,1,1,1,1,1}};
97     PetscReal delta[6]    = {0,0,0,0,0,0};
98     PetscReal betasub[6]  = {0.0370,0.0851,0.1521,0.2562,0.4512,1.0};
99     ierr = SNESMSRegister(SNESMSVLTP61,6,3,3.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
100   }
101   PetscFunctionReturn(0);
102 }
103 
104 /*@C
105    SNESMSRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister().
106 
107    Not Collective
108 
109    Level: advanced
110 
111 .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegister()
112 @*/
113 PetscErrorCode SNESMSRegisterDestroy(void)
114 {
115   PetscErrorCode    ierr;
116   SNESMSTableauLink link;
117 
118   PetscFunctionBegin;
119   while ((link = SNESMSTableauList)) {
120     SNESMSTableau t = &link->tab;
121     SNESMSTableauList = link->next;
122 
123     ierr = PetscFree3(t->gamma,t->delta,t->betasub);CHKERRQ(ierr);
124     ierr = PetscFree(t->name);CHKERRQ(ierr);
125     ierr = PetscFree(link);CHKERRQ(ierr);
126   }
127   SNESMSRegisterAllCalled = PETSC_FALSE;
128   PetscFunctionReturn(0);
129 }
130 
131 /*@C
132   SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
133   from SNESInitializePackage().
134 
135   Level: developer
136 
137 .seealso: PetscInitialize()
138 @*/
139 PetscErrorCode SNESMSInitializePackage(void)
140 {
141   PetscErrorCode ierr;
142 
143   PetscFunctionBegin;
144   if (SNESMSPackageInitialized) PetscFunctionReturn(0);
145   SNESMSPackageInitialized = PETSC_TRUE;
146 
147   ierr = SNESMSRegisterAll();CHKERRQ(ierr);
148   ierr = PetscRegisterFinalize(SNESMSFinalizePackage);CHKERRQ(ierr);
149   PetscFunctionReturn(0);
150 }
151 
152 /*@C
153   SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
154   called from PetscFinalize().
155 
156   Level: developer
157 
158 .seealso: PetscFinalize()
159 @*/
160 PetscErrorCode SNESMSFinalizePackage(void)
161 {
162   PetscErrorCode ierr;
163 
164   PetscFunctionBegin;
165   SNESMSPackageInitialized = PETSC_FALSE;
166 
167   ierr = SNESMSRegisterDestroy();CHKERRQ(ierr);
168   PetscFunctionReturn(0);
169 }
170 
171 /*@C
172    SNESMSRegister - register a multistage scheme
173 
174    Not Collective, but the same schemes should be registered on all processes on which they will be used
175 
176    Input Parameters:
177 +  name - identifier for method
178 .  nstages - number of stages
179 .  nregisters - number of registers used by low-storage implementation
180 .  gamma - coefficients, see Ketcheson's paper
181 .  delta - coefficients, see Ketcheson's paper
182 -  betasub - subdiagonal of Shu-Osher form
183 
184    Notes:
185    The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
186 
187    Level: advanced
188 
189 .seealso: SNESMS
190 @*/
191 PetscErrorCode SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
192 {
193   PetscErrorCode    ierr;
194   SNESMSTableauLink link;
195   SNESMSTableau     t;
196 
197   PetscFunctionBegin;
198   PetscValidCharPointer(name,1);
199   if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage");
200   if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form");
201   PetscValidPointer(gamma,4);
202   PetscValidPointer(delta,5);
203   PetscValidPointer(betasub,6);
204 
205   ierr          = SNESMSInitializePackage();CHKERRQ(ierr);
206   ierr          = PetscNew(&link);CHKERRQ(ierr);
207   t             = &link->tab;
208   ierr          = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
209   t->nstages    = nstages;
210   t->nregisters = nregisters;
211   t->stability  = stability;
212 
213   ierr = PetscMalloc3(nstages*nregisters,&t->gamma,nstages,&t->delta,nstages,&t->betasub);CHKERRQ(ierr);
214   ierr = PetscArraycpy(t->gamma,gamma,nstages*nregisters);CHKERRQ(ierr);
215   ierr = PetscArraycpy(t->delta,delta,nstages);CHKERRQ(ierr);
216   ierr = PetscArraycpy(t->betasub,betasub,nstages);CHKERRQ(ierr);
217 
218   link->next        = SNESMSTableauList;
219   SNESMSTableauList = link;
220   PetscFunctionReturn(0);
221 }
222 
223 /*
224   X - initial state, updated in-place.
225   F - residual, computed at the initial X on input
226 */
227 static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
228 {
229   PetscErrorCode  ierr;
230   SNES_MS         *ms    = (SNES_MS*)snes->data;
231   SNESMSTableau   t      = ms->tableau;
232   const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
233   Vec             S1,S2,S3,Y;
234   PetscInt        i,nstages = t->nstages;;
235 
236 
237   PetscFunctionBegin;
238   Y    = snes->work[0];
239   S1   = X;
240   S2   = snes->work[1];
241   S3   = snes->work[2];
242   ierr = VecZeroEntries(S2);CHKERRQ(ierr);
243   ierr = VecCopy(X,S3);CHKERRQ(ierr);
244   for (i=0; i<nstages; i++) {
245     Vec         Ss[4];
246     PetscScalar scoeff[4];
247 
248     Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;
249 
250     scoeff[0] = gamma[0*nstages+i]-(PetscReal)1.0;
251     scoeff[1] = gamma[1*nstages+i];
252     scoeff[2] = gamma[2*nstages+i];
253     scoeff[3] = -betasub[i]*ms->damping;
254 
255     ierr = VecAXPY(S2,delta[i],S1);CHKERRQ(ierr);
256     if (i>0) {
257       ierr = SNESComputeFunction(snes,S1,F);CHKERRQ(ierr);
258     }
259     ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr);
260     ierr = VecMAXPY(S1,4,scoeff,Ss);CHKERRQ(ierr);
261   }
262   PetscFunctionReturn(0);
263 }
264 
265 static PetscErrorCode SNESSolve_MS(SNES snes)
266 {
267   SNES_MS        *ms = (SNES_MS*)snes->data;
268   Vec            X   = snes->vec_sol,F = snes->vec_func;
269   PetscReal      fnorm;
270   PetscInt       i;
271   PetscErrorCode ierr;
272 
273   PetscFunctionBegin;
274   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
275 
276   ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
277   snes->reason = SNES_CONVERGED_ITERATING;
278   ierr         = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
279   snes->iter   = 0;
280   snes->norm   = 0.;
281   ierr         = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
282   if (!snes->vec_func_init_set) {
283     ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
284   } else snes->vec_func_init_set = PETSC_FALSE;
285 
286   if (snes->jacobian) {         /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
287     ierr = SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
288     SNESCheckJacobianDomainerror(snes);
289   }
290   if (ms->norms) {
291     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
292     SNESCheckFunctionNorm(snes,fnorm);
293     /* Monitor convergence */
294     ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
295     snes->iter = 0;
296     snes->norm = fnorm;
297     ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
298     ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
299     ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
300 
301     /* Test for convergence */
302     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
303     if (snes->reason) PetscFunctionReturn(0);
304   }
305 
306   /* Call general purpose update function */
307   if (snes->ops->update) {
308     ierr = (*snes->ops->update)(snes,snes->iter);CHKERRQ(ierr);
309   }
310   for (i = 0; i < snes->max_its; i++) {
311     ierr = SNESMSStep_3Sstar(snes,X,F);CHKERRQ(ierr);
312 
313     if (i+1 < snes->max_its || ms->norms) {
314       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
315     }
316 
317     if (ms->norms) {
318       ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
319       SNESCheckFunctionNorm(snes,fnorm);
320 
321       /* Monitor convergence */
322       ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
323       snes->iter = i+1;
324       snes->norm = fnorm;
325       ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
326       ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
327       ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
328 
329       /* Test for convergence */
330       ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
331       if (snes->reason) PetscFunctionReturn(0);
332     }
333 
334     /* Call general purpose update function */
335     if (snes->ops->update) {
336       ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
337     }
338   }
339   if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
340   PetscFunctionReturn(0);
341 }
342 
343 static PetscErrorCode SNESSetUp_MS(SNES snes)
344 {
345   SNES_MS        *ms = (SNES_MS*)snes->data;
346   PetscErrorCode ierr;
347 
348   PetscFunctionBegin;
349   if (!ms->tableau) {ierr = SNESMSSetType(snes,SNESMSDefault);CHKERRQ(ierr);}
350   ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr);
351   ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
352   PetscFunctionReturn(0);
353 }
354 
355 static PetscErrorCode SNESReset_MS(SNES snes)
356 {
357 
358   PetscFunctionBegin;
359   PetscFunctionReturn(0);
360 }
361 
362 static PetscErrorCode SNESDestroy_MS(SNES snes)
363 {
364   PetscErrorCode ierr;
365 
366   PetscFunctionBegin;
367   ierr = PetscFree(snes->data);CHKERRQ(ierr);
368   ierr = PetscObjectComposeFunction((PetscObject)snes,"",NULL);CHKERRQ(ierr);
369   PetscFunctionReturn(0);
370 }
371 
372 static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
373 {
374   PetscBool      iascii;
375   PetscErrorCode ierr;
376   SNES_MS        *ms = (SNES_MS*)snes->data;
377 
378   PetscFunctionBegin;
379   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
380   if (iascii) {
381     SNESMSTableau tab = ms->tableau;
382     ierr = PetscViewerASCIIPrintf(viewer,"  multi-stage method type: %s\n",tab ? tab->name : "not yet set");CHKERRQ(ierr);
383   }
384   PetscFunctionReturn(0);
385 }
386 
387 static PetscErrorCode SNESSetFromOptions_MS(PetscOptionItems *PetscOptionsObject,SNES snes)
388 {
389   SNES_MS        *ms = (SNES_MS*)snes->data;
390   PetscErrorCode ierr;
391 
392   PetscFunctionBegin;
393   ierr = PetscOptionsHead(PetscOptionsObject,"SNES MS options");CHKERRQ(ierr);
394   {
395     SNESMSTableauLink link;
396     PetscInt          count,choice;
397     PetscBool         flg;
398     const char        **namelist;
399     char              mstype[256];
400 
401     ierr = PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));CHKERRQ(ierr);
402     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
403     ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr);
404     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
405     ierr = PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);CHKERRQ(ierr);
406     ierr = SNESMSSetType(snes,flg ? namelist[choice] : mstype);CHKERRQ(ierr);
407     ierr = PetscFree(namelist);CHKERRQ(ierr);
408     ierr = PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,NULL);CHKERRQ(ierr);
409     ierr = PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);CHKERRQ(ierr);
410   }
411   ierr = PetscOptionsTail();CHKERRQ(ierr);
412   PetscFunctionReturn(0);
413 }
414 
415 PetscErrorCode  SNESMSSetType_MS(SNES snes,SNESMSType mstype)
416 {
417   PetscErrorCode    ierr;
418   SNES_MS           *ms = (SNES_MS*)snes->data;
419   SNESMSTableauLink link;
420   PetscBool         match;
421 
422   PetscFunctionBegin;
423   if (ms->tableau) {
424     ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
425     if (match) PetscFunctionReturn(0);
426   }
427   for (link = SNESMSTableauList; link; link=link->next) {
428     ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
429     if (match) {
430       ierr        = SNESReset_MS(snes);CHKERRQ(ierr);
431       ms->tableau = &link->tab;
432       PetscFunctionReturn(0);
433     }
434   }
435   SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
436   PetscFunctionReturn(0);
437 }
438 
439 /*@C
440   SNESMSSetType - Set the type of multistage smoother
441 
442   Logically collective
443 
444   Input Parameter:
445 +  snes - nonlinear solver context
446 -  mstype - type of multistage method
447 
448   Level: beginner
449 
450 .seealso: SNESMSGetType(), SNESMS
451 @*/
452 PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
453 {
454   PetscErrorCode ierr;
455 
456   PetscFunctionBegin;
457   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
458   ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));CHKERRQ(ierr);
459   PetscFunctionReturn(0);
460 }
461 
462 /* -------------------------------------------------------------------------- */
463 /*MC
464       SNESMS - multi-stage smoothers
465 
466       Options Database:
467 
468 +     -snes_ms_type - type of multi-stage smoother
469 -     -snes_ms_damping - damping for multi-stage method
470 
471       Notes:
472       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
473       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
474 
475       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
476 
477       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
478 
479       References:
480 +   1. -   Ketcheson (2010) Runge Kutta methods with minimum storage implementations.
481 .   2. -   Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
482 -   3. -   Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
483 
484       Level: beginner
485 
486 .seealso:  SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
487 
488 M*/
489 PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
490 {
491   PetscErrorCode ierr;
492   SNES_MS        *ms;
493 
494   PetscFunctionBegin;
495   ierr = SNESMSInitializePackage();CHKERRQ(ierr);
496 
497   snes->ops->setup          = SNESSetUp_MS;
498   snes->ops->solve          = SNESSolve_MS;
499   snes->ops->destroy        = SNESDestroy_MS;
500   snes->ops->setfromoptions = SNESSetFromOptions_MS;
501   snes->ops->view           = SNESView_MS;
502   snes->ops->reset          = SNESReset_MS;
503 
504   snes->usesnpc = PETSC_FALSE;
505   snes->usesksp = PETSC_TRUE;
506 
507   snes->alwayscomputesfinalresidual = PETSC_FALSE;
508 
509   ierr        = PetscNewLog(snes,&ms);CHKERRQ(ierr);
510   snes->data  = (void*)ms;
511   ms->damping = 0.9;
512   ms->norms   = PETSC_FALSE;
513 
514   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);CHKERRQ(ierr);
515   PetscFunctionReturn(0);
516 }
517 
518