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:   PetscFunctionBegin;
 43:   if (SNESMSRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
 44:   SNESMSRegisterAllCalled = PETSC_TRUE;

 46:   {
 47:     PetscReal alpha[1] = {1.0};
 48:     PetscCall(SNESMSRegister(SNESMSEULER, 1, 1, 1.0, NULL, NULL, alpha));
 49:   }

 51:   {
 52:     PetscReal gamma[3][6] = {
 53:       {0.0000000000000000E+00, -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02,  -1.4204296130641869E-01},
 54:       {1.0000000000000000E+00, 1.1111025767083920E+00,  5.6150921583923230E-01,  7.4151723494934041E-01,  3.1714538168600587E-01,  4.6479276238548706E-01 },
 55:       {0.0000000000000000E+00, 0.0000000000000000E+00,  0.0000000000000000E+00,  6.7968174970583317E-01,  -4.1755042846051737E-03, -1.9115668129923846E-01}
 56:     };
 57:     PetscReal delta[6]   = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
 58:     PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
 59:     PetscCall(SNESMSRegister(SNESMSM62, 6, 3, 1.0, &gamma[0][0], delta, betasub));
 60:   }

 62:   { /* Jameson (1983) */
 63:     PetscReal alpha[4] = {0.25, 0.5, 0.55, 1.0};
 64:     PetscCall(SNESMSRegister(SNESMSJAMESON83, 4, 1, 1.0, NULL, NULL, alpha));
 65:   }

 67:   { /* Van Leer, Tai, and Powell (1989) 1 stage, order 1 */
 68:     PetscReal alpha[1] = {1.0};
 69:     PetscCall(SNESMSRegister(SNESMSVLTP11, 1, 1, 0.5, NULL, NULL, alpha));
 70:   }
 71:   { /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
 72:     PetscReal alpha[2] = {0.3333, 1.0};
 73:     PetscCall(SNESMSRegister(SNESMSVLTP21, 2, 1, 1.0, NULL, NULL, alpha));
 74:   }
 75:   { /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
 76:     PetscReal alpha[3] = {0.1481, 0.4000, 1.0};
 77:     PetscCall(SNESMSRegister(SNESMSVLTP31, 3, 1, 1.5, NULL, NULL, alpha));
 78:   }
 79:   { /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
 80:     PetscReal alpha[4] = {0.0833, 0.2069, 0.4265, 1.0};
 81:     PetscCall(SNESMSRegister(SNESMSVLTP41, 4, 1, 2.0, NULL, NULL, alpha));
 82:   }
 83:   { /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
 84:     PetscReal alpha[5] = {0.0533, 0.1263, 0.2375, 0.4414, 1.0};
 85:     PetscCall(SNESMSRegister(SNESMSVLTP51, 5, 1, 2.5, NULL, NULL, alpha));
 86:   }
 87:   { /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
 88:     PetscReal alpha[6] = {0.0370, 0.0851, 0.1521, 0.2562, 0.4512, 1.0};
 89:     PetscCall(SNESMSRegister(SNESMSVLTP61, 6, 1, 3.0, NULL, NULL, alpha));
 90:   }
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

 94: /*@C
 95:    SNESMSRegisterDestroy - Frees the list of schemes that were registered by `SNESMSRegister()`.

 97:    Not Collective

 99:    Level: developer

101: .seealso: `SNESMSRegister()`, `SNESMSRegisterAll()`
102: @*/
103: PetscErrorCode SNESMSRegisterDestroy(void)
104: {
105:   SNESMSTableauLink link;

107:   PetscFunctionBegin;
108:   while ((link = SNESMSTableauList)) {
109:     SNESMSTableau t   = &link->tab;
110:     SNESMSTableauList = link->next;

112:     PetscCall(PetscFree(t->name));
113:     PetscCall(PetscFree(t->gamma));
114:     PetscCall(PetscFree(t->delta));
115:     PetscCall(PetscFree(t->betasub));
116:     PetscCall(PetscFree(link));
117:   }
118:   SNESMSRegisterAllCalled = PETSC_FALSE;
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: /*@C
123:   SNESMSInitializePackage - This function initializes everything in the `SNESMS` package. It is called
124:   from `SNESInitializePackage()`.

126:   Level: developer

128: .seealso: `PetscInitialize()`
129: @*/
130: PetscErrorCode SNESMSInitializePackage(void)
131: {
132:   PetscFunctionBegin;
133:   if (SNESMSPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
134:   SNESMSPackageInitialized = PETSC_TRUE;

136:   PetscCall(SNESMSRegisterAll());
137:   PetscCall(PetscRegisterFinalize(SNESMSFinalizePackage));
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: /*@C
142:   SNESMSFinalizePackage - This function destroys everything in the `SNESMS` package. It is
143:   called from `PetscFinalize()`.

145:   Level: developer

147: .seealso: `PetscFinalize()`
148: @*/
149: PetscErrorCode SNESMSFinalizePackage(void)
150: {
151:   PetscFunctionBegin;
152:   SNESMSPackageInitialized = PETSC_FALSE;

154:   PetscCall(SNESMSRegisterDestroy());
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: /*@C
159:    SNESMSRegister - register a multistage scheme for `SNESMS`

161:    Logically Collective

163:    Input Parameters:
164: +  name - identifier for method
165: .  nstages - number of stages
166: .  nregisters - number of registers used by low-storage implementation
167: .  stability - scaled stability region
168: .  gamma - coefficients, see Ketcheson's paper
169: .  delta - coefficients, see Ketcheson's paper
170: -  betasub - subdiagonal of Shu-Osher form

172:    Level: advanced

174:    Notes:
175:    The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.

177:    Many multistage schemes are of the form
178:    $ X_0 = X^{(n)}
179:    $ X_k = X_0 + \alpha_k * F(X_{k-1}), k = 1,\ldots,s
180:    $ X^{(n+1)} = X_s
181:    These methods can be registered with
182: .vb
183:    SNESMSRegister("name",s,1,stability,NULL,NULL,alpha);
184: .ve

186: .seealso: `SNESMS`
187: @*/
188: PetscErrorCode SNESMSRegister(SNESMSType name, PetscInt nstages, PetscInt nregisters, PetscReal stability, const PetscReal gamma[], const PetscReal delta[], const PetscReal betasub[])
189: {
190:   SNESMSTableauLink link;
191:   SNESMSTableau     t;

193:   PetscFunctionBegin;
195:   PetscCheck(nstages >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have at least one stage");
196:   if (gamma || delta) {
197:     PetscCheck(nregisters == 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support for methods written in 3-register form");
200:   } else {
201:     PetscCheck(nregisters == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support for methods written in 1-register form");
202:   }

205:   PetscCall(SNESMSInitializePackage());
206:   PetscCall(PetscNew(&link));
207:   t = &link->tab;
208:   PetscCall(PetscStrallocpy(name, &t->name));
209:   t->nstages    = nstages;
210:   t->nregisters = nregisters;
211:   t->stability  = stability;

213:   if (gamma && delta) {
214:     PetscCall(PetscMalloc1(nstages * nregisters, &t->gamma));
215:     PetscCall(PetscMalloc1(nstages, &t->delta));
216:     PetscCall(PetscArraycpy(t->gamma, gamma, nstages * nregisters));
217:     PetscCall(PetscArraycpy(t->delta, delta, nstages));
218:   }
219:   PetscCall(PetscMalloc1(nstages, &t->betasub));
220:   PetscCall(PetscArraycpy(t->betasub, betasub, nstages));

222:   link->next        = SNESMSTableauList;
223:   SNESMSTableauList = link;
224:   PetscFunctionReturn(PETSC_SUCCESS);
225: }

227: /*
228:   X - initial state, updated in-place.
229:   F - residual, computed at the initial X on input
230: */
231: static PetscErrorCode SNESMSStep_3Sstar(SNES snes, Vec X, Vec F)
232: {
233:   SNES_MS         *ms      = (SNES_MS *)snes->data;
234:   SNESMSTableau    t       = ms->tableau;
235:   const PetscInt   nstages = t->nstages;
236:   const PetscReal *gamma = t->gamma, *delta = t->delta, *betasub = t->betasub;
237:   Vec              S1 = X, S2 = snes->work[1], S3 = snes->work[2], Y = snes->work[0], S1copy = snes->work[3];

239:   PetscFunctionBegin;
240:   PetscCall(VecZeroEntries(S2));
241:   PetscCall(VecCopy(X, S3));
242:   for (PetscInt i = 0; i < nstages; i++) {
243:     Vec               Ss[]     = {S1copy, S2, S3, Y};
244:     const PetscScalar scoeff[] = {gamma[0 * nstages + i] - 1, gamma[1 * nstages + i], gamma[2 * nstages + i], -betasub[i] * ms->damping};

246:     PetscCall(VecAXPY(S2, delta[i], S1));
247:     if (i > 0) PetscCall(SNESComputeFunction(snes, S1, F));
248:     PetscCall(KSPSolve(snes->ksp, F, Y));
249:     PetscCall(VecCopy(S1, S1copy));
250:     PetscCall(VecMAXPY(S1, 4, scoeff, Ss));
251:   }
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

255: /*
256:   X - initial state, updated in-place.
257:   F - residual, computed at the initial X on input
258: */
259: static PetscErrorCode SNESMSStep_Basic(SNES snes, Vec X, Vec F)
260: {
261:   SNES_MS         *ms    = (SNES_MS *)snes->data;
262:   SNESMSTableau    tab   = ms->tableau;
263:   const PetscReal *alpha = tab->betasub, h = ms->damping;
264:   PetscInt         i, nstages              = tab->nstages;
265:   Vec              X0 = snes->work[0];

267:   PetscFunctionBegin;
268:   PetscCall(VecCopy(X, X0));
269:   for (i = 0; i < nstages; i++) {
270:     if (i > 0) PetscCall(SNESComputeFunction(snes, X, F));
271:     PetscCall(KSPSolve(snes->ksp, F, X));
272:     PetscCall(VecAYPX(X, -alpha[i] * h, X0));
273:   }
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: static PetscErrorCode SNESMSStep_Step(SNES snes, Vec X, Vec F)
278: {
279:   SNES_MS      *ms  = (SNES_MS *)snes->data;
280:   SNESMSTableau tab = ms->tableau;

282:   PetscFunctionBegin;
283:   if (tab->gamma && tab->delta) {
284:     PetscCall(SNESMSStep_3Sstar(snes, X, F));
285:   } else {
286:     PetscCall(SNESMSStep_Basic(snes, X, F));
287:   }
288:   PetscFunctionReturn(PETSC_SUCCESS);
289: }

291: static PetscErrorCode SNESMSStep_Norms(SNES snes, PetscInt iter, Vec F)
292: {
293:   SNES_MS  *ms = (SNES_MS *)snes->data;
294:   PetscReal fnorm;

296:   PetscFunctionBegin;
297:   if (ms->norms) {
298:     PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
299:     SNESCheckFunctionNorm(snes, fnorm);
300:     /* Monitor convergence */
301:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
302:     snes->iter = iter;
303:     snes->norm = fnorm;
304:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
305:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
306:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
307:     /* Test for convergence */
308:     PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
309:   } else if (iter > 0) {
310:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
311:     snes->iter = iter;
312:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
313:   }
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: static PetscErrorCode SNESSolve_MS(SNES snes)
318: {
319:   SNES_MS *ms = (SNES_MS *)snes->data;
320:   Vec      X = snes->vec_sol, F = snes->vec_func;
321:   PetscInt i;

323:   PetscFunctionBegin;
324:   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);
325:   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));

327:   snes->reason = SNES_CONVERGED_ITERATING;
328:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
329:   snes->iter = 0;
330:   snes->norm = 0;
331:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

333:   if (!snes->vec_func_init_set) {
334:     PetscCall(SNESComputeFunction(snes, X, F));
335:   } else snes->vec_func_init_set = PETSC_FALSE;

337:   PetscCall(SNESMSStep_Norms(snes, 0, F));
338:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

340:   for (i = 0; i < snes->max_its; i++) {
341:     /* Call general purpose update function */
342:     PetscTryTypeMethod(snes, update, snes->iter);

344:     if (i == 0 && snes->jacobian) {
345:       /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
346:       PetscCall(SNESComputeJacobian(snes, snes->vec_sol, snes->jacobian, snes->jacobian_pre));
347:       SNESCheckJacobianDomainerror(snes);
348:       PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
349:     }

351:     PetscCall(SNESMSStep_Step(snes, X, F));

353:     if (i < snes->max_its - 1 || ms->norms) PetscCall(SNESComputeFunction(snes, X, F));

355:     PetscCall(SNESMSStep_Norms(snes, i + 1, F));
356:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
357:   }
358:   if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
359:   PetscFunctionReturn(PETSC_SUCCESS);
360: }

362: static PetscErrorCode SNESSetUp_MS(SNES snes)
363: {
364:   SNES_MS      *ms    = (SNES_MS *)snes->data;
365:   SNESMSTableau tab   = ms->tableau;
366:   PetscInt      nwork = tab->nregisters + 1; // +1 because VecMAXPY() in SNESMSStep_3Sstar()
367:                                              // needs an extra work vec

369:   PetscFunctionBegin;
370:   PetscCall(SNESSetWorkVecs(snes, nwork));
371:   PetscCall(SNESSetUpMatrices(snes));
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: static PetscErrorCode SNESReset_MS(SNES snes)
376: {
377:   PetscFunctionBegin;
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: static PetscErrorCode SNESDestroy_MS(SNES snes)
382: {
383:   PetscFunctionBegin;
384:   PetscCall(SNESReset_MS(snes));
385:   PetscCall(PetscFree(snes->data));
386:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", NULL));
387:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", NULL));
388:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", NULL));
389:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", NULL));
390:   PetscFunctionReturn(PETSC_SUCCESS);
391: }

393: static PetscErrorCode SNESView_MS(SNES snes, PetscViewer viewer)
394: {
395:   PetscBool     iascii;
396:   SNES_MS      *ms  = (SNES_MS *)snes->data;
397:   SNESMSTableau tab = ms->tableau;

399:   PetscFunctionBegin;
400:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
401:   if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  multi-stage method type: %s\n", tab->name));
402:   PetscFunctionReturn(PETSC_SUCCESS);
403: }

405: static PetscErrorCode SNESSetFromOptions_MS(SNES snes, PetscOptionItems *PetscOptionsObject)
406: {
407:   SNES_MS *ms = (SNES_MS *)snes->data;

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;

419:     PetscCall(SNESMSGetType(snes, &mstype));
420:     for (link = SNESMSTableauList, count = 0; link; link = link->next, count++)
421:       ;
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));
430:     PetscCall(PetscOptionsBool("-snes_ms_norms", "Compute norms for monitoring", "none", ms->norms, &ms->norms, NULL));
431:   }
432:   PetscOptionsHeadEnd();
433:   PetscFunctionReturn(PETSC_SUCCESS);
434: }

436: static PetscErrorCode SNESMSGetType_MS(SNES snes, SNESMSType *mstype)
437: {
438:   SNES_MS      *ms  = (SNES_MS *)snes->data;
439:   SNESMSTableau tab = ms->tableau;

441:   PetscFunctionBegin;
442:   *mstype = tab->name;
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: static PetscErrorCode SNESMSSetType_MS(SNES snes, SNESMSType mstype)
447: {
448:   SNES_MS          *ms = (SNES_MS *)snes->data;
449:   SNESMSTableauLink link;
450:   PetscBool         match;

452:   PetscFunctionBegin;
453:   if (ms->tableau) {
454:     PetscCall(PetscStrcmp(ms->tableau->name, mstype, &match));
455:     if (match) PetscFunctionReturn(PETSC_SUCCESS);
456:   }
457:   for (link = SNESMSTableauList; link; link = link->next) {
458:     PetscCall(PetscStrcmp(link->tab.name, mstype, &match));
459:     if (match) {
460:       if (snes->setupcalled) PetscCall(SNESReset_MS(snes));
461:       ms->tableau = &link->tab;
462:       if (snes->setupcalled) PetscCall(SNESSetUp_MS(snes));
463:       PetscFunctionReturn(PETSC_SUCCESS);
464:     }
465:   }
466:   SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mstype);
467: }

469: /*@C
470:   SNESMSGetType - Get the type of multistage smoother `SNESMS`

472:   Not Collective

474:   Input Parameter:
475: .  snes - nonlinear solver context

477:   Output Parameter:
478: .  mstype - type of multistage method

480:   Level: advanced

482: .seealso: `SNESMS`, `SNESMSSetType()`, `SNESMSType`, `SNESMS`
483: @*/
484: PetscErrorCode SNESMSGetType(SNES snes, SNESMSType *mstype)
485: {
486:   PetscFunctionBegin;
489:   PetscUseMethod(snes, "SNESMSGetType_C", (SNES, SNESMSType *), (snes, mstype));
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

493: /*@C
494:   SNESMSSetType - Set the type of multistage smoother `SNESMS`

496:   Logically Collective

498:   Input Parameters:
499: +  snes - nonlinear solver context
500: -  mstype - type of multistage method

502:   Level: advanced

504: .seealso: `SNESMS`, `SNESMSGetType()`, `SNESMSType`, `SNESMS`
505: @*/
506: PetscErrorCode SNESMSSetType(SNES snes, SNESMSType mstype)
507: {
508:   PetscFunctionBegin;
511:   PetscTryMethod(snes, "SNESMSSetType_C", (SNES, SNESMSType), (snes, mstype));
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

515: static PetscErrorCode SNESMSGetDamping_MS(SNES snes, PetscReal *damping)
516: {
517:   SNES_MS *ms = (SNES_MS *)snes->data;

519:   PetscFunctionBegin;
520:   *damping = ms->damping;
521:   PetscFunctionReturn(PETSC_SUCCESS);
522: }

524: static PetscErrorCode SNESMSSetDamping_MS(SNES snes, PetscReal damping)
525: {
526:   SNES_MS *ms = (SNES_MS *)snes->data;

528:   PetscFunctionBegin;
529:   ms->damping = damping;
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

533: /*@C
534:   SNESMSGetDamping - Get the damping parameter of `SNESMS` multistage scheme

536:   Not Collective

538:   Input Parameter:
539: .  snes - nonlinear solver context

541:   Output Parameter:
542: .  damping - damping parameter

544:   Level: advanced

546: .seealso: `SNESMSSetDamping()`, `SNESMS`
547: @*/
548: PetscErrorCode SNESMSGetDamping(SNES snes, PetscReal *damping)
549: {
550:   PetscFunctionBegin;
553:   PetscUseMethod(snes, "SNESMSGetDamping_C", (SNES, PetscReal *), (snes, damping));
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: /*@C
558:   SNESMSSetDamping - Set the damping parameter for a `SNESMS` multistage scheme

560:   Logically Collective

562:   Input Parameters:
563: +  snes - nonlinear solver context
564: -  damping - damping parameter

566:   Level: advanced

568: .seealso: `SNESMSGetDamping()`, `SNESMS`
569: @*/
570: PetscErrorCode SNESMSSetDamping(SNES snes, PetscReal damping)
571: {
572:   PetscFunctionBegin;
575:   PetscTryMethod(snes, "SNESMSSetDamping_C", (SNES, PetscReal), (snes, damping));
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: /*MC
580:       SNESMS - multi-stage smoothers

582:       Options Database Keys:
583: +     -snes_ms_type - type of multi-stage smoother
584: -     -snes_ms_damping - damping for multi-stage method

586:       Notes:
587:       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
588:       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).

590:       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.

592:       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with `SNESMSRegister()`.

594:       References:
595: +     * - Ketcheson (2010) Runge Kutta methods with minimum storage implementations (https://doi.org/10.1016/j.jcp.2009.11.006).
596: .     * - 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).
597: .     * - Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes (https://doi.org/10.1006/jcph.1997.5772).
598: -     * - 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).

600:       Level: advanced

602: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESMS`, `SNESFAS`, `KSPCHEBYSHEV`

604: M*/
605: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
606: {
607:   SNES_MS *ms;

609:   PetscFunctionBegin;
610:   PetscCall(SNESMSInitializePackage());

612:   snes->ops->setup          = SNESSetUp_MS;
613:   snes->ops->solve          = SNESSolve_MS;
614:   snes->ops->destroy        = SNESDestroy_MS;
615:   snes->ops->setfromoptions = SNESSetFromOptions_MS;
616:   snes->ops->view           = SNESView_MS;
617:   snes->ops->reset          = SNESReset_MS;

619:   snes->usesnpc = PETSC_FALSE;
620:   snes->usesksp = PETSC_TRUE;

622:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

624:   PetscCall(PetscNew(&ms));
625:   snes->data  = (void *)ms;
626:   ms->damping = 0.9;
627:   ms->norms   = PETSC_FALSE;

629:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", SNESMSGetType_MS));
630:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", SNESMSSetType_MS));
631:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", SNESMSGetDamping_MS));
632:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", SNESMSSetDamping_MS));

634:   PetscCall(SNESMSSetType(snes, SNESMSDefault));
635:   PetscFunctionReturn(PETSC_SUCCESS);
636: }