Actual source code: ms.c
1: #include <petsc/private/snesimpl.h>
3: static SNESMSType SNESMSDefault = SNESMSM62;
4: static PetscBool SNESMSRegisterAllCalled;
5: static PetscBool SNESMSPackageInitialized;
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: };
18: typedef struct _SNESMSTableauLink *SNESMSTableauLink;
19: struct _SNESMSTableauLink {
20: struct _SNESMSTableau tab;
21: SNESMSTableauLink next;
22: };
23: static SNESMSTableauLink SNESMSTableauList;
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;
31: /*@C
32: SNESMSRegisterAll - Registers all of the multi-stage methods in SNESMS
34: Not Collective, but should be called by all processes which will need the schemes to be registered
36: Level: advanced
38: .seealso: SNESMSRegisterDestroy()
39: @*/
40: PetscErrorCode SNESMSRegisterAll(void)
41: {
42: if (SNESMSRegisterAllCalled) return 0;
43: SNESMSRegisterAllCalled = PETSC_TRUE;
45: {
46: PetscReal alpha[1] = {1.0};
47: SNESMSRegister(SNESMSEULER,1,1,1.0,NULL,NULL,alpha);
48: }
50: {
51: PetscReal gamma[3][6] = {
52: {0.0000000000000000E+00, -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02, -1.4204296130641869E-01},
53: {1.0000000000000000E+00, 1.1111025767083920E+00, 5.6150921583923230E-01, 7.4151723494934041E-01, 3.1714538168600587E-01, 4.6479276238548706E-01},
54: {0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 6.7968174970583317E-01, -4.1755042846051737E-03, -1.9115668129923846E-01}
55: };
56: PetscReal delta[6] = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
57: PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
58: SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);
59: }
61: { /* Jameson (1983) */
62: PetscReal alpha[4] = {0.25, 0.5, 0.55, 1.0};
63: SNESMSRegister(SNESMSJAMESON83,4,1,1.0,NULL,NULL,alpha);
64: }
66: { /* Van Leer, Tai, and Powell (1989) 1 stage, order 1 */
67: PetscReal alpha[1] = {1.0};
68: SNESMSRegister(SNESMSVLTP11,1,1,0.5,NULL,NULL,alpha);
69: }
70: { /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
71: PetscReal alpha[2] = {0.3333, 1.0};
72: SNESMSRegister(SNESMSVLTP21,2,1,1.0,NULL,NULL,alpha);
73: }
74: { /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
75: PetscReal alpha[3] = {0.1481, 0.4000, 1.0};
76: SNESMSRegister(SNESMSVLTP31,3,1,1.5,NULL,NULL,alpha);
77: }
78: { /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
79: PetscReal alpha[4] = {0.0833, 0.2069, 0.4265, 1.0};
80: SNESMSRegister(SNESMSVLTP41,4,1,2.0,NULL,NULL,alpha);
81: }
82: { /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
83: PetscReal alpha[5] = {0.0533, 0.1263, 0.2375, 0.4414,1.0};
84: SNESMSRegister(SNESMSVLTP51,5,1,2.5,NULL,NULL,alpha);
85: }
86: { /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
87: PetscReal alpha[6] = {0.0370, 0.0851, 0.1521, 0.2562, 0.4512, 1.0};
88: SNESMSRegister(SNESMSVLTP61,6,1,3.0,NULL,NULL,alpha);
89: }
90: return 0;
91: }
93: /*@C
94: SNESMSRegisterDestroy - Frees the list of schemes that were registered by SNESMSRegister().
96: Not Collective
98: Level: advanced
100: .seealso: SNESMSRegister(), SNESMSRegisterAll()
101: @*/
102: PetscErrorCode SNESMSRegisterDestroy(void)
103: {
104: SNESMSTableauLink link;
106: while ((link = SNESMSTableauList)) {
107: SNESMSTableau t = &link->tab;
108: SNESMSTableauList = link->next;
110: PetscFree(t->name);
111: PetscFree(t->gamma);
112: PetscFree(t->delta);
113: PetscFree(t->betasub);
114: PetscFree(link);
115: }
116: SNESMSRegisterAllCalled = PETSC_FALSE;
117: return 0;
118: }
120: /*@C
121: SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
122: from SNESInitializePackage().
124: Level: developer
126: .seealso: PetscInitialize()
127: @*/
128: PetscErrorCode SNESMSInitializePackage(void)
129: {
130: if (SNESMSPackageInitialized) return 0;
131: SNESMSPackageInitialized = PETSC_TRUE;
133: SNESMSRegisterAll();
134: PetscRegisterFinalize(SNESMSFinalizePackage);
135: return 0;
136: }
138: /*@C
139: SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
140: called from PetscFinalize().
142: Level: developer
144: .seealso: PetscFinalize()
145: @*/
146: PetscErrorCode SNESMSFinalizePackage(void)
147: {
148: SNESMSPackageInitialized = PETSC_FALSE;
150: SNESMSRegisterDestroy();
151: return 0;
152: }
154: /*@C
155: SNESMSRegister - register a multistage scheme
157: Not Collective, but the same schemes should be registered on all processes on which they will be used
159: Input Parameters:
160: + name - identifier for method
161: . nstages - number of stages
162: . nregisters - number of registers used by low-storage implementation
163: . stability - scaled stability region
164: . gamma - coefficients, see Ketcheson's paper
165: . delta - coefficients, see Ketcheson's paper
166: - betasub - subdiagonal of Shu-Osher form
168: Notes:
169: The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
171: Many multistage schemes are of the form
172: $ X_0 = X^{(n)}
173: $ X_k = X_0 + \alpha_k * F(X_{k-1}), k = 1,\ldots,s
174: $ X^{(n+1)} = X_s
175: These methods can be registered with
176: .vb
177: SNESMSRegister("name",s,1,stability,NULL,NULL,alpha);
178: .ve
180: Level: advanced
182: .seealso: SNESMS
183: @*/
184: PetscErrorCode SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
185: {
186: SNESMSTableauLink link;
187: SNESMSTableau t;
191: if (gamma || delta) {
195: } else {
197: }
200: SNESMSInitializePackage();
201: PetscNew(&link);
202: t = &link->tab;
203: PetscStrallocpy(name,&t->name);
204: t->nstages = nstages;
205: t->nregisters = nregisters;
206: t->stability = stability;
208: if (gamma && delta) {
209: PetscMalloc1(nstages*nregisters,&t->gamma);
210: PetscMalloc1(nstages,&t->delta);
211: PetscArraycpy(t->gamma,gamma,nstages*nregisters);
212: PetscArraycpy(t->delta,delta,nstages);
213: }
214: PetscMalloc1(nstages,&t->betasub);
215: PetscArraycpy(t->betasub,betasub,nstages);
217: link->next = SNESMSTableauList;
218: SNESMSTableauList = link;
219: return 0;
220: }
222: /*
223: X - initial state, updated in-place.
224: F - residual, computed at the initial X on input
225: */
226: static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
227: {
228: SNES_MS *ms = (SNES_MS*)snes->data;
229: SNESMSTableau t = ms->tableau;
230: const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
231: Vec S1,S2,S3,Y;
232: PetscInt i,nstages = t->nstages;
234: Y = snes->work[0];
235: S1 = X;
236: S2 = snes->work[1];
237: S3 = snes->work[2];
238: VecZeroEntries(S2);
239: VecCopy(X,S3);
240: for (i = 0; i < nstages; i++) {
241: Vec Ss[4];
242: PetscScalar scoeff[4];
244: Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;
246: scoeff[0] = gamma[0*nstages+i] - 1;
247: scoeff[1] = gamma[1*nstages+i];
248: scoeff[2] = gamma[2*nstages+i];
249: scoeff[3] = -betasub[i]*ms->damping;
251: VecAXPY(S2,delta[i],S1);
252: if (i > 0) {
253: SNESComputeFunction(snes,S1,F);
254: }
255: KSPSolve(snes->ksp,F,Y);
256: VecMAXPY(S1,4,scoeff,Ss);
257: }
258: return 0;
259: }
261: /*
262: X - initial state, updated in-place.
263: F - residual, computed at the initial X on input
264: */
265: static PetscErrorCode SNESMSStep_Basic(SNES snes,Vec X,Vec F)
266: {
267: SNES_MS *ms = (SNES_MS*)snes->data;
268: SNESMSTableau tab = ms->tableau;
269: const PetscReal *alpha = tab->betasub, h = ms->damping;
270: PetscInt i,nstages = tab->nstages;
271: Vec X0 = snes->work[0];
273: VecCopy(X,X0);
274: for (i = 0; i < nstages; i++) {
275: if (i > 0) {
276: SNESComputeFunction(snes,X,F);
277: }
278: KSPSolve(snes->ksp,F,X);
279: VecAYPX(X,-alpha[i]*h,X0);
280: }
281: return 0;
282: }
284: static PetscErrorCode SNESMSStep_Step(SNES snes,Vec X,Vec F)
285: {
286: SNES_MS *ms = (SNES_MS*)snes->data;
287: SNESMSTableau tab = ms->tableau;
289: if (tab->gamma && tab->delta) {
290: SNESMSStep_3Sstar(snes,X,F);
291: } else {
292: SNESMSStep_Basic(snes,X,F);
293: }
294: return 0;
295: }
297: static PetscErrorCode SNESMSStep_Norms(SNES snes,PetscInt iter,Vec F)
298: {
299: SNES_MS *ms = (SNES_MS*)snes->data;
300: PetscReal fnorm;
302: if (ms->norms) {
303: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
304: SNESCheckFunctionNorm(snes,fnorm);
305: /* Monitor convergence */
306: PetscObjectSAWsTakeAccess((PetscObject)snes);
307: snes->iter = iter;
308: snes->norm = fnorm;
309: PetscObjectSAWsGrantAccess((PetscObject)snes);
310: SNESLogConvergenceHistory(snes,snes->norm,0);
311: SNESMonitor(snes,snes->iter,snes->norm);
312: /* Test for convergence */
313: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
314: } else if (iter > 0) {
315: PetscObjectSAWsTakeAccess((PetscObject)snes);
316: snes->iter = iter;
317: PetscObjectSAWsGrantAccess((PetscObject)snes);
318: }
319: return 0;
320: }
322: static PetscErrorCode SNESSolve_MS(SNES snes)
323: {
324: SNES_MS *ms = (SNES_MS*)snes->data;
325: Vec X = snes->vec_sol,F = snes->vec_func;
326: PetscInt i;
329: PetscCitationsRegister(SNESCitation,&SNEScite);
331: snes->reason = SNES_CONVERGED_ITERATING;
332: PetscObjectSAWsTakeAccess((PetscObject)snes);
333: snes->iter = 0;
334: snes->norm = 0;
335: PetscObjectSAWsGrantAccess((PetscObject)snes);
337: if (!snes->vec_func_init_set) {
338: SNESComputeFunction(snes,X,F);
339: } else snes->vec_func_init_set = PETSC_FALSE;
341: SNESMSStep_Norms(snes,0,F);
342: if (snes->reason) return 0;
344: for (i = 0; i < snes->max_its; i++) {
346: /* Call general purpose update function */
347: if (snes->ops->update) {
348: (*snes->ops->update)(snes,snes->iter);
349: }
351: if (i == 0 && snes->jacobian) {
352: /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
353: SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);
354: SNESCheckJacobianDomainerror(snes);
355: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
356: }
358: SNESMSStep_Step(snes,X,F);
360: if (i < snes->max_its-1 || ms->norms) {
361: SNESComputeFunction(snes,X,F);
362: }
364: SNESMSStep_Norms(snes,i+1,F);
365: if (snes->reason) return 0;
366: }
367: if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
368: return 0;
369: }
371: static PetscErrorCode SNESSetUp_MS(SNES snes)
372: {
373: SNES_MS *ms = (SNES_MS*)snes->data;
374: SNESMSTableau tab = ms->tableau;
375: PetscInt nwork = tab->nregisters;
377: SNESSetWorkVecs(snes,nwork);
378: SNESSetUpMatrices(snes);
379: return 0;
380: }
382: static PetscErrorCode SNESReset_MS(SNES snes)
383: {
384: return 0;
385: }
387: static PetscErrorCode SNESDestroy_MS(SNES snes)
388: {
389: SNESReset_MS(snes);
390: PetscFree(snes->data);
391: PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetType_C",NULL);
392: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",NULL);
393: PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetDamping_C",NULL);
394: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetDamping_C",NULL);
395: return 0;
396: }
398: static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
399: {
400: PetscBool iascii;
401: SNES_MS *ms = (SNES_MS*)snes->data;
402: SNESMSTableau tab = ms->tableau;
404: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
405: if (iascii) {
406: PetscViewerASCIIPrintf(viewer," multi-stage method type: %s\n",tab->name);
407: }
408: return 0;
409: }
411: static PetscErrorCode SNESSetFromOptions_MS(PetscOptionItems *PetscOptionsObject,SNES snes)
412: {
413: SNES_MS *ms = (SNES_MS*)snes->data;
415: PetscOptionsHead(PetscOptionsObject,"SNES MS options");
416: {
417: SNESMSTableauLink link;
418: PetscInt count,choice;
419: PetscBool flg;
420: const char **namelist;
421: SNESMSType mstype;
422: PetscReal damping;
424: SNESMSGetType(snes,&mstype);
425: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
426: PetscMalloc1(count,(char***)&namelist);
427: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
428: PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);
429: if (flg) SNESMSSetType(snes,namelist[choice]);
430: PetscFree(namelist);
431: SNESMSGetDamping(snes,&damping);
432: PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",damping,&damping,&flg);
433: if (flg) SNESMSSetDamping(snes,damping);
434: PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);
435: }
436: PetscOptionsTail();
437: return 0;
438: }
440: static PetscErrorCode SNESMSGetType_MS(SNES snes,SNESMSType *mstype)
441: {
442: SNES_MS *ms = (SNES_MS*)snes->data;
443: SNESMSTableau tab = ms->tableau;
445: *mstype = tab->name;
446: return 0;
447: }
449: static PetscErrorCode SNESMSSetType_MS(SNES snes,SNESMSType mstype)
450: {
451: SNES_MS *ms = (SNES_MS*)snes->data;
452: SNESMSTableauLink link;
453: PetscBool match;
455: if (ms->tableau) {
456: PetscStrcmp(ms->tableau->name,mstype,&match);
457: if (match) return 0;
458: }
459: for (link = SNESMSTableauList; link; link=link->next) {
460: PetscStrcmp(link->tab.name,mstype,&match);
461: if (match) {
462: if (snes->setupcalled) SNESReset_MS(snes);
463: ms->tableau = &link->tab;
464: if (snes->setupcalled) SNESSetUp_MS(snes);
465: return 0;
466: }
467: }
468: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
469: }
471: /*@C
472: SNESMSGetType - Get the type of multistage smoother
474: Not collective
476: Input Parameter:
477: . snes - nonlinear solver context
479: Output Parameter:
480: . mstype - type of multistage method
482: Level: beginner
484: .seealso: SNESMSSetType(), SNESMSType, SNESMS
485: @*/
486: PetscErrorCode SNESMSGetType(SNES snes,SNESMSType *mstype)
487: {
490: PetscUseMethod(snes,"SNESMSGetType_C",(SNES,SNESMSType*),(snes,mstype));
491: return 0;
492: }
494: /*@C
495: SNESMSSetType - Set the type of multistage smoother
497: Logically collective
499: Input Parameters:
500: + snes - nonlinear solver context
501: - mstype - type of multistage method
503: Level: beginner
505: .seealso: SNESMSGetType(), SNESMSType, SNESMS
506: @*/
507: PetscErrorCode SNESMSSetType(SNES snes,SNESMSType mstype)
508: {
511: PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,mstype));
512: return 0;
513: }
515: static PetscErrorCode SNESMSGetDamping_MS(SNES snes,PetscReal *damping)
516: {
517: SNES_MS *ms = (SNES_MS*)snes->data;
519: *damping = ms->damping;
520: return 0;
521: }
523: static PetscErrorCode SNESMSSetDamping_MS(SNES snes,PetscReal damping)
524: {
525: SNES_MS *ms = (SNES_MS*)snes->data;
527: ms->damping = damping;
528: return 0;
529: }
531: /*@C
532: SNESMSGetDamping - Get the damping parameter
534: Not collective
536: Input Parameter:
537: . snes - nonlinear solver context
539: Output Parameter:
540: . damping - damping parameter
542: Level: advanced
544: .seealso: SNESMSSetDamping(), SNESMS
545: @*/
546: PetscErrorCode SNESMSGetDamping(SNES snes,PetscReal *damping)
547: {
550: PetscUseMethod(snes,"SNESMSGetDamping_C",(SNES,PetscReal*),(snes,damping));
551: return 0;
552: }
554: /*@C
555: SNESMSSetDamping - Set the damping parameter
557: Logically collective
559: Input Parameters:
560: + snes - nonlinear solver context
561: - damping - damping parameter
563: Level: advanced
565: .seealso: SNESMSGetDamping(), SNESMS
566: @*/
567: PetscErrorCode SNESMSSetDamping(SNES snes,PetscReal damping)
568: {
571: PetscTryMethod(snes,"SNESMSSetDamping_C",(SNES,PetscReal),(snes,damping));
572: return 0;
573: }
575: /* -------------------------------------------------------------------------- */
576: /*MC
577: SNESMS - multi-stage smoothers
579: Options Database:
581: + -snes_ms_type - type of multi-stage smoother
582: - -snes_ms_damping - damping for multi-stage method
584: Notes:
585: These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
586: FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
588: Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
590: The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
592: References:
593: + * - Ketcheson (2010) Runge Kutta methods with minimum storage implementations (https://doi.org/10.1016/j.jcp.2009.11.006).
594: . * - 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).
595: . * - Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes (https://doi.org/10.1006/jcph.1997.5772).
596: - * - 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).
598: Level: beginner
600: .seealso: SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
602: M*/
603: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
604: {
605: SNES_MS *ms;
607: SNESMSInitializePackage();
609: snes->ops->setup = SNESSetUp_MS;
610: snes->ops->solve = SNESSolve_MS;
611: snes->ops->destroy = SNESDestroy_MS;
612: snes->ops->setfromoptions = SNESSetFromOptions_MS;
613: snes->ops->view = SNESView_MS;
614: snes->ops->reset = SNESReset_MS;
616: snes->usesnpc = PETSC_FALSE;
617: snes->usesksp = PETSC_TRUE;
619: snes->alwayscomputesfinalresidual = PETSC_FALSE;
621: PetscNewLog(snes,&ms);
622: snes->data = (void*)ms;
623: ms->damping = 0.9;
624: ms->norms = PETSC_FALSE;
626: PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetType_C",SNESMSGetType_MS);
627: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);
628: PetscObjectComposeFunction((PetscObject)snes,"SNESMSGetDamping_C",SNESMSGetDamping_MS);
629: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetDamping_C",SNESMSSetDamping_MS);
631: SNESMSSetType(snes,SNESMSDefault);
632: return 0;
633: }