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 @*/
SNESMSRegisterAll(void)40 PetscErrorCode SNESMSRegisterAll(void)
41 {
42 PetscErrorCode ierr;
43
44 PetscFunctionBegin;
45 if (SNESMSRegisterAllCalled) PetscFunctionReturn(0);
46 SNESMSRegisterAllCalled = PETSC_TRUE;
47
48 {
49 PetscReal alpha[1] = {1.0};
50 ierr = SNESMSRegister(SNESMSEULER,1,1,1.0,NULL,NULL,alpha);CHKERRQ(ierr);
51 }
52
53 {
54 PetscReal gamma[3][6] = {
55 {0.0000000000000000E+00, -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02, -1.4204296130641869E-01},
56 {1.0000000000000000E+00, 1.1111025767083920E+00, 5.6150921583923230E-01, 7.4151723494934041E-01, 3.1714538168600587E-01, 4.6479276238548706E-01},
57 {0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 6.7968174970583317E-01, -4.1755042846051737E-03, -1.9115668129923846E-01}
58 };
59 PetscReal delta[6] = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
60 PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
61 ierr = SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr);
62 }
63
64 { /* Jameson (1983) */
65 PetscReal alpha[4] = {0.25, 0.5, 0.55, 1.0};
66 ierr = SNESMSRegister(SNESMSJAMESON83,4,1,1.0,NULL,NULL,alpha);CHKERRQ(ierr);
67 }
68
69 { /* Van Leer, Tai, and Powell (1989) 1 stage, order 1 */
70 PetscReal alpha[1] = {1.0};
71 ierr = SNESMSRegister(SNESMSVLTP11,1,1,0.5,NULL,NULL,alpha);CHKERRQ(ierr);
72 }
73 { /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
74 PetscReal alpha[2] = {0.3333, 1.0};
75 ierr = SNESMSRegister(SNESMSVLTP21,2,1,1.0,NULL,NULL,alpha);CHKERRQ(ierr);
76 }
77 { /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
78 PetscReal alpha[3] = {0.1481, 0.4000, 1.0};
79 ierr = SNESMSRegister(SNESMSVLTP31,3,1,1.5,NULL,NULL,alpha);CHKERRQ(ierr);
80 }
81 { /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
82 PetscReal alpha[4] = {0.0833, 0.2069, 0.4265, 1.0};
83 ierr = SNESMSRegister(SNESMSVLTP41,4,1,2.0,NULL,NULL,alpha);CHKERRQ(ierr);
84 }
85 { /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
86 PetscReal alpha[5] = {0.0533, 0.1263, 0.2375, 0.4414,1.0};
87 ierr = SNESMSRegister(SNESMSVLTP51,5,1,2.5,NULL,NULL,alpha);CHKERRQ(ierr);
88 }
89 { /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
90 PetscReal alpha[6] = {0.0370, 0.0851, 0.1521, 0.2562, 0.4512, 1.0};
91 ierr = SNESMSRegister(SNESMSVLTP61,6,1,3.0,NULL,NULL,alpha);CHKERRQ(ierr);
92 }
93 PetscFunctionReturn(0);
94 }
95
96 /*@C
97 SNESMSRegisterDestroy - Frees the list of schemes that were registered by SNESMSRegister().
98
99 Not Collective
100
101 Level: advanced
102
103 .seealso: SNESMSRegister(), SNESMSRegisterAll()
104 @*/
SNESMSRegisterDestroy(void)105 PetscErrorCode SNESMSRegisterDestroy(void)
106 {
107 PetscErrorCode ierr;
108 SNESMSTableauLink link;
109
110 PetscFunctionBegin;
111 while ((link = SNESMSTableauList)) {
112 SNESMSTableau t = &link->tab;
113 SNESMSTableauList = link->next;
114
115 ierr = PetscFree(t->name);CHKERRQ(ierr);
116 ierr = PetscFree(t->gamma);CHKERRQ(ierr);
117 ierr = PetscFree(t->delta);CHKERRQ(ierr);
118 ierr = PetscFree(t->betasub);CHKERRQ(ierr);
119 ierr = PetscFree(link);CHKERRQ(ierr);
120 }
121 SNESMSRegisterAllCalled = PETSC_FALSE;
122 PetscFunctionReturn(0);
123 }
124
125 /*@C
126 SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
127 from SNESInitializePackage().
128
129 Level: developer
130
131 .seealso: PetscInitialize()
132 @*/
SNESMSInitializePackage(void)133 PetscErrorCode SNESMSInitializePackage(void)
134 {
135 PetscErrorCode ierr;
136
137 PetscFunctionBegin;
138 if (SNESMSPackageInitialized) PetscFunctionReturn(0);
139 SNESMSPackageInitialized = PETSC_TRUE;
140
141 ierr = SNESMSRegisterAll();CHKERRQ(ierr);
142 ierr = PetscRegisterFinalize(SNESMSFinalizePackage);CHKERRQ(ierr);
143 PetscFunctionReturn(0);
144 }
145
146 /*@C
147 SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
148 called from PetscFinalize().
149
150 Level: developer
151
152 .seealso: PetscFinalize()
153 @*/
SNESMSFinalizePackage(void)154 PetscErrorCode SNESMSFinalizePackage(void)
155 {
156 PetscErrorCode ierr;
157
158 PetscFunctionBegin;
159 SNESMSPackageInitialized = PETSC_FALSE;
160
161 ierr = SNESMSRegisterDestroy();CHKERRQ(ierr);
162 PetscFunctionReturn(0);
163 }
164
165 /*@C
166 SNESMSRegister - register a multistage scheme
167
168 Not Collective, but the same schemes should be registered on all processes on which they will be used
169
170 Input Parameters:
171 + name - identifier for method
172 . nstages - number of stages
173 . nregisters - number of registers used by low-storage implementation
174 . stability - scaled stability region
175 . gamma - coefficients, see Ketcheson's paper
176 . delta - coefficients, see Ketcheson's paper
177 - betasub - subdiagonal of Shu-Osher form
178
179 Notes:
180 The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
181
182 Many multistage schemes are of the form
183 $ X_0 = X^{(n)}
184 $ X_k = X_0 + \alpha_k * F(X_{k-1}), k = 1,\ldots,s
185 $ X^{(n+1)} = X_s
186 These methods can be registered with
187 .vb
188 SNESMSRegister("name",s,1,stability,NULL,NULL,alpha);
189 .ve
190
191 Level: advanced
192
193 .seealso: SNESMS
194 @*/
SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])195 PetscErrorCode SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
196 {
197 PetscErrorCode ierr;
198 SNESMSTableauLink link;
199 SNESMSTableau t;
200
201 PetscFunctionBegin;
202 PetscValidCharPointer(name,1);
203 if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage");
204 if (gamma || delta) {
205 if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form");
206 PetscValidRealPointer(gamma,4);
207 PetscValidRealPointer(delta,5);
208 } else {
209 if (nregisters != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 1-register form");
210 }
211 PetscValidRealPointer(betasub,6);
212
213 ierr = SNESMSInitializePackage();CHKERRQ(ierr);
214 ierr = PetscNew(&link);CHKERRQ(ierr);
215 t = &link->tab;
216 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
217 t->nstages = nstages;
218 t->nregisters = nregisters;
219 t->stability = stability;
220
221 if (gamma && delta) {
222 ierr = PetscMalloc1(nstages*nregisters,&t->gamma);CHKERRQ(ierr);
223 ierr = PetscMalloc1(nstages,&t->delta);CHKERRQ(ierr);
224 ierr = PetscArraycpy(t->gamma,gamma,nstages*nregisters);CHKERRQ(ierr);
225 ierr = PetscArraycpy(t->delta,delta,nstages);CHKERRQ(ierr);
226 }
227 ierr = PetscMalloc1(nstages,&t->betasub);CHKERRQ(ierr);
228 ierr = PetscArraycpy(t->betasub,betasub,nstages);CHKERRQ(ierr);
229
230 link->next = SNESMSTableauList;
231 SNESMSTableauList = link;
232 PetscFunctionReturn(0);
233 }
234
235 /*
236 X - initial state, updated in-place.
237 F - residual, computed at the initial X on input
238 */
SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)239 static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
240 {
241 PetscErrorCode ierr;
242 SNES_MS *ms = (SNES_MS*)snes->data;
243 SNESMSTableau t = ms->tableau;
244 const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
245 Vec S1,S2,S3,Y;
246 PetscInt i,nstages = t->nstages;
247
248 PetscFunctionBegin;
249 Y = snes->work[0];
250 S1 = X;
251 S2 = snes->work[1];
252 S3 = snes->work[2];
253 ierr = VecZeroEntries(S2);CHKERRQ(ierr);
254 ierr = VecCopy(X,S3);CHKERRQ(ierr);
255 for (i = 0; i < nstages; i++) {
256 Vec Ss[4];
257 PetscScalar scoeff[4];
258
259 Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;
260
261 scoeff[0] = gamma[0*nstages+i] - 1;
262 scoeff[1] = gamma[1*nstages+i];
263 scoeff[2] = gamma[2*nstages+i];
264 scoeff[3] = -betasub[i]*ms->damping;
265
266 ierr = VecAXPY(S2,delta[i],S1);CHKERRQ(ierr);
267 if (i > 0) {
268 ierr = SNESComputeFunction(snes,S1,F);CHKERRQ(ierr);
269 }
270 ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr);
271 ierr = VecMAXPY(S1,4,scoeff,Ss);CHKERRQ(ierr);
272 }
273 PetscFunctionReturn(0);
274 }
275
276 /*
277 X - initial state, updated in-place.
278 F - residual, computed at the initial X on input
279 */
SNESMSStep_Basic(SNES snes,Vec X,Vec F)280 static PetscErrorCode SNESMSStep_Basic(SNES snes,Vec X,Vec F)
281 {
282 PetscErrorCode ierr;
283 SNES_MS *ms = (SNES_MS*)snes->data;
284 SNESMSTableau tab = ms->tableau;
285 const PetscReal *alpha = tab->betasub, h = ms->damping;
286 PetscInt i,nstages = tab->nstages;
287 Vec X0 = snes->work[0];
288
289 PetscFunctionBegin;
290 ierr = VecCopy(X,X0);CHKERRQ(ierr);
291 for (i = 0; i < nstages; i++) {
292 if (i > 0) {
293 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
294 }
295 ierr = KSPSolve(snes->ksp,F,X);CHKERRQ(ierr);
296 ierr = VecAYPX(X,-alpha[i]*h,X0);CHKERRQ(ierr);
297 }
298 PetscFunctionReturn(0);
299 }
300
SNESMSStep_Step(SNES snes,Vec X,Vec F)301 static PetscErrorCode SNESMSStep_Step(SNES snes,Vec X,Vec F)
302 {
303 SNES_MS *ms = (SNES_MS*)snes->data;
304 SNESMSTableau tab = ms->tableau;
305 PetscErrorCode ierr;
306
307 PetscFunctionBegin;
308 if (tab->gamma && tab->delta) {
309 ierr = SNESMSStep_3Sstar(snes,X,F);CHKERRQ(ierr);
310 } else {
311 ierr = SNESMSStep_Basic(snes,X,F);CHKERRQ(ierr);
312 }
313 PetscFunctionReturn(0);
314 }
315
SNESMSStep_Norms(SNES snes,PetscInt iter,Vec F)316 static PetscErrorCode SNESMSStep_Norms(SNES snes,PetscInt iter,Vec F)
317 {
318 SNES_MS *ms = (SNES_MS*)snes->data;
319 PetscReal fnorm;
320 PetscErrorCode ierr;
321
322 PetscFunctionBegin;
323 if (ms->norms) {
324 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */
325 SNESCheckFunctionNorm(snes,fnorm);
326 /* Monitor convergence */
327 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
328 snes->iter = iter;
329 snes->norm = fnorm;
330 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
331 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
332 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
333 /* Test for convergence */
334 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
335 } else if (iter > 0) {
336 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
337 snes->iter = iter;
338 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
339 }
340 PetscFunctionReturn(0);
341 }
342
SNESSolve_MS(SNES snes)343 static PetscErrorCode SNESSolve_MS(SNES snes)
344 {
345 SNES_MS *ms = (SNES_MS*)snes->data;
346 Vec X = snes->vec_sol,F = snes->vec_func;
347 PetscInt i;
348 PetscErrorCode ierr;
349
350 PetscFunctionBegin;
351 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);
352 ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
353
354 snes->reason = SNES_CONVERGED_ITERATING;
355 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
356 snes->iter = 0;
357 snes->norm = 0;
358 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
359
360 if (!snes->vec_func_init_set) {
361 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
362 } else snes->vec_func_init_set = PETSC_FALSE;
363
364 ierr = SNESMSStep_Norms(snes,0,F);CHKERRQ(ierr);
365 if (snes->reason) PetscFunctionReturn(0);
366
367 for (i = 0; i < snes->max_its; i++) {
368
369 /* Call general purpose update function */
370 if (snes->ops->update) {
371 ierr = (*snes->ops->update)(snes,snes->iter);CHKERRQ(ierr);
372 }
373
374 if (i == 0 && snes->jacobian) {
375 /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
376 ierr = SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
377 SNESCheckJacobianDomainerror(snes);
378 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
379 }
380
381 ierr = SNESMSStep_Step(snes,X,F);CHKERRQ(ierr);
382
383 if (i < snes->max_its-1 || ms->norms) {
384 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
385 }
386
387 ierr = SNESMSStep_Norms(snes,i+1,F);CHKERRQ(ierr);
388 if (snes->reason) PetscFunctionReturn(0);
389 }
390 if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
391 PetscFunctionReturn(0);
392 }
393
SNESSetUp_MS(SNES snes)394 static PetscErrorCode SNESSetUp_MS(SNES snes)
395 {
396 SNES_MS *ms = (SNES_MS*)snes->data;
397 SNESMSTableau tab = ms->tableau;
398 PetscInt nwork = tab->nregisters;
399 PetscErrorCode ierr;
400
401 PetscFunctionBegin;
402 ierr = SNESSetWorkVecs(snes,nwork);CHKERRQ(ierr);
403 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
404 PetscFunctionReturn(0);
405 }
406
SNESReset_MS(SNES snes)407 static PetscErrorCode SNESReset_MS(SNES snes)
408 {
409 PetscFunctionBegin;
410 PetscFunctionReturn(0);
411 }
412
SNESDestroy_MS(SNES snes)413 static PetscErrorCode SNESDestroy_MS(SNES snes)
414 {
415 PetscErrorCode ierr;
416
417 PetscFunctionBegin;
418 ierr = SNESReset_MS(snes);CHKERRQ(ierr);
419 ierr = PetscFree(snes->data);CHKERRQ(ierr);
420 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetType_C",NULL);CHKERRQ(ierr);
421 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",NULL);CHKERRQ(ierr);
422 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetDamping_C",NULL);CHKERRQ(ierr);
423 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetDamping_C",NULL);CHKERRQ(ierr);
424 PetscFunctionReturn(0);
425 }
426
SNESView_MS(SNES snes,PetscViewer viewer)427 static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
428 {
429 PetscBool iascii;
430 PetscErrorCode ierr;
431 SNES_MS *ms = (SNES_MS*)snes->data;
432 SNESMSTableau tab = ms->tableau;
433
434 PetscFunctionBegin;
435 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
436 if (iascii) {
437 ierr = PetscViewerASCIIPrintf(viewer," multi-stage method type: %s\n",tab->name);CHKERRQ(ierr);
438 }
439 PetscFunctionReturn(0);
440 }
441
SNESSetFromOptions_MS(PetscOptionItems * PetscOptionsObject,SNES snes)442 static PetscErrorCode SNESSetFromOptions_MS(PetscOptionItems *PetscOptionsObject,SNES snes)
443 {
444 SNES_MS *ms = (SNES_MS*)snes->data;
445 PetscErrorCode ierr;
446
447 PetscFunctionBegin;
448 ierr = PetscOptionsHead(PetscOptionsObject,"SNES MS options");CHKERRQ(ierr);
449 {
450 SNESMSTableauLink link;
451 PetscInt count,choice;
452 PetscBool flg;
453 const char **namelist;
454 SNESMSType mstype;
455 PetscReal damping;
456
457 ierr = SNESMSGetType(snes,&mstype);CHKERRQ(ierr);
458 for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
459 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr);
460 for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
461 ierr = PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);CHKERRQ(ierr);
462 if (flg) {ierr = SNESMSSetType(snes,namelist[choice]);CHKERRQ(ierr);}
463 ierr = PetscFree(namelist);CHKERRQ(ierr);
464 ierr = SNESMSGetDamping(snes,&damping);CHKERRQ(ierr);
465 ierr = PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",damping,&damping,&flg);CHKERRQ(ierr);
466 if (flg) {ierr = SNESMSSetDamping(snes,damping);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
SNESMSGetType_MS(SNES snes,SNESMSType * mstype)473 static PetscErrorCode SNESMSGetType_MS(SNES snes,SNESMSType *mstype)
474 {
475 SNES_MS *ms = (SNES_MS*)snes->data;
476 SNESMSTableau tab = ms->tableau;
477
478 PetscFunctionBegin;
479 *mstype = tab->name;
480 PetscFunctionReturn(0);
481 }
482
SNESMSSetType_MS(SNES snes,SNESMSType mstype)483 static PetscErrorCode SNESMSSetType_MS(SNES snes,SNESMSType mstype)
484 {
485 SNES_MS *ms = (SNES_MS*)snes->data;
486 SNESMSTableauLink link;
487 PetscBool match;
488 PetscErrorCode ierr;
489
490 PetscFunctionBegin;
491 if (ms->tableau) {
492 ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
493 if (match) PetscFunctionReturn(0);
494 }
495 for (link = SNESMSTableauList; link; link=link->next) {
496 ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
497 if (match) {
498 if (snes->setupcalled) {ierr = SNESReset_MS(snes);CHKERRQ(ierr);}
499 ms->tableau = &link->tab;
500 if (snes->setupcalled) {ierr = SNESSetUp_MS(snes);CHKERRQ(ierr);}
501 PetscFunctionReturn(0);
502 }
503 }
504 SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
505 PetscFunctionReturn(0);
506 }
507
508 /*@C
509 SNESMSGetType - Get the type of multistage smoother
510
511 Not collective
512
513 Input Parameter:
514 . snes - nonlinear solver context
515
516 Output Parameter:
517 . mstype - type of multistage method
518
519 Level: beginner
520
521 .seealso: SNESMSSetType(), SNESMSType, SNESMS
522 @*/
SNESMSGetType(SNES snes,SNESMSType * mstype)523 PetscErrorCode SNESMSGetType(SNES snes,SNESMSType *mstype)
524 {
525 PetscErrorCode ierr;
526
527 PetscFunctionBegin;
528 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
529 PetscValidPointer(mstype,2);
530 ierr = PetscUseMethod(snes,"SNESMSGetType_C",(SNES,SNESMSType*),(snes,mstype));CHKERRQ(ierr);
531 PetscFunctionReturn(0);
532 }
533
534 /*@C
535 SNESMSSetType - Set the type of multistage smoother
536
537 Logically collective
538
539 Input Parameter:
540 + snes - nonlinear solver context
541 - mstype - type of multistage method
542
543 Level: beginner
544
545 .seealso: SNESMSGetType(), SNESMSType, SNESMS
546 @*/
SNESMSSetType(SNES snes,SNESMSType mstype)547 PetscErrorCode SNESMSSetType(SNES snes,SNESMSType mstype)
548 {
549 PetscErrorCode ierr;
550
551 PetscFunctionBegin;
552 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
553 PetscValidPointer(mstype,2);
554 ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,mstype));CHKERRQ(ierr);
555 PetscFunctionReturn(0);
556 }
557
SNESMSGetDamping_MS(SNES snes,PetscReal * damping)558 static PetscErrorCode SNESMSGetDamping_MS(SNES snes,PetscReal *damping)
559 {
560 SNES_MS *ms = (SNES_MS*)snes->data;
561
562 PetscFunctionBegin;
563 *damping = ms->damping;
564 PetscFunctionReturn(0);
565 }
566
SNESMSSetDamping_MS(SNES snes,PetscReal damping)567 static PetscErrorCode SNESMSSetDamping_MS(SNES snes,PetscReal damping)
568 {
569 SNES_MS *ms = (SNES_MS*)snes->data;
570
571 PetscFunctionBegin;
572 ms->damping = damping;
573 PetscFunctionReturn(0);
574 }
575
576 /*@C
577 SNESMSGetDamping - Get the damping parameter
578
579 Not collective
580
581 Input Parameter:
582 . snes - nonlinear solver context
583
584 Output Parameter:
585 . damping - damping parameter
586
587 Level: advanced
588
589 .seealso: SNESMSSetDamping(), SNESMS
590 @*/
SNESMSGetDamping(SNES snes,PetscReal * damping)591 PetscErrorCode SNESMSGetDamping(SNES snes,PetscReal *damping)
592 {
593 PetscErrorCode ierr;
594
595 PetscFunctionBegin;
596 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
597 PetscValidPointer(damping,2);
598 ierr = PetscUseMethod(snes,"SNESMSGetDamping_C",(SNES,PetscReal*),(snes,damping));CHKERRQ(ierr);
599 PetscFunctionReturn(0);
600 }
601
602 /*@C
603 SNESMSSetDamping - Set the damping parameter
604
605 Logically collective
606
607 Input Parameter:
608 + snes - nonlinear solver context
609 - damping - damping parameter
610
611 Level: advanced
612
613 .seealso: SNESMSGetDamping(), SNESMS
614 @*/
SNESMSSetDamping(SNES snes,PetscReal damping)615 PetscErrorCode SNESMSSetDamping(SNES snes,PetscReal damping)
616 {
617 PetscErrorCode ierr;
618
619 PetscFunctionBegin;
620 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
621 PetscValidLogicalCollectiveReal(snes,damping,2);
622 ierr = PetscTryMethod(snes,"SNESMSSetDamping_C",(SNES,PetscReal),(snes,damping));CHKERRQ(ierr);
623 PetscFunctionReturn(0);
624 }
625
626 /* -------------------------------------------------------------------------- */
627 /*MC
628 SNESMS - multi-stage smoothers
629
630 Options Database:
631
632 + -snes_ms_type - type of multi-stage smoother
633 - -snes_ms_damping - damping for multi-stage method
634
635 Notes:
636 These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
637 FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
638
639 Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
640
641 The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
642
643 References:
644 + 1. - Ketcheson (2010) Runge Kutta methods with minimum storage implementations (https://doi.org/10.1016/j.jcp.2009.11.006).
645 . 2. - Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method (https://doi.org/10.1016/0096-3003(83)90019-X).
646 . 3. - Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes (https://doi.org/10.1006/jcph.1997.5772).
647 - 4. - Van Leer, Tai, and Powell (1989) Design of optimally smoothing multi-stage schemes for the Euler equations (https://doi.org/10.2514/6.1989-1933).
648
649 Level: beginner
650
651 .seealso: SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
652
653 M*/
SNESCreate_MS(SNES snes)654 PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
655 {
656 PetscErrorCode ierr;
657 SNES_MS *ms;
658
659 PetscFunctionBegin;
660 ierr = SNESMSInitializePackage();CHKERRQ(ierr);
661
662 snes->ops->setup = SNESSetUp_MS;
663 snes->ops->solve = SNESSolve_MS;
664 snes->ops->destroy = SNESDestroy_MS;
665 snes->ops->setfromoptions = SNESSetFromOptions_MS;
666 snes->ops->view = SNESView_MS;
667 snes->ops->reset = SNESReset_MS;
668
669 snes->usesnpc = PETSC_FALSE;
670 snes->usesksp = PETSC_TRUE;
671
672 snes->alwayscomputesfinalresidual = PETSC_FALSE;
673
674 ierr = PetscNewLog(snes,&ms);CHKERRQ(ierr);
675 snes->data = (void*)ms;
676 ms->damping = 0.9;
677 ms->norms = PETSC_FALSE;
678
679 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetType_C",SNESMSGetType_MS);CHKERRQ(ierr);
680 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);CHKERRQ(ierr);
681 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetDamping_C",SNESMSGetDamping_MS);CHKERRQ(ierr);
682 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetDamping_C",SNESMSSetDamping_MS);CHKERRQ(ierr);
683
684 ierr = SNESMSSetType(snes,SNESMSDefault);CHKERRQ(ierr);
685 PetscFunctionReturn(0);
686 }
687