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: Logically Collective
36: Level: developer
38: .seealso: `SNES`, `SNESMS`, `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: developer
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 for `SNESMS`
157: Logically Collective
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 PetscInt nstages = t->nstages;
231: const PetscReal *gamma = t->gamma, *delta = t->delta, *betasub = t->betasub;
232: Vec S1 = X, S2 = snes->work[1], S3 = snes->work[2], Y = snes->work[0], S1copy = snes->work[3];
234: VecZeroEntries(S2);
235: VecCopy(X, S3);
236: for (PetscInt i = 0; i < nstages; i++) {
237: Vec Ss[] = {S1copy, S2, S3, Y};
238: const PetscScalar scoeff[] = {gamma[0 * nstages + i] - 1, gamma[1 * nstages + i], gamma[2 * nstages + i], -betasub[i] * ms->damping};
240: VecAXPY(S2, delta[i], S1);
241: if (i > 0) SNESComputeFunction(snes, S1, F);
242: KSPSolve(snes->ksp, F, Y);
243: VecCopy(S1, S1copy);
244: VecMAXPY(S1, 4, scoeff, Ss);
245: }
246: return 0;
247: }
249: /*
250: X - initial state, updated in-place.
251: F - residual, computed at the initial X on input
252: */
253: static PetscErrorCode SNESMSStep_Basic(SNES snes, Vec X, Vec F)
254: {
255: SNES_MS *ms = (SNES_MS *)snes->data;
256: SNESMSTableau tab = ms->tableau;
257: const PetscReal *alpha = tab->betasub, h = ms->damping;
258: PetscInt i, nstages = tab->nstages;
259: Vec X0 = snes->work[0];
261: VecCopy(X, X0);
262: for (i = 0; i < nstages; i++) {
263: if (i > 0) SNESComputeFunction(snes, X, F);
264: KSPSolve(snes->ksp, F, X);
265: VecAYPX(X, -alpha[i] * h, X0);
266: }
267: return 0;
268: }
270: static PetscErrorCode SNESMSStep_Step(SNES snes, Vec X, Vec F)
271: {
272: SNES_MS *ms = (SNES_MS *)snes->data;
273: SNESMSTableau tab = ms->tableau;
275: if (tab->gamma && tab->delta) {
276: SNESMSStep_3Sstar(snes, X, F);
277: } else {
278: SNESMSStep_Basic(snes, X, F);
279: }
280: return 0;
281: }
283: static PetscErrorCode SNESMSStep_Norms(SNES snes, PetscInt iter, Vec F)
284: {
285: SNES_MS *ms = (SNES_MS *)snes->data;
286: PetscReal fnorm;
288: if (ms->norms) {
289: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
290: SNESCheckFunctionNorm(snes, fnorm);
291: /* Monitor convergence */
292: PetscObjectSAWsTakeAccess((PetscObject)snes);
293: snes->iter = iter;
294: snes->norm = fnorm;
295: PetscObjectSAWsGrantAccess((PetscObject)snes);
296: SNESLogConvergenceHistory(snes, snes->norm, 0);
297: SNESMonitor(snes, snes->iter, snes->norm);
298: /* Test for convergence */
299: PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
300: } else if (iter > 0) {
301: PetscObjectSAWsTakeAccess((PetscObject)snes);
302: snes->iter = iter;
303: PetscObjectSAWsGrantAccess((PetscObject)snes);
304: }
305: return 0;
306: }
308: static PetscErrorCode SNESSolve_MS(SNES snes)
309: {
310: SNES_MS *ms = (SNES_MS *)snes->data;
311: Vec X = snes->vec_sol, F = snes->vec_func;
312: PetscInt i;
315: PetscCitationsRegister(SNESCitation, &SNEScite);
317: snes->reason = SNES_CONVERGED_ITERATING;
318: PetscObjectSAWsTakeAccess((PetscObject)snes);
319: snes->iter = 0;
320: snes->norm = 0;
321: PetscObjectSAWsGrantAccess((PetscObject)snes);
323: if (!snes->vec_func_init_set) {
324: SNESComputeFunction(snes, X, F);
325: } else snes->vec_func_init_set = PETSC_FALSE;
327: SNESMSStep_Norms(snes, 0, F);
328: if (snes->reason) return 0;
330: for (i = 0; i < snes->max_its; i++) {
331: /* Call general purpose update function */
332: PetscTryTypeMethod(snes, update, snes->iter);
334: if (i == 0 && snes->jacobian) {
335: /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
336: SNESComputeJacobian(snes, snes->vec_sol, snes->jacobian, snes->jacobian_pre);
337: SNESCheckJacobianDomainerror(snes);
338: KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre);
339: }
341: SNESMSStep_Step(snes, X, F);
343: if (i < snes->max_its - 1 || ms->norms) SNESComputeFunction(snes, X, F);
345: SNESMSStep_Norms(snes, i + 1, F);
346: if (snes->reason) return 0;
347: }
348: if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
349: return 0;
350: }
352: static PetscErrorCode SNESSetUp_MS(SNES snes)
353: {
354: SNES_MS *ms = (SNES_MS *)snes->data;
355: SNESMSTableau tab = ms->tableau;
356: PetscInt nwork = tab->nregisters + 1; // +1 because VecMAXPY() in SNESMSStep_3Sstar()
357: // needs an extra work vec
359: SNESSetWorkVecs(snes, nwork);
360: SNESSetUpMatrices(snes);
361: return 0;
362: }
364: static PetscErrorCode SNESReset_MS(SNES snes)
365: {
366: return 0;
367: }
369: static PetscErrorCode SNESDestroy_MS(SNES snes)
370: {
371: SNESReset_MS(snes);
372: PetscFree(snes->data);
373: PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", NULL);
374: PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", NULL);
375: PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", NULL);
376: PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", NULL);
377: return 0;
378: }
380: static PetscErrorCode SNESView_MS(SNES snes, PetscViewer viewer)
381: {
382: PetscBool iascii;
383: SNES_MS *ms = (SNES_MS *)snes->data;
384: SNESMSTableau tab = ms->tableau;
386: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
387: if (iascii) PetscViewerASCIIPrintf(viewer, " multi-stage method type: %s\n", tab->name);
388: return 0;
389: }
391: static PetscErrorCode SNESSetFromOptions_MS(SNES snes, PetscOptionItems *PetscOptionsObject)
392: {
393: SNES_MS *ms = (SNES_MS *)snes->data;
395: PetscOptionsHeadBegin(PetscOptionsObject, "SNES MS options");
396: {
397: SNESMSTableauLink link;
398: PetscInt count, choice;
399: PetscBool flg;
400: const char **namelist;
401: SNESMSType mstype;
402: PetscReal damping;
404: SNESMSGetType(snes, &mstype);
405: for (link = SNESMSTableauList, count = 0; link; link = link->next, count++)
406: ;
407: PetscMalloc1(count, (char ***)&namelist);
408: for (link = SNESMSTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
409: PetscOptionsEList("-snes_ms_type", "Multistage smoother type", "SNESMSSetType", (const char *const *)namelist, count, mstype, &choice, &flg);
410: if (flg) SNESMSSetType(snes, namelist[choice]);
411: PetscFree(namelist);
412: SNESMSGetDamping(snes, &damping);
413: PetscOptionsReal("-snes_ms_damping", "Damping for multistage method", "SNESMSSetDamping", damping, &damping, &flg);
414: if (flg) SNESMSSetDamping(snes, damping);
415: PetscOptionsBool("-snes_ms_norms", "Compute norms for monitoring", "none", ms->norms, &ms->norms, NULL);
416: }
417: PetscOptionsHeadEnd();
418: return 0;
419: }
421: static PetscErrorCode SNESMSGetType_MS(SNES snes, SNESMSType *mstype)
422: {
423: SNES_MS *ms = (SNES_MS *)snes->data;
424: SNESMSTableau tab = ms->tableau;
426: *mstype = tab->name;
427: return 0;
428: }
430: static PetscErrorCode SNESMSSetType_MS(SNES snes, SNESMSType mstype)
431: {
432: SNES_MS *ms = (SNES_MS *)snes->data;
433: SNESMSTableauLink link;
434: PetscBool match;
436: if (ms->tableau) {
437: PetscStrcmp(ms->tableau->name, mstype, &match);
438: if (match) return 0;
439: }
440: for (link = SNESMSTableauList; link; link = link->next) {
441: PetscStrcmp(link->tab.name, mstype, &match);
442: if (match) {
443: if (snes->setupcalled) SNESReset_MS(snes);
444: ms->tableau = &link->tab;
445: if (snes->setupcalled) SNESSetUp_MS(snes);
446: return 0;
447: }
448: }
449: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mstype);
450: }
452: /*@C
453: SNESMSGetType - Get the type of multistage smoother `SNESMS`
455: Not collective
457: Input Parameter:
458: . snes - nonlinear solver context
460: Output Parameter:
461: . mstype - type of multistage method
463: Level: advanced
465: .seealso: `SNESMS`, `SNESMSSetType()`, `SNESMSType`, `SNESMS`
466: @*/
467: PetscErrorCode SNESMSGetType(SNES snes, SNESMSType *mstype)
468: {
471: PetscUseMethod(snes, "SNESMSGetType_C", (SNES, SNESMSType *), (snes, mstype));
472: return 0;
473: }
475: /*@C
476: SNESMSSetType - Set the type of multistage smoother `SNESMS`
478: Logically collective
480: Input Parameters:
481: + snes - nonlinear solver context
482: - mstype - type of multistage method
484: Level: advanced
486: .seealso: `SNESMS`, `SNESMSGetType()`, `SNESMSType`, `SNESMS`
487: @*/
488: PetscErrorCode SNESMSSetType(SNES snes, SNESMSType mstype)
489: {
492: PetscTryMethod(snes, "SNESMSSetType_C", (SNES, SNESMSType), (snes, mstype));
493: return 0;
494: }
496: static PetscErrorCode SNESMSGetDamping_MS(SNES snes, PetscReal *damping)
497: {
498: SNES_MS *ms = (SNES_MS *)snes->data;
500: *damping = ms->damping;
501: return 0;
502: }
504: static PetscErrorCode SNESMSSetDamping_MS(SNES snes, PetscReal damping)
505: {
506: SNES_MS *ms = (SNES_MS *)snes->data;
508: ms->damping = damping;
509: return 0;
510: }
512: /*@C
513: SNESMSGetDamping - Get the damping parameter of `SNESMS` multistage scheme
515: Not collective
517: Input Parameter:
518: . snes - nonlinear solver context
520: Output Parameter:
521: . damping - damping parameter
523: Level: advanced
525: .seealso: `SNESMSSetDamping()`, `SNESMS`
526: @*/
527: PetscErrorCode SNESMSGetDamping(SNES snes, PetscReal *damping)
528: {
531: PetscUseMethod(snes, "SNESMSGetDamping_C", (SNES, PetscReal *), (snes, damping));
532: return 0;
533: }
535: /*@C
536: SNESMSSetDamping - Set the damping parameter for a `SNESMS` multistage scheme
538: Logically collective
540: Input Parameters:
541: + snes - nonlinear solver context
542: - damping - damping parameter
544: Level: advanced
546: .seealso: `SNESMSGetDamping()`, `SNESMS`
547: @*/
548: PetscErrorCode SNESMSSetDamping(SNES snes, PetscReal damping)
549: {
552: PetscTryMethod(snes, "SNESMSSetDamping_C", (SNES, PetscReal), (snes, damping));
553: return 0;
554: }
556: /*MC
557: SNESMS - multi-stage smoothers
559: Options Database Keys:
560: + -snes_ms_type - type of multi-stage smoother
561: - -snes_ms_damping - damping for multi-stage method
563: Notes:
564: These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
565: FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
567: Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
569: The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with `SNESMSRegister()`.
571: References:
572: + * - Ketcheson (2010) Runge Kutta methods with minimum storage implementations (https://doi.org/10.1016/j.jcp.2009.11.006).
573: . * - 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).
574: . * - Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes (https://doi.org/10.1006/jcph.1997.5772).
575: - * - 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).
577: Level: advanced
579: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESMS`, `SNESFAS`, `KSPCHEBYSHEV`
581: M*/
582: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
583: {
584: SNES_MS *ms;
586: SNESMSInitializePackage();
588: snes->ops->setup = SNESSetUp_MS;
589: snes->ops->solve = SNESSolve_MS;
590: snes->ops->destroy = SNESDestroy_MS;
591: snes->ops->setfromoptions = SNESSetFromOptions_MS;
592: snes->ops->view = SNESView_MS;
593: snes->ops->reset = SNESReset_MS;
595: snes->usesnpc = PETSC_FALSE;
596: snes->usesksp = PETSC_TRUE;
598: snes->alwayscomputesfinalresidual = PETSC_FALSE;
600: PetscNew(&ms);
601: snes->data = (void *)ms;
602: ms->damping = 0.9;
603: ms->norms = PETSC_FALSE;
605: PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", SNESMSGetType_MS);
606: PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", SNESMSSetType_MS);
607: PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", SNESMSGetDamping_MS);
608: PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", SNESMSSetDamping_MS);
610: SNESMSSetType(snes, SNESMSDefault);
611: return 0;
612: }