xref: /petsc/src/snes/impls/ms/ms.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94)
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   PetscValidRealPointer(gamma,4);
202   PetscValidRealPointer(delta,5);
203   PetscValidRealPointer(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] - 1;
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 SNESMSStep_Norms(SNES snes,PetscInt iter,Vec F)
266 {
267   SNES_MS        *ms = (SNES_MS*)snes->data;
268   PetscReal      fnorm;
269   PetscErrorCode ierr;
270 
271   PetscFunctionBegin;
272   if (ms->norms) {
273     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F||  */
274     SNESCheckFunctionNorm(snes,fnorm);
275     /* Monitor convergence */
276     ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
277     snes->iter = iter;
278     snes->norm = fnorm;
279     ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
280     ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
281     ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
282     /* Test for convergence */
283     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
284   } else if (iter > 0) {
285     ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
286     snes->iter = iter;
287     ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
288   }
289   PetscFunctionReturn(0);
290 }
291 
292 
293 static PetscErrorCode SNESSolve_MS(SNES snes)
294 {
295   SNES_MS        *ms = (SNES_MS*)snes->data;
296   Vec            X   = snes->vec_sol,F = snes->vec_func;
297   PetscInt       i;
298   PetscErrorCode ierr;
299 
300   PetscFunctionBegin;
301   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);
302   ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
303 
304   snes->reason = SNES_CONVERGED_ITERATING;
305   ierr         = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
306   snes->iter   = 0;
307   snes->norm   = 0;
308   ierr         = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
309 
310   if (!snes->vec_func_init_set) {
311     ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
312   } else snes->vec_func_init_set = PETSC_FALSE;
313 
314   ierr = SNESMSStep_Norms(snes,0,F);CHKERRQ(ierr);
315   if (snes->reason) PetscFunctionReturn(0);
316 
317   for (i = 0; i < snes->max_its; i++) {
318 
319     /* Call general purpose update function */
320     if (snes->ops->update) {
321       ierr = (*snes->ops->update)(snes,snes->iter);CHKERRQ(ierr);
322     }
323 
324     if (i == 0 && snes->jacobian) {
325       /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
326       ierr = SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
327       SNESCheckJacobianDomainerror(snes);
328       ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
329     }
330 
331     ierr = SNESMSStep_3Sstar(snes,X,F);CHKERRQ(ierr);
332 
333     if (i < snes->max_its-1 || ms->norms) {
334       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
335     }
336 
337     ierr = SNESMSStep_Norms(snes,i+1,F);CHKERRQ(ierr);
338     if (snes->reason) PetscFunctionReturn(0);
339   }
340   if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
341   PetscFunctionReturn(0);
342 }
343 
344 static PetscErrorCode SNESSetUp_MS(SNES snes)
345 {
346   SNES_MS        *ms = (SNES_MS*)snes->data;
347   PetscErrorCode ierr;
348 
349   PetscFunctionBegin;
350   if (!ms->tableau) {ierr = SNESMSSetType(snes,SNESMSDefault);CHKERRQ(ierr);}
351   ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr);
352   ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
353   PetscFunctionReturn(0);
354 }
355 
356 static PetscErrorCode SNESReset_MS(SNES snes)
357 {
358 
359   PetscFunctionBegin;
360   PetscFunctionReturn(0);
361 }
362 
363 static PetscErrorCode SNESDestroy_MS(SNES snes)
364 {
365   PetscErrorCode ierr;
366 
367   PetscFunctionBegin;
368   ierr = SNESReset_MS(snes);CHKERRQ(ierr);
369   ierr = PetscFree(snes->data);CHKERRQ(ierr);
370   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",NULL);CHKERRQ(ierr);
371   PetscFunctionReturn(0);
372 }
373 
374 static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
375 {
376   PetscBool      iascii;
377   PetscErrorCode ierr;
378   SNES_MS        *ms = (SNES_MS*)snes->data;
379 
380   PetscFunctionBegin;
381   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
382   if (iascii) {
383     SNESMSTableau tab = ms->tableau;
384     ierr = PetscViewerASCIIPrintf(viewer,"  multi-stage method type: %s\n",tab ? tab->name : "not yet set");CHKERRQ(ierr);
385   }
386   PetscFunctionReturn(0);
387 }
388 
389 static PetscErrorCode SNESSetFromOptions_MS(PetscOptionItems *PetscOptionsObject,SNES snes)
390 {
391   SNES_MS        *ms = (SNES_MS*)snes->data;
392   PetscErrorCode ierr;
393 
394   PetscFunctionBegin;
395   ierr = PetscOptionsHead(PetscOptionsObject,"SNES MS options");CHKERRQ(ierr);
396   {
397     SNESMSTableauLink link;
398     PetscInt          count,choice;
399     PetscBool         flg;
400     const char        **namelist;
401     char              mstype[256];
402 
403     ierr = PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));CHKERRQ(ierr);
404     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
405     ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr);
406     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
407     ierr = PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);CHKERRQ(ierr);
408     ierr = SNESMSSetType(snes,flg ? namelist[choice] : mstype);CHKERRQ(ierr);
409     ierr = PetscFree(namelist);CHKERRQ(ierr);
410     ierr = PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,NULL);CHKERRQ(ierr);
411     ierr = PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);CHKERRQ(ierr);
412   }
413   ierr = PetscOptionsTail();CHKERRQ(ierr);
414   PetscFunctionReturn(0);
415 }
416 
417 PetscErrorCode  SNESMSSetType_MS(SNES snes,SNESMSType mstype)
418 {
419   PetscErrorCode    ierr;
420   SNES_MS           *ms = (SNES_MS*)snes->data;
421   SNESMSTableauLink link;
422   PetscBool         match;
423 
424   PetscFunctionBegin;
425   if (ms->tableau) {
426     ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
427     if (match) PetscFunctionReturn(0);
428   }
429   for (link = SNESMSTableauList; link; link=link->next) {
430     ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
431     if (match) {
432       ierr        = SNESReset_MS(snes);CHKERRQ(ierr);
433       ms->tableau = &link->tab;
434       PetscFunctionReturn(0);
435     }
436   }
437   SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
438   PetscFunctionReturn(0);
439 }
440 
441 /*@C
442   SNESMSSetType - Set the type of multistage smoother
443 
444   Logically collective
445 
446   Input Parameter:
447 +  snes - nonlinear solver context
448 -  mstype - type of multistage method
449 
450   Level: beginner
451 
452 .seealso: SNESMSGetType(), SNESMS
453 @*/
454 PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
455 {
456   PetscErrorCode ierr;
457 
458   PetscFunctionBegin;
459   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
460   ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));CHKERRQ(ierr);
461   PetscFunctionReturn(0);
462 }
463 
464 /* -------------------------------------------------------------------------- */
465 /*MC
466       SNESMS - multi-stage smoothers
467 
468       Options Database:
469 
470 +     -snes_ms_type - type of multi-stage smoother
471 -     -snes_ms_damping - damping for multi-stage method
472 
473       Notes:
474       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
475       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
476 
477       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
478 
479       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
480 
481       References:
482 +   1. -   Ketcheson (2010) Runge Kutta methods with minimum storage implementations.
483 .   2. -   Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
484 -   3. -   Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
485 
486       Level: beginner
487 
488 .seealso:  SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
489 
490 M*/
491 PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
492 {
493   PetscErrorCode ierr;
494   SNES_MS        *ms;
495 
496   PetscFunctionBegin;
497   ierr = SNESMSInitializePackage();CHKERRQ(ierr);
498 
499   snes->ops->setup          = SNESSetUp_MS;
500   snes->ops->solve          = SNESSolve_MS;
501   snes->ops->destroy        = SNESDestroy_MS;
502   snes->ops->setfromoptions = SNESSetFromOptions_MS;
503   snes->ops->view           = SNESView_MS;
504   snes->ops->reset          = SNESReset_MS;
505 
506   snes->usesnpc = PETSC_FALSE;
507   snes->usesksp = PETSC_TRUE;
508 
509   snes->alwayscomputesfinalresidual = PETSC_FALSE;
510 
511   ierr        = PetscNewLog(snes,&ms);CHKERRQ(ierr);
512   snes->data  = (void*)ms;
513   ms->damping = 0.9;
514   ms->norms   = PETSC_FALSE;
515 
516   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);CHKERRQ(ierr);
517   PetscFunctionReturn(0);
518 }
519 
520