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:     SNESCheckFunctionDomainError(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) PetscCall(SNESComputeFunction(snes, X, F));
    337:   else snes->vec_func_init_set = PETSC_FALSE;

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

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

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

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

    355:     if (i < snes->max_its - 1 || SNESNeedNorm_Private(snes, i + 1)) PetscCall(SNESComputeFunction(snes, X, F));

    357:     PetscCall(SNESMSStep_Norms(snes, i + 1, F));
    358:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
    359:   }
    360:   PetscFunctionReturn(PETSC_SUCCESS);
    361: }

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

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

    376: static PetscErrorCode SNESDestroy_MS(SNES snes)
    377: {
    378:   PetscFunctionBegin;
    379:   PetscCall(PetscFree(snes->data));
    380:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", NULL));
    381:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", NULL));
    382:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", NULL));
    383:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", NULL));
    384:   PetscFunctionReturn(PETSC_SUCCESS);
    385: }

    387: static PetscErrorCode SNESView_MS(SNES snes, PetscViewer viewer)
    388: {
    389:   PetscBool     isascii;
    390:   SNES_MS      *ms  = (SNES_MS *)snes->data;
    391:   SNESMSTableau tab = ms->tableau;

    393:   PetscFunctionBegin;
    394:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    395:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  multi-stage method type: %s\n", tab->name));
    396:   PetscFunctionReturn(PETSC_SUCCESS);
    397: }

    399: static PetscErrorCode SNESSetFromOptions_MS(SNES snes, PetscOptionItems PetscOptionsObject)
    400: {
    401:   PetscFunctionBegin;
    402:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES MS options");
    403:   {
    404:     SNESMSTableauLink link;
    405:     PetscInt          count, choice;
    406:     PetscBool         flg;
    407:     const char      **namelist;
    408:     SNESMSType        mstype;
    409:     PetscReal         damping;
    410:     PetscBool         norms = PETSC_FALSE;

    412:     PetscCall(SNESMSGetType(snes, &mstype));
    413:     for (link = SNESMSTableauList, count = 0; link; link = link->next, count++);
    414:     PetscCall(PetscMalloc1(count, (char ***)&namelist));
    415:     for (link = SNESMSTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
    416:     PetscCall(PetscOptionsEList("-snes_ms_type", "Multistage smoother type", "SNESMSSetType", (const char *const *)namelist, count, mstype, &choice, &flg));
    417:     if (flg) PetscCall(SNESMSSetType(snes, namelist[choice]));
    418:     PetscCall(PetscFree(namelist));
    419:     PetscCall(SNESMSGetDamping(snes, &damping));
    420:     PetscCall(PetscOptionsReal("-snes_ms_damping", "Damping for multistage method", "SNESMSSetDamping", damping, &damping, &flg));
    421:     if (flg) PetscCall(SNESMSSetDamping(snes, damping));

    423:     /* deprecated option */
    424:     PetscCall(PetscOptionsDeprecated("-snes_ms_norms", NULL, "3.20", "Use -snes_norm_schedule always"));
    425:     PetscCall(PetscOptionsBool("-snes_ms_norms", NULL, NULL, norms, &norms, NULL));
    426:     if (norms) PetscCall(SNESSetNormSchedule(snes, SNES_NORM_ALWAYS));
    427:   }
    428:   PetscOptionsHeadEnd();
    429:   PetscFunctionReturn(PETSC_SUCCESS);
    430: }

    432: static PetscErrorCode SNESMSGetType_MS(SNES snes, SNESMSType *mstype)
    433: {
    434:   SNES_MS      *ms  = (SNES_MS *)snes->data;
    435:   SNESMSTableau tab = ms->tableau;

    437:   PetscFunctionBegin;
    438:   *mstype = tab->name;
    439:   PetscFunctionReturn(PETSC_SUCCESS);
    440: }

    442: static PetscErrorCode SNESMSSetType_MS(SNES snes, SNESMSType mstype)
    443: {
    444:   SNES_MS          *ms = (SNES_MS *)snes->data;
    445:   SNESMSTableauLink link;
    446:   PetscBool         match;

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

    464: /*@
    465:   SNESMSGetType - Get the type of multistage smoother `SNESMS`

    467:   Not Collective

    469:   Input Parameter:
    470: . snes - nonlinear solver context

    472:   Output Parameter:
    473: . mstype - type of multistage method

    475:   Level: advanced

    477: .seealso: [](ch_snes), `SNESMS`, `SNESMSSetType()`, `SNESMSType`
    478: @*/
    479: PetscErrorCode SNESMSGetType(SNES snes, SNESMSType *mstype)
    480: {
    481:   PetscFunctionBegin;
    483:   PetscAssertPointer(mstype, 2);
    484:   PetscUseMethod(snes, "SNESMSGetType_C", (SNES, SNESMSType *), (snes, mstype));
    485:   PetscFunctionReturn(PETSC_SUCCESS);
    486: }

    488: /*@
    489:   SNESMSSetType - Set the type of multistage smoother `SNESMS`

    491:   Logically Collective

    493:   Input Parameters:
    494: + snes   - nonlinear solver context
    495: - mstype - type of multistage method

    497:   Level: advanced

    499: .seealso: [](ch_snes), `SNESMS`, `SNESMSGetType()`, `SNESMSType`
    500: @*/
    501: PetscErrorCode SNESMSSetType(SNES snes, SNESMSType mstype)
    502: {
    503:   PetscFunctionBegin;
    505:   PetscAssertPointer(mstype, 2);
    506:   PetscTryMethod(snes, "SNESMSSetType_C", (SNES, SNESMSType), (snes, mstype));
    507:   PetscFunctionReturn(PETSC_SUCCESS);
    508: }

    510: static PetscErrorCode SNESMSGetDamping_MS(SNES snes, PetscReal *damping)
    511: {
    512:   SNES_MS *ms = (SNES_MS *)snes->data;

    514:   PetscFunctionBegin;
    515:   *damping = ms->damping;
    516:   PetscFunctionReturn(PETSC_SUCCESS);
    517: }

    519: static PetscErrorCode SNESMSSetDamping_MS(SNES snes, PetscReal damping)
    520: {
    521:   SNES_MS *ms = (SNES_MS *)snes->data;

    523:   PetscFunctionBegin;
    524:   ms->damping = damping;
    525:   PetscFunctionReturn(PETSC_SUCCESS);
    526: }

    528: /*@
    529:   SNESMSGetDamping - Get the damping parameter of `SNESMS` multistage scheme

    531:   Not Collective

    533:   Input Parameter:
    534: . snes - nonlinear solver context

    536:   Output Parameter:
    537: . damping - damping parameter

    539:   Level: advanced

    541: .seealso: [](ch_snes), `SNESMSSetDamping()`, `SNESMS`
    542: @*/
    543: PetscErrorCode SNESMSGetDamping(SNES snes, PetscReal *damping)
    544: {
    545:   PetscFunctionBegin;
    547:   PetscAssertPointer(damping, 2);
    548:   PetscUseMethod(snes, "SNESMSGetDamping_C", (SNES, PetscReal *), (snes, damping));
    549:   PetscFunctionReturn(PETSC_SUCCESS);
    550: }

    552: /*@
    553:   SNESMSSetDamping - Set the damping parameter for a `SNESMS` multistage scheme

    555:   Logically Collective

    557:   Input Parameters:
    558: + snes    - nonlinear solver context
    559: - damping - damping parameter

    561:   Level: advanced

    563: .seealso: [](ch_snes), `SNESMSGetDamping()`, `SNESMS`
    564: @*/
    565: PetscErrorCode SNESMSSetDamping(SNES snes, PetscReal damping)
    566: {
    567:   PetscFunctionBegin;
    570:   PetscTryMethod(snes, "SNESMSSetDamping_C", (SNES, PetscReal), (snes, damping));
    571:   PetscFunctionReturn(PETSC_SUCCESS);
    572: }

    574: /*MC
    575:       SNESMS - multi-stage smoothers

    577:       Options Database Keys:
    578: +     -snes_ms_type    - type of multi-stage smoother
    579: -     -snes_ms_damping - damping for multi-stage method

    581:       Level: advanced

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

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

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

    591:       See {cite}`ketcheson2010runge`, {cite}`jameson1983`, {cite}`pierce1997preconditioned`, and {cite}`va1981design`

    593: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESMS`, `SNESFAS`, `KSPCHEBYSHEV`, `SNESMSSetDamping()`, `SNESMSGetDamping()`,
    594:           `SNESMSSetType()`, `SNESMSGetType()`
    595: M*/

    597: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
    598: {
    599:   SNES_MS *ms;

    601:   PetscFunctionBegin;
    602:   PetscCall(SNESMSInitializePackage());

    604:   snes->ops->setup          = SNESSetUp_MS;
    605:   snes->ops->solve          = SNESSolve_MS;
    606:   snes->ops->destroy        = SNESDestroy_MS;
    607:   snes->ops->setfromoptions = SNESSetFromOptions_MS;
    608:   snes->ops->view           = SNESView_MS;

    610:   snes->usesnpc = PETSC_FALSE;
    611:   snes->usesksp = PETSC_TRUE;

    613:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

    615:   PetscCall(SNESParametersInitialize(snes));

    617:   PetscCall(PetscNew(&ms));
    618:   snes->data  = (void *)ms;
    619:   ms->damping = 0.9;

    621:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", SNESMSGetType_MS));
    622:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", SNESMSSetType_MS));
    623:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", SNESMSGetDamping_MS));
    624:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", SNESMSSetDamping_MS));

    626:   PetscCall(SNESMSSetType(snes, SNESMSDefault));
    627:   PetscFunctionReturn(PETSC_SUCCESS);
    628: }