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: } SNES_MS;
30: /*@C
31: SNESMSRegisterAll - Registers all of the multi-stage methods in `SNESMS`
33: Logically Collective
35: Level: developer
37: .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegisterDestroy()`
38: @*/
39: PetscErrorCode SNESMSRegisterAll(void)
40: {
41: PetscFunctionBegin;
42: if (SNESMSRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
43: SNESMSRegisterAllCalled = PETSC_TRUE;
45: {
46: PetscReal alpha[1] = {1.0};
47: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(SNESMSRegister(SNESMSVLTP61, 6, 1, 3.0, NULL, NULL, alpha));
89: }
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: /*@C
94: SNESMSRegisterDestroy - Frees the list of schemes that were registered by `SNESMSRegister()`.
96: Logically Collective
98: Level: developer
100: .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegister()`, `SNESMSRegisterAll()`
101: @*/
102: PetscErrorCode SNESMSRegisterDestroy(void)
103: {
104: SNESMSTableauLink link;
106: PetscFunctionBegin;
107: while ((link = SNESMSTableauList)) {
108: SNESMSTableau t = &link->tab;
109: SNESMSTableauList = link->next;
111: PetscCall(PetscFree(t->name));
112: PetscCall(PetscFree(t->gamma));
113: PetscCall(PetscFree(t->delta));
114: PetscCall(PetscFree(t->betasub));
115: PetscCall(PetscFree(link));
116: }
117: SNESMSRegisterAllCalled = PETSC_FALSE;
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /*@C
122: SNESMSInitializePackage - This function initializes everything in the `SNESMS` package. It is called
123: from `SNESInitializePackage()`.
125: Level: developer
127: .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegister()`, `SNESMSRegisterAll()`, `PetscInitialize()`
128: @*/
129: PetscErrorCode SNESMSInitializePackage(void)
130: {
131: PetscFunctionBegin;
132: if (SNESMSPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
133: SNESMSPackageInitialized = PETSC_TRUE;
135: PetscCall(SNESMSRegisterAll());
136: PetscCall(PetscRegisterFinalize(SNESMSFinalizePackage));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*@C
141: SNESMSFinalizePackage - This function destroys everything in the `SNESMS` package. It is
142: called from `PetscFinalize()`.
144: Level: developer
146: .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegister()`, `SNESMSRegisterAll()`, `SNESMSInitializePackage()`, `PetscFinalize()`
147: @*/
148: PetscErrorCode SNESMSFinalizePackage(void)
149: {
150: PetscFunctionBegin;
151: SNESMSPackageInitialized = PETSC_FALSE;
153: PetscCall(SNESMSRegisterDestroy());
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: /*@C
158: SNESMSRegister - register a multistage scheme for `SNESMS`
160: Logically Collective, No Fortran Support
162: Input Parameters:
163: + name - identifier for method
164: . nstages - number of stages
165: . nregisters - number of registers used by low-storage implementation
166: . stability - scaled stability region
167: . gamma - coefficients, see Ketcheson's paper {cite}`ketcheson2010runge`
168: . delta - coefficients, see Ketcheson's paper {cite}`ketcheson2010runge`
169: - betasub - subdiagonal of Shu-Osher form
171: Level: advanced
173: Notes:
174: The notation is described in {cite}`ketcheson2010runge` Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
176: Many multistage schemes are of the form
178: $$
179: \begin{align*}
180: X_0 = X^{(n)} \\
181: X_k = X_0 + \alpha_k * F(X_{k-1}), k = 1,\ldots,s \\
182: X^{(n+1)} = X_s
183: \end{align*}
184: $$
186: These methods can be registered with
187: .vb
188: SNESMSRegister("name",s,1,stability,NULL,NULL,alpha);
189: .ve
191: .seealso: [](ch_snes), `SNES`, `SNESMS`
192: @*/
193: PetscErrorCode SNESMSRegister(SNESMSType name, PetscInt nstages, PetscInt nregisters, PetscReal stability, const PetscReal gamma[], const PetscReal delta[], const PetscReal betasub[])
194: {
195: SNESMSTableauLink link;
196: SNESMSTableau t;
198: PetscFunctionBegin;
199: PetscAssertPointer(name, 1);
200: PetscCheck(nstages >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have at least one stage");
201: if (gamma || delta) {
202: PetscCheck(nregisters == 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support for methods written in 3-register form");
203: PetscAssertPointer(gamma, 5);
204: PetscAssertPointer(delta, 6);
205: } else {
206: PetscCheck(nregisters == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support for methods written in 1-register form");
207: }
208: PetscAssertPointer(betasub, 7);
210: PetscCall(SNESMSInitializePackage());
211: PetscCall(PetscNew(&link));
212: t = &link->tab;
213: PetscCall(PetscStrallocpy(name, &t->name));
214: t->nstages = nstages;
215: t->nregisters = nregisters;
216: t->stability = stability;
218: if (gamma && delta) {
219: PetscCall(PetscMalloc1(nstages * nregisters, &t->gamma));
220: PetscCall(PetscMalloc1(nstages, &t->delta));
221: PetscCall(PetscArraycpy(t->gamma, gamma, nstages * nregisters));
222: PetscCall(PetscArraycpy(t->delta, delta, nstages));
223: }
224: PetscCall(PetscMalloc1(nstages, &t->betasub));
225: PetscCall(PetscArraycpy(t->betasub, betasub, nstages));
227: link->next = SNESMSTableauList;
228: SNESMSTableauList = link;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*
233: X - initial state, updated in-place.
234: F - residual, computed at the initial X on input
235: */
236: static PetscErrorCode SNESMSStep_3Sstar(SNES snes, Vec X, Vec F)
237: {
238: SNES_MS *ms = (SNES_MS *)snes->data;
239: SNESMSTableau t = ms->tableau;
240: const PetscInt nstages = t->nstages;
241: const PetscReal *gamma = t->gamma, *delta = t->delta, *betasub = t->betasub;
242: Vec S1 = X, S2 = snes->work[1], S3 = snes->work[2], Y = snes->work[0], S1copy = snes->work[3];
244: PetscFunctionBegin;
245: PetscCall(VecZeroEntries(S2));
246: PetscCall(VecCopy(X, S3));
247: for (PetscInt i = 0; i < nstages; i++) {
248: Vec Ss[] = {S1copy, S2, S3, Y};
249: const PetscScalar scoeff[] = {gamma[0 * nstages + i] - 1, gamma[1 * nstages + i], gamma[2 * nstages + i], -betasub[i] * ms->damping};
251: PetscCall(VecAXPY(S2, delta[i], S1));
252: if (i > 0) PetscCall(SNESComputeFunction(snes, S1, F));
253: PetscCall(KSPSolve(snes->ksp, F, Y));
254: PetscCall(VecCopy(S1, S1copy));
255: PetscCall(VecMAXPY(S1, 4, scoeff, Ss));
256: }
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*
261: X - initial state, updated in-place.
262: F - residual, computed at the initial X on input
263: */
264: static PetscErrorCode SNESMSStep_Basic(SNES snes, Vec X, Vec F)
265: {
266: SNES_MS *ms = (SNES_MS *)snes->data;
267: SNESMSTableau tab = ms->tableau;
268: const PetscReal *alpha = tab->betasub, h = ms->damping;
269: PetscInt i, nstages = tab->nstages;
270: Vec X0 = snes->work[0];
272: PetscFunctionBegin;
273: PetscCall(VecCopy(X, X0));
274: for (i = 0; i < nstages; i++) {
275: if (i > 0) PetscCall(SNESComputeFunction(snes, X, F));
276: PetscCall(KSPSolve(snes->ksp, F, X));
277: PetscCall(VecAYPX(X, -alpha[i] * h, X0));
278: }
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: static PetscErrorCode SNESMSStep_Step(SNES snes, Vec X, Vec F)
283: {
284: SNES_MS *ms = (SNES_MS *)snes->data;
285: SNESMSTableau tab = ms->tableau;
287: PetscFunctionBegin;
288: if (tab->gamma && tab->delta) {
289: PetscCall(SNESMSStep_3Sstar(snes, X, F));
290: } else {
291: PetscCall(SNESMSStep_Basic(snes, X, F));
292: }
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode SNESMSStep_Norms(SNES snes, PetscInt iter, Vec F)
297: {
298: PetscReal fnorm;
300: PetscFunctionBegin;
301: if (SNESNeedNorm_Private(snes, iter)) {
302: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
303: SNESCheckFunctionNorm(snes, fnorm);
304: /* Monitor convergence */
305: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
306: snes->iter = iter;
307: snes->norm = fnorm;
308: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
309: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
310: /* Test for convergence */
311: PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
312: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
313: } else if (iter > 0) {
314: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
315: snes->iter = iter;
316: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
317: }
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: static PetscErrorCode SNESSolve_MS(SNES snes)
322: {
323: Vec X = snes->vec_sol, F = snes->vec_func;
324: PetscInt i;
326: PetscFunctionBegin;
327: PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
328: PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
330: snes->reason = SNES_CONVERGED_ITERATING;
331: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
332: snes->iter = 0;
333: snes->norm = 0;
334: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
336: if (!snes->vec_func_init_set) {
337: PetscCall(SNESComputeFunction(snes, X, F));
338: } else snes->vec_func_init_set = PETSC_FALSE;
340: PetscCall(SNESMSStep_Norms(snes, 0, F));
341: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
343: for (i = 0; i < snes->max_its; i++) {
344: /* Call general purpose update function */
345: PetscTryTypeMethod(snes, update, snes->iter);
347: if (i == 0 && snes->jacobian) {
348: /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
349: PetscCall(SNESComputeJacobian(snes, snes->vec_sol, snes->jacobian, snes->jacobian_pre));
350: SNESCheckJacobianDomainerror(snes);
351: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
352: }
354: PetscCall(SNESMSStep_Step(snes, X, F));
356: if (i < snes->max_its - 1 || SNESNeedNorm_Private(snes, i + 1)) PetscCall(SNESComputeFunction(snes, X, F));
358: PetscCall(SNESMSStep_Norms(snes, i + 1, F));
359: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
360: }
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode SNESSetUp_MS(SNES snes)
365: {
366: SNES_MS *ms = (SNES_MS *)snes->data;
367: SNESMSTableau tab = ms->tableau;
368: PetscInt nwork = tab->nregisters + 1; // +1 because VecMAXPY() in SNESMSStep_3Sstar()
369: // needs an extra work vec
371: PetscFunctionBegin;
372: PetscCall(SNESSetWorkVecs(snes, nwork));
373: PetscCall(SNESSetUpMatrices(snes));
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: static PetscErrorCode SNESReset_MS(SNES snes)
378: {
379: PetscFunctionBegin;
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: static PetscErrorCode SNESDestroy_MS(SNES snes)
384: {
385: PetscFunctionBegin;
386: PetscCall(SNESReset_MS(snes));
387: PetscCall(PetscFree(snes->data));
388: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", NULL));
389: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", NULL));
390: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", NULL));
391: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", NULL));
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: static PetscErrorCode SNESView_MS(SNES snes, PetscViewer viewer)
396: {
397: PetscBool iascii;
398: SNES_MS *ms = (SNES_MS *)snes->data;
399: SNESMSTableau tab = ms->tableau;
401: PetscFunctionBegin;
402: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
403: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " multi-stage method type: %s\n", tab->name));
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: static PetscErrorCode SNESSetFromOptions_MS(SNES snes, PetscOptionItems *PetscOptionsObject)
408: {
409: PetscFunctionBegin;
410: PetscOptionsHeadBegin(PetscOptionsObject, "SNES MS options");
411: {
412: SNESMSTableauLink link;
413: PetscInt count, choice;
414: PetscBool flg;
415: const char **namelist;
416: SNESMSType mstype;
417: PetscReal damping;
418: PetscBool norms = PETSC_FALSE;
420: PetscCall(SNESMSGetType(snes, &mstype));
421: for (link = SNESMSTableauList, count = 0; link; link = link->next, count++);
422: PetscCall(PetscMalloc1(count, (char ***)&namelist));
423: for (link = SNESMSTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
424: PetscCall(PetscOptionsEList("-snes_ms_type", "Multistage smoother type", "SNESMSSetType", (const char *const *)namelist, count, mstype, &choice, &flg));
425: if (flg) PetscCall(SNESMSSetType(snes, namelist[choice]));
426: PetscCall(PetscFree(namelist));
427: PetscCall(SNESMSGetDamping(snes, &damping));
428: PetscCall(PetscOptionsReal("-snes_ms_damping", "Damping for multistage method", "SNESMSSetDamping", damping, &damping, &flg));
429: if (flg) PetscCall(SNESMSSetDamping(snes, damping));
431: /* deprecated option */
432: PetscCall(PetscOptionsDeprecated("-snes_ms_norms", NULL, "3.20", "Use -snes_norm_schedule always"));
433: PetscCall(PetscOptionsBool("-snes_ms_norms", NULL, NULL, norms, &norms, NULL));
434: if (norms) PetscCall(SNESSetNormSchedule(snes, SNES_NORM_ALWAYS));
435: }
436: PetscOptionsHeadEnd();
437: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
446: *mstype = tab->name;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: static PetscErrorCode SNESMSSetType_MS(SNES snes, SNESMSType mstype)
451: {
452: SNES_MS *ms = (SNES_MS *)snes->data;
453: SNESMSTableauLink link;
454: PetscBool match;
456: PetscFunctionBegin;
457: if (ms->tableau) {
458: PetscCall(PetscStrcmp(ms->tableau->name, mstype, &match));
459: if (match) PetscFunctionReturn(PETSC_SUCCESS);
460: }
461: for (link = SNESMSTableauList; link; link = link->next) {
462: PetscCall(PetscStrcmp(link->tab.name, mstype, &match));
463: if (match) {
464: if (snes->setupcalled) PetscCall(SNESReset_MS(snes));
465: ms->tableau = &link->tab;
466: if (snes->setupcalled) PetscCall(SNESSetUp_MS(snes));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
469: }
470: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mstype);
471: }
473: /*@
474: SNESMSGetType - Get the type of multistage smoother `SNESMS`
476: Not Collective
478: Input Parameter:
479: . snes - nonlinear solver context
481: Output Parameter:
482: . mstype - type of multistage method
484: Level: advanced
486: .seealso: [](ch_snes), `SNESMS`, `SNESMSSetType()`, `SNESMSType`
487: @*/
488: PetscErrorCode SNESMSGetType(SNES snes, SNESMSType *mstype)
489: {
490: PetscFunctionBegin;
492: PetscAssertPointer(mstype, 2);
493: PetscUseMethod(snes, "SNESMSGetType_C", (SNES, SNESMSType *), (snes, mstype));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: /*@
498: SNESMSSetType - Set the type of multistage smoother `SNESMS`
500: Logically Collective
502: Input Parameters:
503: + snes - nonlinear solver context
504: - mstype - type of multistage method
506: Level: advanced
508: .seealso: [](ch_snes), `SNESMS`, `SNESMSGetType()`, `SNESMSType`
509: @*/
510: PetscErrorCode SNESMSSetType(SNES snes, SNESMSType mstype)
511: {
512: PetscFunctionBegin;
514: PetscAssertPointer(mstype, 2);
515: PetscTryMethod(snes, "SNESMSSetType_C", (SNES, SNESMSType), (snes, mstype));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode SNESMSGetDamping_MS(SNES snes, PetscReal *damping)
520: {
521: SNES_MS *ms = (SNES_MS *)snes->data;
523: PetscFunctionBegin;
524: *damping = ms->damping;
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: static PetscErrorCode SNESMSSetDamping_MS(SNES snes, PetscReal damping)
529: {
530: SNES_MS *ms = (SNES_MS *)snes->data;
532: PetscFunctionBegin;
533: ms->damping = damping;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@
538: SNESMSGetDamping - Get the damping parameter of `SNESMS` multistage scheme
540: Not Collective
542: Input Parameter:
543: . snes - nonlinear solver context
545: Output Parameter:
546: . damping - damping parameter
548: Level: advanced
550: .seealso: [](ch_snes), `SNESMSSetDamping()`, `SNESMS`
551: @*/
552: PetscErrorCode SNESMSGetDamping(SNES snes, PetscReal *damping)
553: {
554: PetscFunctionBegin;
556: PetscAssertPointer(damping, 2);
557: PetscUseMethod(snes, "SNESMSGetDamping_C", (SNES, PetscReal *), (snes, damping));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: /*@
562: SNESMSSetDamping - Set the damping parameter for a `SNESMS` multistage scheme
564: Logically Collective
566: Input Parameters:
567: + snes - nonlinear solver context
568: - damping - damping parameter
570: Level: advanced
572: .seealso: [](ch_snes), `SNESMSGetDamping()`, `SNESMS`
573: @*/
574: PetscErrorCode SNESMSSetDamping(SNES snes, PetscReal damping)
575: {
576: PetscFunctionBegin;
579: PetscTryMethod(snes, "SNESMSSetDamping_C", (SNES, PetscReal), (snes, damping));
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: /*MC
584: SNESMS - multi-stage smoothers
586: Options Database Keys:
587: + -snes_ms_type - type of multi-stage smoother
588: - -snes_ms_damping - damping for multi-stage method
590: Level: advanced
592: Notes:
593: These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
594: FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
596: Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
598: The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with `SNESMSRegister()`.
600: See {cite}`ketcheson2010runge`, {cite}`jameson1983`, {cite}`pierce1997preconditioned`, and {cite}`va1981design`
602: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESMS`, `SNESFAS`, `KSPCHEBYSHEV`, `SNESMSSetDamping()`, `SNESMSGetDamping()`,
603: `SNESMSSetType()`, `SNESMSGetType()`
604: M*/
606: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
607: {
608: SNES_MS *ms;
610: PetscFunctionBegin;
611: PetscCall(SNESMSInitializePackage());
613: snes->ops->setup = SNESSetUp_MS;
614: snes->ops->solve = SNESSolve_MS;
615: snes->ops->destroy = SNESDestroy_MS;
616: snes->ops->setfromoptions = SNESSetFromOptions_MS;
617: snes->ops->view = SNESView_MS;
618: snes->ops->reset = SNESReset_MS;
620: snes->usesnpc = PETSC_FALSE;
621: snes->usesksp = PETSC_TRUE;
623: snes->alwayscomputesfinalresidual = PETSC_FALSE;
625: PetscCall(SNESParametersInitialize(snes));
627: PetscCall(PetscNew(&ms));
628: snes->data = (void *)ms;
629: ms->damping = 0.9;
631: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", SNESMSGetType_MS));
632: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", SNESMSSetType_MS));
633: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", SNESMSGetDamping_MS));
634: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", SNESMSSetDamping_MS));
636: PetscCall(SNESMSSetType(snes, SNESMSDefault));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }