Actual source code: glle.c

      1: #include <../src/ts/impls/implicit/glle/glle.h>
      2: #include <petscdm.h>
      3: #include <petscblaslapack.h>

      5: static const char       *TSGLLEErrorDirections[] = {"FORWARD", "BACKWARD", "TSGLLEErrorDirection", "TSGLLEERROR_", NULL};
      6: static PetscFunctionList TSGLLEList;
      7: static PetscFunctionList TSGLLEAcceptList;
      8: static PetscBool         TSGLLEPackageInitialized;
      9: static PetscBool         TSGLLERegisterAllCalled;

     11: /* This function is pure */
     12: static PetscScalar Factorial(PetscInt n)
     13: {
     14:   PetscInt i;
     15:   if (n < 12) { /* Can compute with 32-bit integers */
     16:     PetscInt f = 1;
     17:     for (i = 2; i <= n; i++) f *= i;
     18:     return (PetscScalar)f;
     19:   } else {
     20:     PetscScalar f = 1.;
     21:     for (i = 2; i <= n; i++) f *= (PetscScalar)i;
     22:     return f;
     23:   }
     24: }

     26: /* This function is pure */
     27: static PetscScalar CPowF(PetscScalar c, PetscInt p)
     28: {
     29:   return PetscPowRealInt(PetscRealPart(c), p) / Factorial(p);
     30: }

     32: static PetscErrorCode TSGLLEGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
     33: {
     34:   TS_GLLE *gl = (TS_GLLE *)ts->data;

     36:   PetscFunctionBegin;
     37:   if (Z) {
     38:     if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Z", Z));
     39:     else *Z = gl->Z;
     40:   }
     41:   if (Ydotstage) {
     42:     if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage));
     43:     else *Ydotstage = gl->Ydot[gl->stage];
     44:   }
     45:   PetscFunctionReturn(PETSC_SUCCESS);
     46: }

     48: static PetscErrorCode TSGLLERestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
     49: {
     50:   PetscFunctionBegin;
     51:   if (Z) {
     52:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Z", Z));
     53:   }
     54:   if (Ydotstage) {
     55:     if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage));
     56:   }
     57:   PetscFunctionReturn(PETSC_SUCCESS);
     58: }

     60: static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, PetscCtx ctx)
     61: {
     62:   PetscFunctionBegin;
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
     67: {
     68:   TS  ts = (TS)ctx;
     69:   Vec Ydot, Ydot_c;

     71:   PetscFunctionBegin;
     72:   PetscCall(TSGLLEGetVecs(ts, fine, NULL, &Ydot));
     73:   PetscCall(TSGLLEGetVecs(ts, coarse, NULL, &Ydot_c));
     74:   PetscCall(MatRestrict(restrct, Ydot, Ydot_c));
     75:   PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c));
     76:   PetscCall(TSGLLERestoreVecs(ts, fine, NULL, &Ydot));
     77:   PetscCall(TSGLLERestoreVecs(ts, coarse, NULL, &Ydot_c));
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, PetscCtx ctx)
     82: {
     83:   PetscFunctionBegin;
     84:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

     87: static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
     88: {
     89:   TS  ts = (TS)ctx;
     90:   Vec Ydot, Ydot_s;

     92:   PetscFunctionBegin;
     93:   PetscCall(TSGLLEGetVecs(ts, dm, NULL, &Ydot));
     94:   PetscCall(TSGLLEGetVecs(ts, subdm, NULL, &Ydot_s));

     96:   PetscCall(VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));
     97:   PetscCall(VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));

     99:   PetscCall(TSGLLERestoreVecs(ts, dm, NULL, &Ydot));
    100:   PetscCall(TSGLLERestoreVecs(ts, subdm, NULL, &Ydot_s));
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: static PetscErrorCode TSGLLESchemeCreate(PetscInt p, PetscInt q, PetscInt r, PetscInt s, const PetscScalar *c, const PetscScalar *a, const PetscScalar *b, const PetscScalar *u, const PetscScalar *v, TSGLLEScheme *inscheme)
    105: {
    106:   TSGLLEScheme scheme;
    107:   PetscInt     j;

    109:   PetscFunctionBegin;
    110:   PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scheme order must be positive");
    111:   PetscCheck(r >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one item must be carried between steps");
    112:   PetscCheck(s >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one stage is required");
    113:   PetscAssertPointer(inscheme, 10);
    114:   *inscheme = NULL;
    115:   PetscCall(PetscNew(&scheme));
    116:   scheme->p = p;
    117:   scheme->q = q;
    118:   scheme->r = r;
    119:   scheme->s = s;

    121:   PetscCall(PetscMalloc5(s, &scheme->c, s * s, &scheme->a, r * s, &scheme->b, r * s, &scheme->u, r * r, &scheme->v));
    122:   PetscCall(PetscArraycpy(scheme->c, c, s));
    123:   for (j = 0; j < s * s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j];
    124:   for (j = 0; j < r * s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j];
    125:   for (j = 0; j < s * r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j];
    126:   for (j = 0; j < r * r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j];

    128:   PetscCall(PetscMalloc6(r, &scheme->alpha, r, &scheme->beta, r, &scheme->gamma, 3 * s, &scheme->phi, 3 * r, &scheme->psi, r, &scheme->stage_error));
    129:   {
    130:     PetscInt     i, j, k, ss = s + 2;
    131:     PetscBLASInt m, n, one = 1, *ipiv, lwork, info, ldb;
    132:     PetscReal    rcond, *sing, *workreal;
    133:     PetscScalar *ImV, *H, *bmat, *workscalar, *c = scheme->c, *a = scheme->a, *b = scheme->b, *u = scheme->u, *v = scheme->v;
    134:     PetscBLASInt rank;

    136:     PetscCall(PetscBLASIntCast(4 * ((s + 3) * 3 + 3), &lwork));
    137:     PetscCall(PetscMalloc7(PetscSqr(r), &ImV, 3 * s, &H, 3 * ss, &bmat, lwork, &workscalar, 5 * (3 + r), &workreal, r + s, &sing, r + s, &ipiv));

    139:     /* column-major input */
    140:     for (i = 0; i < r - 1; i++) {
    141:       for (j = 0; j < r - 1; j++) ImV[i + j * r] = 1.0 * (i == j) - v[(i + 1) * r + j + 1];
    142:     }
    143:     /* Build right-hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */
    144:     for (i = 1; i < r; i++) {
    145:       scheme->alpha[i] = 1. / Factorial(p + 1 - i);
    146:       for (j = 0; j < s; j++) scheme->alpha[i] -= b[i * s + j] * CPowF(c[j], p);
    147:     }
    148:     PetscCall(PetscBLASIntCast(r - 1, &m));
    149:     PetscCall(PetscBLASIntCast(r, &n));
    150:     PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&m, &one, ImV, &n, ipiv, scheme->alpha + 1, &n, &info));
    151:     PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GESV");
    152:     PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization");

    154:     /* Build right-hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */
    155:     for (i = 1; i < r; i++) {
    156:       scheme->beta[i] = 1. / Factorial(p + 2 - i) - scheme->alpha[i];
    157:       for (j = 0; j < s; j++) scheme->beta[i] -= b[i * s + j] * CPowF(c[j], p + 1);
    158:     }
    159:     PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->beta + 1, &n, &info));
    160:     PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS");
    161:     PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen");

    163:     /* Build stage_error vector
    164:            xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha;
    165:     */
    166:     for (i = 0; i < s; i++) {
    167:       scheme->stage_error[i] = CPowF(c[i], p + 1);
    168:       for (j = 0; j < s; j++) scheme->stage_error[i] -= a[i * s + j] * CPowF(c[j], p);
    169:       for (j = 1; j < r; j++) scheme->stage_error[i] += u[i * r + j] * scheme->alpha[j];
    170:     }

    172:     /* alpha[0] (epsilon in B,J,W 2007)
    173:            epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha;
    174:     */
    175:     scheme->alpha[0] = 1. / Factorial(p + 1);
    176:     for (j = 0; j < s; j++) scheme->alpha[0] -= b[0 * s + j] * CPowF(c[j], p);
    177:     for (j = 1; j < r; j++) scheme->alpha[0] += v[0 * r + j] * scheme->alpha[j];

    179:     /* right-hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */
    180:     for (i = 1; i < r; i++) {
    181:       scheme->gamma[i] = (i == 1 ? -1. : 0) * scheme->alpha[0];
    182:       for (j = 0; j < s; j++) scheme->gamma[i] += b[i * s + j] * scheme->stage_error[j];
    183:     }
    184:     PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->gamma + 1, &n, &info));
    185:     PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS");
    186:     PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen");

    188:     /* beta[0] (rho in B,J,W 2007)
    189:         e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ...
    190:             + glm.V(1,2:end)*e.beta;% - e.epsilon;
    191:     % Note: The paper (B,J,W 2007) includes the last term in their definition
    192:     * */
    193:     scheme->beta[0] = 1. / Factorial(p + 2);
    194:     for (j = 0; j < s; j++) scheme->beta[0] -= b[0 * s + j] * CPowF(c[j], p + 1);
    195:     for (j = 1; j < r; j++) scheme->beta[0] += v[0 * r + j] * scheme->beta[j];

    197:     /* gamma[0] (sigma in B,J,W 2007)
    198:     *   e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma;
    199:     * */
    200:     scheme->gamma[0] = 0.0;
    201:     for (j = 0; j < s; j++) scheme->gamma[0] += b[0 * s + j] * scheme->stage_error[j];
    202:     for (j = 1; j < r; j++) scheme->gamma[0] += v[0 * s + j] * scheme->gamma[j];

    204:     /* Assemble H
    205:     *    % " PetscInt_FMT "etermine the error estimators phi
    206:        H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ...
    207:                [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]';
    208:     % Paper has formula above without the 0, but that term must be left
    209:     % out to satisfy the conditions they propose and to make the
    210:     % example schemes work
    211:     e.H = H;
    212:     e.phi = (H \ [1 0 0;1 1 0;0 0 -1])';
    213:     e.psi = -e.phi*C;
    214:     * */
    215:     for (j = 0; j < s; j++) {
    216:       H[0 + j * 3] = CPowF(c[j], p);
    217:       H[1 + j * 3] = CPowF(c[j], p + 1);
    218:       H[2 + j * 3] = scheme->stage_error[j];
    219:       for (k = 1; k < r; k++) {
    220:         H[0 + j * 3] += CPowF(c[j], k - 1) * scheme->alpha[k];
    221:         H[1 + j * 3] += CPowF(c[j], k - 1) * scheme->beta[k];
    222:         H[2 + j * 3] -= CPowF(c[j], k - 1) * scheme->gamma[k];
    223:       }
    224:     }
    225:     bmat[0 + 0 * ss] = 1.;
    226:     bmat[0 + 1 * ss] = 0.;
    227:     bmat[0 + 2 * ss] = 0.;
    228:     bmat[1 + 0 * ss] = 1.;
    229:     bmat[1 + 1 * ss] = 1.;
    230:     bmat[1 + 2 * ss] = 0.;
    231:     bmat[2 + 0 * ss] = 0.;
    232:     bmat[2 + 1 * ss] = 0.;
    233:     bmat[2 + 2 * ss] = -1.;
    234:     m                = 3;
    235:     PetscCall(PetscBLASIntCast(s, &n));
    236:     PetscCall(PetscBLASIntCast(ss, &ldb));
    237:     rcond = 1e-12;
    238: #if defined(PETSC_USE_COMPLEX)
    239:     /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
    240:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, workreal, &info));
    241: #else
    242:     /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
    243:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, &info));
    244: #endif
    245:     PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELSS");
    246:     PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "SVD failed to converge");

    248:     for (j = 0; j < 3; j++) {
    249:       for (k = 0; k < s; k++) scheme->phi[k + j * s] = bmat[k + j * ss];
    250:     }

    252:     /* the other part of the error estimator, psi in B,J,W 2007 */
    253:     scheme->psi[0 * r + 0] = 0.;
    254:     scheme->psi[1 * r + 0] = 0.;
    255:     scheme->psi[2 * r + 0] = 0.;
    256:     for (j = 1; j < r; j++) {
    257:       scheme->psi[0 * r + j] = 0.;
    258:       scheme->psi[1 * r + j] = 0.;
    259:       scheme->psi[2 * r + j] = 0.;
    260:       for (k = 0; k < s; k++) {
    261:         scheme->psi[0 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[0 * s + k];
    262:         scheme->psi[1 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[1 * s + k];
    263:         scheme->psi[2 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[2 * s + k];
    264:       }
    265:     }
    266:     PetscCall(PetscFree7(ImV, H, bmat, workscalar, workreal, sing, ipiv));
    267:   }
    268:   /* Check which properties are satisfied */
    269:   scheme->stiffly_accurate = PETSC_TRUE;
    270:   if (scheme->c[s - 1] != 1.) scheme->stiffly_accurate = PETSC_FALSE;
    271:   for (j = 0; j < s; j++)
    272:     if (a[(s - 1) * s + j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE;
    273:   for (j = 0; j < r; j++)
    274:     if (u[(s - 1) * r + j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE;
    275:   scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */
    276:   for (j = 0; j < s - 1; j++)
    277:     if (r > 1 && b[1 * s + j] != 0.) scheme->fsal = PETSC_FALSE;
    278:   if (b[1 * s + r - 1] != 1.) scheme->fsal = PETSC_FALSE;
    279:   for (j = 0; j < r; j++)
    280:     if (r > 1 && v[1 * r + j] != 0.) scheme->fsal = PETSC_FALSE;

    282:   *inscheme = scheme;
    283:   PetscFunctionReturn(PETSC_SUCCESS);
    284: }

    286: static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc)
    287: {
    288:   PetscFunctionBegin;
    289:   PetscCall(PetscFree5(sc->c, sc->a, sc->b, sc->u, sc->v));
    290:   PetscCall(PetscFree6(sc->alpha, sc->beta, sc->gamma, sc->phi, sc->psi, sc->stage_error));
    291:   PetscCall(PetscFree(sc));
    292:   PetscFunctionReturn(PETSC_SUCCESS);
    293: }

    295: static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
    296: {
    297:   PetscInt i;

    299:   PetscFunctionBegin;
    300:   for (i = 0; i < gl->nschemes; i++) {
    301:     if (gl->schemes[i]) PetscCall(TSGLLESchemeDestroy(gl->schemes[i]));
    302:   }
    303:   PetscCall(PetscFree(gl->schemes));
    304:   gl->nschemes = 0;
    305:   PetscCall(PetscMemzero(gl->type_name, sizeof(gl->type_name)));
    306:   PetscFunctionReturn(PETSC_SUCCESS);
    307: }

    309: static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer, PetscInt m, PetscInt n, const PetscScalar a[], const char name[])
    310: {
    311:   PetscBool isascii;
    312:   PetscInt  i, j;

    314:   PetscFunctionBegin;
    315:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    316:   if (isascii) {
    317:     PetscCall(PetscViewerASCIIPrintf(viewer, "%30s = [", name));
    318:     for (i = 0; i < m; i++) {
    319:       if (i) PetscCall(PetscViewerASCIIPrintf(viewer, "%30s   [", ""));
    320:       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
    321:       for (j = 0; j < n; j++) PetscCall(PetscViewerASCIIPrintf(viewer, " %12.8g", (double)PetscRealPart(a[i * n + j])));
    322:       PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
    323:       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    324:     }
    325:   }
    326:   PetscFunctionReturn(PETSC_SUCCESS);
    327: }

    329: static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer)
    330: {
    331:   PetscBool isascii;

    333:   PetscFunctionBegin;
    334:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    335:   PetscCheck(isascii, PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported", ((PetscObject)viewer)->type_name);
    336:   PetscCall(PetscViewerASCIIPrintf(viewer, "GL scheme p,q,r,s = %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "\n", sc->p, sc->q, sc->r, sc->s));
    337:   PetscCall(PetscViewerASCIIPushTab(viewer));
    338:   PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s,  FSAL: %s\n", sc->stiffly_accurate ? "yes" : "no", sc->fsal ? "yes" : "no"));
    339:   PetscCall(PetscViewerASCIIPrintf(viewer, "Leading error constants: %10.3e  %10.3e  %10.3e\n", (double)PetscRealPart(sc->alpha[0]), (double)PetscRealPart(sc->beta[0]), (double)PetscRealPart(sc->gamma[0])));
    340:   PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->c, "Abscissas c"));
    341:   if (view_details) {
    342:     PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->s, sc->a, "A"));
    343:     PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->s, sc->b, "B"));
    344:     PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->r, sc->u, "U"));
    345:     PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->r, sc->v, "V"));

    347:     PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->s, sc->phi, "Error estimate phi"));
    348:     PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->r, sc->psi, "Error estimate psi"));
    349:     PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->alpha, "Modify alpha"));
    350:     PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->beta, "Modify beta"));
    351:     PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->gamma, "Modify gamma"));
    352:     PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->stage_error, "Stage error xi"));
    353:   }
    354:   PetscCall(PetscViewerASCIIPopTab(viewer));
    355:   PetscFunctionReturn(PETSC_SUCCESS);
    356: }

    358: static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[])
    359: {
    360:   PetscInt i;

    362:   PetscFunctionBegin;
    363:   PetscCheck(sc->r <= 64 && sc->s <= 64, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Ridiculous number of stages or items passed between stages");
    364:   /* build error vectors*/
    365:   for (i = 0; i < 3; i++) {
    366:     PetscScalar phih[64];
    367:     PetscInt    j;
    368:     for (j = 0; j < sc->s; j++) phih[j] = sc->phi[i * sc->s + j] * h;
    369:     PetscCall(VecZeroEntries(hm[i]));
    370:     PetscCall(VecMAXPY(hm[i], sc->s, phih, Ydot));
    371:     PetscCall(VecMAXPY(hm[i], sc->r, &sc->psi[i * sc->r], Xold));
    372:   }
    373:   PetscFunctionReturn(PETSC_SUCCESS);
    374: }

    376: static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
    377: {
    378:   PetscScalar brow[32], vrow[32];
    379:   PetscInt    i, j, r, s;

    381:   PetscFunctionBegin;
    382:   /* Build the new solution from (X,Ydot) */
    383:   r = sc->r;
    384:   s = sc->s;
    385:   for (i = 0; i < r; i++) {
    386:     PetscCall(VecZeroEntries(X[i]));
    387:     for (j = 0; j < s; j++) brow[j] = h * sc->b[i * s + j];
    388:     PetscCall(VecMAXPY(X[i], s, brow, Ydot));
    389:     for (j = 0; j < r; j++) vrow[j] = sc->v[i * r + j];
    390:     PetscCall(VecMAXPY(X[i], r, vrow, Xold));
    391:   }
    392:   PetscFunctionReturn(PETSC_SUCCESS);
    393: }

    395: static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
    396: {
    397:   PetscScalar brow[32], vrow[32];
    398:   PetscReal   ratio;
    399:   PetscInt    i, j, p, r, s;

    401:   PetscFunctionBegin;
    402:   /* Build the new solution from (X,Ydot) */
    403:   p     = sc->p;
    404:   r     = sc->r;
    405:   s     = sc->s;
    406:   ratio = next_h / h;
    407:   for (i = 0; i < r; i++) {
    408:     PetscCall(VecZeroEntries(X[i]));
    409:     for (j = 0; j < s; j++) {
    410:       brow[j] = h * (PetscPowRealInt(ratio, i) * sc->b[i * s + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->phi[0 * s + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->phi[1 * s + j] + sc->gamma[i] * sc->phi[2 * s + j]));
    411:     }
    412:     PetscCall(VecMAXPY(X[i], s, brow, Ydot));
    413:     for (j = 0; j < r; j++) {
    414:       vrow[j] = (PetscPowRealInt(ratio, i) * sc->v[i * r + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->psi[0 * r + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->psi[1 * r + j] + sc->gamma[i] * sc->psi[2 * r + j]));
    415:     }
    416:     PetscCall(VecMAXPY(X[i], r, vrow, Xold));
    417:   }
    418:   if (r < next_sc->r) {
    419:     PetscCheck(r + 1 == next_sc->r, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot accommodate jump in r greater than 1");
    420:     PetscCall(VecZeroEntries(X[r]));
    421:     for (j = 0; j < s; j++) brow[j] = h * PetscPowRealInt(ratio, p + 1) * sc->phi[0 * s + j];
    422:     PetscCall(VecMAXPY(X[r], s, brow, Ydot));
    423:     for (j = 0; j < r; j++) vrow[j] = PetscPowRealInt(ratio, p + 1) * sc->psi[0 * r + j];
    424:     PetscCall(VecMAXPY(X[r], r, vrow, Xold));
    425:   }
    426:   PetscFunctionReturn(PETSC_SUCCESS);
    427: }

    429: static PetscErrorCode TSGLLECreate_IRKS(TS ts)
    430: {
    431:   TS_GLLE *gl = (TS_GLLE *)ts->data;

    433:   PetscFunctionBegin;
    434:   gl->Destroy               = TSGLLEDestroy_Default;
    435:   gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default;
    436:   gl->CompleteStep          = TSGLLECompleteStep_RescaleAndModify;
    437:   PetscCall(PetscMalloc1(10, &gl->schemes));
    438:   gl->nschemes = 0;

    440:   {
    441:     /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3
    442:     * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE.
    443:     * irks(0.3,0,[.3,1],[1],1)
    444:     * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2)
    445:     * but doing so would sacrifice the error estimator.
    446:     */
    447:     const PetscScalar c[2]    = {3. / 10., 1.};
    448:     const PetscScalar a[2][2] = {
    449:       {3. / 10., 0       },
    450:       {7. / 10., 3. / 10.}
    451:     };
    452:     const PetscScalar b[2][2] = {
    453:       {7. / 10., 3. / 10.},
    454:       {0,        1       }
    455:     };
    456:     const PetscScalar u[2][2] = {
    457:       {1, 0},
    458:       {1, 0}
    459:     };
    460:     const PetscScalar v[2][2] = {
    461:       {1, 0},
    462:       {0, 0}
    463:     };
    464:     PetscCall(TSGLLESchemeCreate(1, 1, 2, 2, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
    465:   }

    467:   {
    468:     /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */
    469:     /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */
    470:     const PetscScalar c[3]    = {1. / 3., 2. / 3., 1};
    471:     const PetscScalar a[3][3] = {
    472:       {4. / 9.,              0,                     0      },
    473:       {1.03750643704090e+00, 4. / 9.,               0      },
    474:       {7.67024779410304e-01, -3.81140216918943e-01, 4. / 9.}
    475:     };
    476:     const PetscScalar b[3][3] = {
    477:       {0.767024779410304,  -0.381140216918943, 4. / 9.          },
    478:       {0.000000000000000,  0.000000000000000,  1.000000000000000},
    479:       {-2.075048385225385, 0.621728385225383,  1.277197204924873}
    480:     };
    481:     const PetscScalar u[3][3] = {
    482:       {1.0000000000000000, -0.1111111111111109, -0.0925925925925922},
    483:       {1.0000000000000000, -0.8152842148186744, -0.4199095530877056},
    484:       {1.0000000000000000, 0.1696709930641948,  0.0539741070314165 }
    485:     };
    486:     const PetscScalar v[3][3] = {
    487:       {1.0000000000000000, 0.1696709930641948, 0.0539741070314165},
    488:       {0.000000000000000,  0.000000000000000,  0.000000000000000 },
    489:       {0.000000000000000,  0.176122795075129,  0.000000000000000 }
    490:     };
    491:     PetscCall(TSGLLESchemeCreate(2, 2, 3, 3, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
    492:   }
    493:   {
    494:     /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */
    495:     const PetscScalar c[4]    = {0.25, 0.5, 0.75, 1.0};
    496:     const PetscScalar a[4][4] = {
    497:       {9. / 40.,             0,                     0,                 0       },
    498:       {2.11286958887701e-01, 9. / 40.,              0,                 0       },
    499:       {9.46338294287584e-01, -3.42942861246094e-01, 9. / 40.,          0       },
    500:       {0.521490453970721,    -0.662474225622980,    0.490476425459734, 9. / 40.}
    501:     };
    502:     const PetscScalar b[4][4] = {
    503:       {0.521490453970721,  -0.662474225622980, 0.490476425459734,  9. / 40.         },
    504:       {0.000000000000000,  0.000000000000000,  0.000000000000000,  1.000000000000000},
    505:       {-0.084677029310348, 1.390757514776085,  -1.568157386206001, 2.023192696767826},
    506:       {0.465383797936408,  1.478273530625148,  -1.930836081010182, 1.644872111193354}
    507:     };
    508:     const PetscScalar u[4][4] = {
    509:       {1.00000000000000000, 0.02500000000001035,  -0.02499999999999053, -0.00442708333332865},
    510:       {1.00000000000000000, 0.06371304111232945,  -0.04032173972189845, -0.01389438413189452},
    511:       {1.00000000000000000, -0.07839543304147778, 0.04738685705116663,  0.02032603595928376 },
    512:       {1.00000000000000000, 0.42550734619251651,  0.10800718022400080,  -0.01726712647760034}
    513:     };
    514:     const PetscScalar v[4][4] = {
    515:       {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034},
    516:       {0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000   },
    517:       {0.000000000000000,   -1.761115796027561,  -0.521284157173780,  0.258249384305463   },
    518:       {0.000000000000000,   -1.657693358744728,  -1.052227765232394,  0.521284157173780   }
    519:     };
    520:     PetscCall(TSGLLESchemeCreate(3, 3, 4, 4, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
    521:   }
    522:   {
    523:     /* p=q=4, r=s=5:
    524:           irks(3/11,0,[1:5]/5, [0.1715   -0.1238    0.6617],...
    525:           [ -0.0812    0.4079    1.0000
    526:              1.0000         0         0
    527:              0.8270    1.0000         0])
    528:     */
    529:     const PetscScalar c[5]    = {0.2, 0.4, 0.6, 0.8, 1.0};
    530:     const PetscScalar a[5][5] = {
    531:       {2.72727272727352e-01,  0.00000000000000e+00,  0.00000000000000e+00,  0.00000000000000e+00, 0.00000000000000e+00},
    532:       {-1.03980153733431e-01, 2.72727272727405e-01,  0.00000000000000e+00,  0.00000000000000e+00, 0.00000000000000e+00},
    533:       {-1.58615400341492e+00, 7.44168951881122e-01,  2.72727272727309e-01,  0.00000000000000e+00, 0.00000000000000e+00},
    534:       {-8.73658042865628e-01, 5.37884671894595e-01,  -1.63298538799523e-01, 2.72727272726996e-01, 0.00000000000000e+00},
    535:       {2.95489397443992e-01,  -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01}
    536:     };
    537:     const PetscScalar b[5][5] = {
    538:       {2.95489397443992e-01,  -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00,  2.72727272727288e-01},
    539:       {0.00000000000000e+00,  1.11022302462516e-16,  -2.22044604925031e-16, 0.00000000000000e+00,  1.00000000000000e+00},
    540:       {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00,  6.32331093108427e-01},
    541:       {8.35690179937017e+00,  -2.26640927349732e+00, 6.86647884973826e+00,  -5.22595158025740e+00, 4.50893068837431e+00},
    542:       {1.27656267027479e+01,  2.80882153840821e+00,  8.91173096522890e+00,  -1.07936444078906e+01, 4.82534148988854e+00}
    543:     };
    544:     const PetscScalar u[5][5] = {
    545:       {1.00000000000000e+00, -7.27272727273551e-02, -3.45454545454419e-02, -4.12121212119565e-03, -2.96969696964014e-04},
    546:       {1.00000000000000e+00, 2.31252881006154e-01,  -8.29487834416481e-03, -9.07191207681020e-03, -1.70378403743473e-03},
    547:       {1.00000000000000e+00, 1.16925777880663e+00,  3.59268562942635e-02,  -4.09013451730615e-02, -1.02411119670164e-02},
    548:       {1.00000000000000e+00, 1.02634463704356e+00,  1.59375044913405e-01,  1.89673015035370e-03,  -4.89987231897569e-03},
    549:       {1.00000000000000e+00, 1.27746320298021e+00,  2.37186008132728e-01,  -8.28694373940065e-02, -5.34396510196430e-02}
    550:     };
    551:     const PetscScalar v[5][5] = {
    552:       {1.00000000000000e+00, 1.27746320298021e+00,  2.37186008132728e-01,  -8.28694373940065e-02, -5.34396510196430e-02},
    553:       {0.00000000000000e+00, -1.77635683940025e-15, -1.99840144432528e-15, -9.99200722162641e-16, -3.33066907387547e-16},
    554:       {0.00000000000000e+00, 4.37280081906924e+00,  5.49221645016377e-02,  -8.88913177394943e-02, 1.12879077989154e-01 },
    555:       {0.00000000000000e+00, -1.22399504837280e+01, -5.21287338448645e+00, -8.03952325565291e-01, 4.60298678047147e-01 },
    556:       {0.00000000000000e+00, -1.85178762883829e+01, -5.21411849862624e+00, -1.04283436528809e+00, 7.49030161063651e-01 }
    557:     };
    558:     PetscCall(TSGLLESchemeCreate(4, 4, 5, 5, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
    559:   }
    560:   {
    561:     /* p=q=5, r=s=6;
    562:        irks(1/3,0,[1:6]/6,...
    563:           [-0.0489    0.4228   -0.8814    0.9021],...
    564:           [-0.3474   -0.6617    0.6294    0.2129
    565:             0.0044   -0.4256   -0.1427   -0.8936
    566:            -0.8267    0.4821    0.1371   -0.2557
    567:            -0.4426   -0.3855   -0.7514    0.3014])
    568:     */
    569:     const PetscScalar c[6]    = {1. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1.};
    570:     const PetscScalar a[6][6] = {
    571:       {3.33333333333940e-01,  0,                     0,                     0,                     0,                    0                   },
    572:       {-8.64423857333350e-02, 3.33333333332888e-01,  0,                     0,                     0,                    0                   },
    573:       {-2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01,  0,                     0,                    0                   },
    574:       {-4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01,  0,                    0                   },
    575:       {-6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01,  -4.48352364517632e-01, 3.33333333328483e-01, 0                   },
    576:       {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00,  -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}
    577:     };
    578:     const PetscScalar b[6][6] = {
    579:       {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00,  -2.23969848002481e+00, 6.62807710124007e-01,  3.33333333335440e-01 },
    580:       {-8.88178419700125e-16, 4.44089209850063e-16,  -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00,  1.00000000000001e+00 },
    581:       {-2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01,  2.56943874812797e+01,  -3.06702268304488e+01, 6.68067773510103e+00 },
    582:       {5.47971245256474e+01,  6.80366875868284e+01,  -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01,  -1.17819043489036e+01},
    583:       {-2.33332114788869e+02, 6.12942539462634e+01,  -4.91850135865944e+01, 1.82716844135480e+02,  -1.29788173979395e+02, 3.09968095651099e+01 },
    584:       {-1.72049132343751e+02, 8.60194713593999e+00,  7.98154219170200e-01,  1.50371386053218e+02,  -1.18515423962066e+02, 2.50898277784663e+01 }
    585:     };
    586:     const PetscScalar u[6][6] = {
    587:       {1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06},
    588:       {1.00000000000000e+00, 8.64423857327162e-02,  -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04},
    589:       {1.00000000000000e+00, 4.57135912953434e+00,  1.06514719719137e+00,  1.33517564218007e-01,  1.11365952968659e-02,  6.12382756769504e-04 },
    590:       {1.00000000000000e+00, 9.23391519753404e+00,  2.22431212392095e+00,  2.91823807741891e-01,  2.52058456411084e-02,  1.22800542949647e-03 },
    591:       {1.00000000000000e+00, 1.48175480533865e+01,  3.73439117461835e+00,  5.14648336541804e-01,  4.76430038853402e-02,  2.56798515502156e-03 },
    592:       {1.00000000000000e+00, 1.50512347758335e+01,  4.10099701165164e+00,  5.66039141003603e-01,  3.91213893800891e-02,  -2.99136269067853e-03}
    593:     };
    594:     const PetscScalar v[6][6] = {
    595:       {1.00000000000000e+00, 1.50512347758335e+01,  4.10099701165164e+00,  5.66039141003603e-01,  3.91213893800891e-02,  -2.99136269067853e-03},
    596:       {0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16},
    597:       {0.00000000000000e+00, 1.22250171233141e+01,  -1.77150760606169e+00, 3.54516769879390e-01,  6.22298845883398e-01,  2.31647447450276e-01 },
    598:       {0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01,  6.55727990241799e-02,  1.63175368287079e-01 },
    599:       {0.00000000000000e+00, 1.37297394708005e+02,  -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01,  9.16629423682464e-01 },
    600:       {0.00000000000000e+00, 1.05703241119022e+02,  -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00 }
    601:     };
    602:     PetscCall(TSGLLESchemeCreate(5, 5, 6, 6, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
    603:   }
    604:   PetscFunctionReturn(PETSC_SUCCESS);
    605: }

    607: /*@
    608:   TSGLLESetType - sets the class of general linear method, `TSGLLE` to use for time-stepping

    610:   Collective

    612:   Input Parameters:
    613: + ts   - the `TS` context
    614: - type - a method, currently only `TSGLLE_IRKS` is available

    616:   Options Database Key:
    617: . -ts_gl_type (irks) - sets the method

    619:   Level: intermediate

    621: .seealso: [](ch_ts), `TS`, `TSGLLEType`, `TSGLLE`, `TSGLLERegister()`, `TSGLLE_IRKS`, `TSGLLEGetAcceptType()`,
    622:           `TSGLLESetAcceptType()`, `TSGLLEAcceptType()`
    623: @*/
    624: PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type)
    625: {
    626:   PetscFunctionBegin;
    628:   PetscAssertPointer(type, 2);
    629:   PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type));
    630:   PetscFunctionReturn(PETSC_SUCCESS);
    631: }

    633: /*@C
    634:   TSGLLESetAcceptType - sets the acceptance test for `TSGLLE`

    636:   Logically Collective

    638:   Input Parameters:
    639: + ts   - the `TS` context
    640: - type - the type

    642:   Options Database Key:
    643: . -ts_gl_accept_type (always) - sets the method used to determine whether to accept or reject a step

    645:   Level: intermediate

    647:   Notes:
    648:   Time integrators that need to control error must have the option to reject a time step based
    649:   on local error estimates. This function allows different schemes to be set.

    651: .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt`
    652: @*/
    653: PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type)
    654: {
    655:   PetscFunctionBegin;
    657:   PetscAssertPointer(type, 2);
    658:   PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type));
    659:   PetscFunctionReturn(PETSC_SUCCESS);
    660: }

    662: /*@
    663:   TSGLLEGetAdapt - gets the `TSGLLEAdapt` object from the `TS`

    665:   Not Collective

    667:   Input Parameter:
    668: . ts - the `TS` context

    670:   Output Parameter:
    671: . adapt - the `TSGLLEAdapt` context

    673:   Level: advanced

    675:   Note:
    676:   This allows the user set options on the `TSGLLEAdapt` object. Usually it is better to do this
    677:   using the options database, so this function is rarely needed.

    679: .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegister()`
    680: @*/
    681: PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt)
    682: {
    683:   PetscFunctionBegin;
    685:   PetscAssertPointer(adapt, 2);
    686:   PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt));
    687:   PetscFunctionReturn(PETSC_SUCCESS);
    688: }

    690: static PetscErrorCode TSGLLEAccept_Always(TS ts, PetscReal tleft, PetscReal h, const PetscReal enorms[], PetscBool *accept)
    691: {
    692:   PetscFunctionBegin;
    693:   *accept = PETSC_TRUE;
    694:   PetscFunctionReturn(PETSC_SUCCESS);
    695: }

    697: static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
    698: {
    699:   TS_GLLE     *gl = (TS_GLLE *)ts->data;
    700:   PetscScalar *x, *w;
    701:   PetscInt     n, i;

    703:   PetscFunctionBegin;
    704:   PetscCall(VecGetArray(gl->X[0], &x));
    705:   PetscCall(VecGetArray(gl->W, &w));
    706:   PetscCall(VecGetLocalSize(gl->W, &n));
    707:   for (i = 0; i < n; i++) w[i] = 1. / (gl->wrms_atol + gl->wrms_rtol * PetscAbsScalar(x[i]));
    708:   PetscCall(VecRestoreArray(gl->X[0], &x));
    709:   PetscCall(VecRestoreArray(gl->W, &w));
    710:   PetscFunctionReturn(PETSC_SUCCESS);
    711: }

    713: static PetscErrorCode TSGLLEVecNormWRMS(TS ts, Vec X, PetscReal *nrm)
    714: {
    715:   TS_GLLE     *gl = (TS_GLLE *)ts->data;
    716:   PetscScalar *x, *w;
    717:   PetscReal    sum = 0.0, gsum;
    718:   PetscInt     n, N, i;

    720:   PetscFunctionBegin;
    721:   PetscCall(VecGetArray(X, &x));
    722:   PetscCall(VecGetArray(gl->W, &w));
    723:   PetscCall(VecGetLocalSize(gl->W, &n));
    724:   for (i = 0; i < n; i++) sum += PetscAbsScalar(PetscSqr(x[i] * w[i]));
    725:   PetscCall(VecRestoreArray(X, &x));
    726:   PetscCall(VecRestoreArray(gl->W, &w));
    727:   PetscCallMPI(MPIU_Allreduce(&sum, &gsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
    728:   PetscCall(VecGetSize(gl->W, &N));
    729:   *nrm = PetscSqrtReal(gsum / (1. * N));
    730:   PetscFunctionReturn(PETSC_SUCCESS);
    731: }

    733: static PetscErrorCode TSGLLESetType_GLLE(TS ts, TSGLLEType type)
    734: {
    735:   PetscBool same;
    736:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
    737:   PetscErrorCode (*r)(TS);

    739:   PetscFunctionBegin;
    740:   if (gl->type_name[0]) {
    741:     PetscCall(PetscStrcmp(gl->type_name, type, &same));
    742:     if (same) PetscFunctionReturn(PETSC_SUCCESS);
    743:     PetscCall((*gl->Destroy)(gl));
    744:   }

    746:   PetscCall(PetscFunctionListFind(TSGLLEList, type, &r));
    747:   PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLE type \"%s\" given", type);
    748:   PetscCall((*r)(ts));
    749:   PetscCall(PetscStrncpy(gl->type_name, type, sizeof(gl->type_name)));
    750:   PetscFunctionReturn(PETSC_SUCCESS);
    751: }

    753: static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type)
    754: {
    755:   TSGLLEAcceptFn *r;
    756:   TS_GLLE        *gl = (TS_GLLE *)ts->data;

    758:   PetscFunctionBegin;
    759:   PetscCall(PetscFunctionListFind(TSGLLEAcceptList, type, &r));
    760:   PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAccept type \"%s\" given", type);
    761:   gl->Accept = r;
    762:   PetscCall(PetscStrncpy(gl->accept_name, type, sizeof(gl->accept_name)));
    763:   PetscFunctionReturn(PETSC_SUCCESS);
    764: }

    766: static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt)
    767: {
    768:   TS_GLLE *gl = (TS_GLLE *)ts->data;

    770:   PetscFunctionBegin;
    771:   if (!gl->adapt) {
    772:     PetscCall(TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts), &gl->adapt));
    773:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)gl->adapt, (PetscObject)ts, 1));
    774:   }
    775:   *adapt = gl->adapt;
    776:   PetscFunctionReturn(PETSC_SUCCESS);
    777: }

    779: static PetscErrorCode TSGLLEChooseNextScheme(TS ts, PetscReal h, const PetscReal hmnorm[], PetscInt *next_scheme, PetscReal *next_h, PetscBool *finish)
    780: {
    781:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
    782:   PetscInt  i, n, cur_p, cur, next_sc, candidates[64], orders[64];
    783:   PetscReal errors[64], costs[64], tleft;

    785:   PetscFunctionBegin;
    786:   cur   = -1;
    787:   cur_p = gl->schemes[gl->current_scheme]->p;
    788:   tleft = ts->max_time - (ts->ptime + ts->time_step);
    789:   for (i = 0, n = 0; i < gl->nschemes; i++) {
    790:     TSGLLEScheme sc = gl->schemes[i];
    791:     if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
    792:     if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[0];
    793:     else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[1];
    794:     else if (sc->p == cur_p + 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * (hmnorm[2] + hmnorm[3]);
    795:     else continue;
    796:     candidates[n] = i;
    797:     orders[n]     = PetscMin(sc->p, sc->q); /* order of global truncation error */
    798:     costs[n]      = sc->s;                  /* estimate the cost as the number of stages */
    799:     if (i == gl->current_scheme) cur = n;
    800:     n++;
    801:   }
    802:   PetscCheck(cur >= 0 && gl->nschemes > cur, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Current scheme not found in scheme list");
    803:   PetscCall(TSGLLEAdaptChoose(gl->adapt, n, orders, errors, costs, cur, h, tleft, &next_sc, next_h, finish));
    804:   *next_scheme = candidates[next_sc];
    805:   PetscCall(PetscInfo(ts, "Adapt chose scheme %" PetscInt_FMT " (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ") with step size %6.2e, finish=%s\n", *next_scheme, gl->schemes[*next_scheme]->p, gl->schemes[*next_scheme]->q,
    806:                       gl->schemes[*next_scheme]->r, gl->schemes[*next_scheme]->s, (double)*next_h, PetscBools[*finish]));
    807:   PetscFunctionReturn(PETSC_SUCCESS);
    808: }

    810: static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s)
    811: {
    812:   TS_GLLE *gl = (TS_GLLE *)ts->data;

    814:   PetscFunctionBegin;
    815:   *max_r = gl->schemes[gl->nschemes - 1]->r;
    816:   *max_s = gl->schemes[gl->nschemes - 1]->s;
    817:   PetscFunctionReturn(PETSC_SUCCESS);
    818: }

    820: static PetscErrorCode TSSolve_GLLE(TS ts)
    821: {
    822:   TS_GLLE            *gl = (TS_GLLE *)ts->data;
    823:   PetscInt            i, k, its, lits, max_r, max_s;
    824:   PetscBool           final_step, finish;
    825:   SNESConvergedReason snesreason;

    827:   PetscFunctionBegin;
    828:   PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

    830:   PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
    831:   PetscCall(VecCopy(ts->vec_sol, gl->X[0]));
    832:   for (i = 1; i < max_r; i++) PetscCall(VecZeroEntries(gl->X[i]));
    833:   PetscCall(TSGLLEUpdateWRMS(ts));

    835:   if (0) {
    836:     /* Find consistent initial data for DAE */
    837:     gl->stage_time = ts->ptime + ts->time_step;
    838:     gl->scoeff     = 1.;
    839:     gl->stage      = 0;

    841:     PetscCall(VecCopy(ts->vec_sol, gl->Z));
    842:     PetscCall(VecCopy(ts->vec_sol, gl->Y));
    843:     PetscCall(SNESSolve(ts->snes, NULL, gl->Y));
    844:     PetscCall(SNESGetIterationNumber(ts->snes, &its));
    845:     PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
    846:     PetscCall(SNESGetConvergedReason(ts->snes, &snesreason));

    848:     ts->snes_its += its;
    849:     ts->ksp_its += lits;
    850:     if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) {
    851:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
    852:       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures));
    853:       PetscFunctionReturn(PETSC_SUCCESS);
    854:     }
    855:   }

    857:   PetscCheck(gl->current_scheme >= 0, PETSC_COMM_SELF, PETSC_ERR_ORDER, "A starting scheme has not been provided");

    859:   for (k = 0, final_step = PETSC_FALSE, finish = PETSC_FALSE; k < ts->max_steps && !finish; k++) {
    860:     PetscInt           j, r, s, next_scheme = 0, rejections;
    861:     PetscReal          h, hmnorm[4], enorm[3], next_h;
    862:     PetscBool          accept;
    863:     const PetscScalar *c, *a, *u;
    864:     Vec               *X, *Ydot, Y;
    865:     TSGLLEScheme       scheme = gl->schemes[gl->current_scheme];

    867:     r    = scheme->r;
    868:     s    = scheme->s;
    869:     c    = scheme->c;
    870:     a    = scheme->a;
    871:     u    = scheme->u;
    872:     h    = ts->time_step;
    873:     X    = gl->X;
    874:     Ydot = gl->Ydot;
    875:     Y    = gl->Y;

    877:     if (ts->ptime > ts->max_time) break;

    879:     /*
    880:       We only call PreStep at the start of each STEP, not each STAGE.  This is because it is
    881:       possible to fail (have to restart a step) after multiple stages.
    882:     */
    883:     PetscCall(TSPreStep(ts));

    885:     rejections = 0;
    886:     while (1) {
    887:       for (i = 0; i < s; i++) {
    888:         PetscScalar shift;
    889:         gl->scoeff     = 1. / PetscRealPart(a[i * s + i]);
    890:         shift          = gl->scoeff / ts->time_step;
    891:         gl->stage      = i;
    892:         gl->stage_time = ts->ptime + PetscRealPart(c[i]) * h;

    894:         /*
    895:         * Stage equation: Y = h A Y' + U X
    896:         * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
    897:         * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
    898:         * Then y'_i = z + 1/(h a_ii) y_i
    899:         */
    900:         PetscCall(VecZeroEntries(gl->Z));
    901:         for (j = 0; j < r; j++) PetscCall(VecAXPY(gl->Z, -shift * u[i * r + j], X[j]));
    902:         for (j = 0; j < i; j++) PetscCall(VecAXPY(gl->Z, -shift * h * a[i * s + j], Ydot[j]));
    903:         /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */

    905:         /* Compute an estimate of Y to start Newton iteration */
    906:         if (gl->extrapolate) {
    907:           if (i == 0) {
    908:             /* Linear extrapolation on the first stage */
    909:             PetscCall(VecWAXPY(Y, c[i] * h, X[1], X[0]));
    910:           } else {
    911:             /* Linear extrapolation from the last stage */
    912:             PetscCall(VecAXPY(Y, (c[i] - c[i - 1]) * h, Ydot[i - 1]));
    913:           }
    914:         } else if (i == 0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
    915:           PetscCall(VecCopy(X[0], Y));
    916:         }

    918:         /* Solve this stage (Ydot[i] is computed during function evaluation) */
    919:         PetscCall(SNESSolve(ts->snes, NULL, Y));
    920:         PetscCall(SNESGetIterationNumber(ts->snes, &its));
    921:         PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
    922:         PetscCall(SNESGetConvergedReason(ts->snes, &snesreason));
    923:         ts->snes_its += its;
    924:         ts->ksp_its += lits;
    925:         if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) {
    926:           ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
    927:           PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures));
    928:           PetscFunctionReturn(PETSC_SUCCESS);
    929:         }
    930:       }

    932:       gl->stage_time = ts->ptime + ts->time_step;

    934:       PetscCall((*gl->EstimateHigherMoments)(scheme, h, Ydot, gl->X, gl->himom));
    935:       /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */
    936:       for (i = 0; i < 3; i++) PetscCall(TSGLLEVecNormWRMS(ts, gl->himom[i], &hmnorm[i + 1]));
    937:       enorm[0] = PetscRealPart(scheme->alpha[0]) * hmnorm[1];
    938:       enorm[1] = PetscRealPart(scheme->beta[0]) * hmnorm[2];
    939:       enorm[2] = PetscRealPart(scheme->gamma[0]) * hmnorm[3];
    940:       PetscCall((*gl->Accept)(ts, ts->max_time - gl->stage_time, h, enorm, &accept));
    941:       if (accept) goto accepted;
    942:       rejections++;
    943:       PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " (t=%g) not accepted, rejections=%" PetscInt_FMT "\n", k, (double)gl->stage_time, rejections));
    944:       if (rejections > gl->max_step_rejections) break;
    945:       /*
    946:         There are lots of reasons why a step might be rejected, including solvers not converging and other factors that
    947:         TSGLLEChooseNextScheme does not support.  Additionally, the error estimates may be very screwed up, so I'm not
    948:         convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor
    949:         (the adaptor interface probably has to change).  Here we make an arbitrary and naive choice.  This assumes that
    950:         steps were written in Nordsieck form.  The "correct" method would be to re-complete the previous time step with
    951:         the correct "next" step size.  It is unclear to me whether the present ad-hoc method of rescaling X is stable.
    952:       */
    953:       h *= 0.5;
    954:       for (i = 1; i < scheme->r; i++) PetscCall(VecScale(X[i], PetscPowRealInt(0.5, i)));
    955:     }
    956:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Time step %" PetscInt_FMT " (t=%g) not accepted after %" PetscInt_FMT " failures", k, (double)gl->stage_time, rejections);

    958:   accepted:
    959:     /* This term is not error, but it *would* be the leading term for a lower order method */
    960:     PetscCall(TSGLLEVecNormWRMS(ts, gl->X[scheme->r - 1], &hmnorm[0]));
    961:     /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */

    963:     PetscCall(PetscInfo(ts, "Last moment norm %10.2e, estimated error norms %10.2e %10.2e %10.2e\n", (double)hmnorm[0], (double)enorm[0], (double)enorm[1], (double)enorm[2]));
    964:     if (!final_step) {
    965:       PetscCall(TSGLLEChooseNextScheme(ts, h, hmnorm, &next_scheme, &next_h, &final_step));
    966:     } else {
    967:       /* Dummy values to complete the current step in a consistent manner */
    968:       next_scheme = gl->current_scheme;
    969:       next_h      = h;
    970:       finish      = PETSC_TRUE;
    971:     }

    973:     X        = gl->Xold;
    974:     gl->Xold = gl->X;
    975:     gl->X    = X;
    976:     PetscCall((*gl->CompleteStep)(scheme, h, gl->schemes[next_scheme], next_h, Ydot, gl->Xold, gl->X));

    978:     PetscCall(TSGLLEUpdateWRMS(ts));

    980:     /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */
    981:     PetscCall(VecCopy(gl->X[0], ts->vec_sol));
    982:     ts->ptime += h;
    983:     ts->steps++;

    985:     PetscCall(TSPostEvaluate(ts));
    986:     PetscCall(TSPostStep(ts));
    987:     PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

    989:     gl->current_scheme = next_scheme;
    990:     ts->time_step      = next_h;
    991:   }
    992:   PetscFunctionReturn(PETSC_SUCCESS);
    993: }

    995: static PetscErrorCode TSReset_GLLE(TS ts)
    996: {
    997:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    998:   PetscInt max_r, max_s;

    1000:   PetscFunctionBegin;
    1001:   if (gl->setupcalled) {
    1002:     PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
    1003:     PetscCall(VecDestroyVecs(max_r, &gl->Xold));
    1004:     PetscCall(VecDestroyVecs(max_r, &gl->X));
    1005:     PetscCall(VecDestroyVecs(max_s, &gl->Ydot));
    1006:     PetscCall(VecDestroyVecs(3, &gl->himom));
    1007:     PetscCall(VecDestroy(&gl->W));
    1008:     PetscCall(VecDestroy(&gl->Y));
    1009:     PetscCall(VecDestroy(&gl->Z));
    1010:   }
    1011:   gl->setupcalled = PETSC_FALSE;
    1012:   PetscFunctionReturn(PETSC_SUCCESS);
    1013: }

    1015: static PetscErrorCode TSDestroy_GLLE(TS ts)
    1016: {
    1017:   TS_GLLE *gl = (TS_GLLE *)ts->data;

    1019:   PetscFunctionBegin;
    1020:   PetscCall(TSReset_GLLE(ts));
    1021:   if (ts->dm) {
    1022:     PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
    1023:     PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
    1024:   }
    1025:   PetscCall(TSGLLEAdaptDestroy(&gl->adapt));
    1026:   if (gl->Destroy) PetscCall((*gl->Destroy)(gl));
    1027:   PetscCall(PetscFree(ts->data));
    1028:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL));
    1029:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL));
    1030:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL));
    1031:   PetscFunctionReturn(PETSC_SUCCESS);
    1032: }

    1034: /*
    1035:     This defines the nonlinear equation that is to be solved with SNES
    1036:     g(x) = f(t,x,z+shift*x) = 0
    1037: */
    1038: static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts)
    1039: {
    1040:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    1041:   Vec      Z, Ydot;
    1042:   DM       dm, dmsave;

    1044:   PetscFunctionBegin;
    1045:   PetscCall(SNESGetDM(snes, &dm));
    1046:   PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
    1047:   PetscCall(VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z));
    1048:   dmsave = ts->dm;
    1049:   ts->dm = dm;
    1050:   PetscCall(TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE));
    1051:   ts->dm = dmsave;
    1052:   PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
    1053:   PetscFunctionReturn(PETSC_SUCCESS);
    1054: }

    1056: static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts)
    1057: {
    1058:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    1059:   Vec      Z, Ydot;
    1060:   DM       dm, dmsave;

    1062:   PetscFunctionBegin;
    1063:   PetscCall(SNESGetDM(snes, &dm));
    1064:   PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
    1065:   dmsave = ts->dm;
    1066:   ts->dm = dm;
    1067:   /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
    1068:   PetscCall(TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE));
    1069:   ts->dm = dmsave;
    1070:   PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
    1071:   PetscFunctionReturn(PETSC_SUCCESS);
    1072: }

    1074: static PetscErrorCode TSSetUp_GLLE(TS ts)
    1075: {
    1076:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    1077:   PetscInt max_r, max_s;
    1078:   DM       dm;

    1080:   PetscFunctionBegin;
    1081:   gl->setupcalled = PETSC_TRUE;
    1082:   PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
    1083:   PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->X));
    1084:   PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold));
    1085:   PetscCall(VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot));
    1086:   PetscCall(VecDuplicateVecs(ts->vec_sol, 3, &gl->himom));
    1087:   PetscCall(VecDuplicate(ts->vec_sol, &gl->W));
    1088:   PetscCall(VecDuplicate(ts->vec_sol, &gl->Y));
    1089:   PetscCall(VecDuplicate(ts->vec_sol, &gl->Z));

    1091:   /* Default acceptance tests and adaptivity */
    1092:   if (!gl->Accept) PetscCall(TSGLLESetAcceptType(ts, TSGLLEACCEPT_ALWAYS));
    1093:   if (!gl->adapt) PetscCall(TSGLLEGetAdapt(ts, &gl->adapt));

    1095:   if (gl->current_scheme < 0) {
    1096:     PetscInt i;
    1097:     for (i = 0;; i++) {
    1098:       if (gl->schemes[i]->p == gl->start_order) break;
    1099:       PetscCheck(i + 1 != gl->nschemes, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No schemes available with requested start order %" PetscInt_FMT, i);
    1100:     }
    1101:     gl->current_scheme = i;
    1102:   }
    1103:   PetscCall(TSGetDM(ts, &dm));
    1104:   PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
    1105:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
    1106:   PetscFunctionReturn(PETSC_SUCCESS);
    1107: }

    1109: static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems PetscOptionsObject)
    1110: {
    1111:   TS_GLLE *gl         = (TS_GLLE *)ts->data;
    1112:   char     tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify";

    1114:   PetscFunctionBegin;
    1115:   PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options");
    1116:   {
    1117:     PetscBool flg;
    1118:     PetscCall(PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg));
    1119:     if (flg || !gl->type_name[0]) PetscCall(TSGLLESetType(ts, tname));
    1120:     PetscCall(PetscOptionsInt("-ts_gl_max_step_rejections", "Maximum number of times to attempt a step", "None", gl->max_step_rejections, &gl->max_step_rejections, NULL));
    1121:     PetscCall(PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL));
    1122:     PetscCall(PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL));
    1123:     PetscCall(PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL));
    1124:     PetscCall(PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL));
    1125:     PetscCall(PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL));
    1126:     PetscCall(PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL));
    1127:     PetscCall(PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL));
    1128:     PetscCall(PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg));
    1129:     if (flg) {
    1130:       PetscBool match1, match2;
    1131:       PetscCall(PetscStrcmp(completef, "rescale", &match1));
    1132:       PetscCall(PetscStrcmp(completef, "rescale-and-modify", &match2));
    1133:       if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale;
    1134:       else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
    1135:       else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef);
    1136:     }
    1137:     {
    1138:       char type[256] = TSGLLEACCEPT_ALWAYS;
    1139:       PetscCall(PetscOptionsFList("-ts_gl_accept_type", "Method to use for determining whether to accept a step", "TSGLLESetAcceptType", TSGLLEAcceptList, gl->accept_name[0] ? gl->accept_name : type, type, sizeof(type), &flg));
    1140:       if (flg || !gl->accept_name[0]) PetscCall(TSGLLESetAcceptType(ts, type));
    1141:     }
    1142:     {
    1143:       TSGLLEAdapt adapt;
    1144:       PetscCall(TSGLLEGetAdapt(ts, &adapt));
    1145:       PetscCall(TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject));
    1146:     }
    1147:   }
    1148:   PetscOptionsHeadEnd();
    1149:   PetscFunctionReturn(PETSC_SUCCESS);
    1150: }

    1152: static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer)
    1153: {
    1154:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
    1155:   PetscInt  i;
    1156:   PetscBool isascii, details;

    1158:   PetscFunctionBegin;
    1159:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    1160:   if (isascii) {
    1161:     PetscCall(PetscViewerASCIIPrintf(viewer, "  min order %" PetscInt_FMT ", max order %" PetscInt_FMT ", current order %" PetscInt_FMT "\n", gl->min_order, gl->max_order, gl->schemes[gl->current_scheme]->p));
    1162:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction]));
    1163:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Extrapolation: %s\n", gl->extrapolate ? "yes" : "no"));
    1164:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)"));
    1165:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1166:     PetscCall(TSGLLEAdaptView(gl->adapt, viewer));
    1167:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1168:     PetscCall(PetscViewerASCIIPrintf(viewer, "  type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)"));
    1169:     PetscCall(PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes));
    1170:     details = PETSC_FALSE;
    1171:     PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL));
    1172:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1173:     for (i = 0; i < gl->nschemes; i++) PetscCall(TSGLLESchemeView(gl->schemes[i], details, viewer));
    1174:     if (gl->View) PetscCall((*gl->View)(gl, viewer));
    1175:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1176:   }
    1177:   PetscFunctionReturn(PETSC_SUCCESS);
    1178: }

    1180: /*@C
    1181:   TSGLLERegister -  adds a `TSGLLE` implementation

    1183:   Not Collective, No Fortran Support

    1185:   Input Parameters:
    1186: + sname    - name of user-defined general linear scheme
    1187: - function - routine to create method context

    1189:   Level: advanced

    1191:   Note:
    1192:   `TSGLLERegister()` may be called multiple times to add several user-defined families.

    1194:   Example Usage:
    1195: .vb
    1196:   TSGLLERegister("my_scheme", MySchemeCreate);
    1197: .ve

    1199:   Then, your scheme can be chosen with the procedural interface via
    1200: .vb
    1201:   TSGLLESetType(ts, "my_scheme")
    1202: .ve
    1203:   or at runtime via the option
    1204: .vb
    1205:   -ts_gl_type my_scheme
    1206: .ve

    1208: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLESetType()`
    1209: @*/
    1210: PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS))
    1211: {
    1212:   PetscFunctionBegin;
    1213:   PetscCall(TSGLLEInitializePackage());
    1214:   PetscCall(PetscFunctionListAdd(&TSGLLEList, sname, function));
    1215:   PetscFunctionReturn(PETSC_SUCCESS);
    1216: }

    1218: /*@C
    1219:   TSGLLEAcceptRegister -  adds a `TSGLLE` acceptance scheme

    1221:   Not Collective

    1223:   Input Parameters:
    1224: + sname    - name of user-defined acceptance scheme
    1225: - function - routine to create method context, see `TSGLLEAcceptFn` for the calling sequence

    1227:   Level: advanced

    1229:   Note:
    1230:   `TSGLLEAcceptRegister()` may be called multiple times to add several user-defined families.

    1232:   Example Usage:
    1233: .vb
    1234:   TSGLLEAcceptRegister("my_scheme", MySchemeCreate);
    1235: .ve

    1237:   Then, your scheme can be chosen with the procedural interface via
    1238: .vb
    1239:   TSGLLESetAcceptType(ts, "my_scheme")
    1240: .ve
    1241:   or at runtime via the option `-ts_gl_accept_type my_scheme`

    1243: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLEAcceptFn`
    1244: @*/
    1245: PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFn *function)
    1246: {
    1247:   PetscFunctionBegin;
    1248:   PetscCall(PetscFunctionListAdd(&TSGLLEAcceptList, sname, function));
    1249:   PetscFunctionReturn(PETSC_SUCCESS);
    1250: }

    1252: /*@C
    1253:   TSGLLERegisterAll - Registers all of the general linear methods in `TSGLLE`

    1255:   Not Collective

    1257:   Level: advanced

    1259: .seealso: [](ch_ts), `TSGLLE`, `TSGLLERegisterDestroy()`
    1260: @*/
    1261: PetscErrorCode TSGLLERegisterAll(void)
    1262: {
    1263:   PetscFunctionBegin;
    1264:   if (TSGLLERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
    1265:   TSGLLERegisterAllCalled = PETSC_TRUE;

    1267:   PetscCall(TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS));
    1268:   PetscCall(TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always));
    1269:   PetscFunctionReturn(PETSC_SUCCESS);
    1270: }

    1272: /*@C
    1273:   TSGLLEInitializePackage - This function initializes everything in the `TSGLLE` package. It is called
    1274:   from `TSInitializePackage()`.

    1276:   Level: developer

    1278: .seealso: [](ch_ts), `PetscInitialize()`, `TSInitializePackage()`, `TSGLLEFinalizePackage()`
    1279: @*/
    1280: PetscErrorCode TSGLLEInitializePackage(void)
    1281: {
    1282:   PetscFunctionBegin;
    1283:   if (TSGLLEPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
    1284:   TSGLLEPackageInitialized = PETSC_TRUE;
    1285:   PetscCall(TSGLLERegisterAll());
    1286:   PetscCall(PetscRegisterFinalize(TSGLLEFinalizePackage));
    1287:   PetscFunctionReturn(PETSC_SUCCESS);
    1288: }

    1290: /*@C
    1291:   TSGLLEFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
    1292:   called from `PetscFinalize()`.

    1294:   Level: developer

    1296: .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEInitializePackage()`, `TSInitializePackage()`
    1297: @*/
    1298: PetscErrorCode TSGLLEFinalizePackage(void)
    1299: {
    1300:   PetscFunctionBegin;
    1301:   PetscCall(PetscFunctionListDestroy(&TSGLLEList));
    1302:   PetscCall(PetscFunctionListDestroy(&TSGLLEAcceptList));
    1303:   TSGLLEPackageInitialized = PETSC_FALSE;
    1304:   TSGLLERegisterAllCalled  = PETSC_FALSE;
    1305:   PetscFunctionReturn(PETSC_SUCCESS);
    1306: }

    1308: /*MC
    1309:   TSGLLE - DAE solver using implicit General Linear methods {cite}`butcher_2007` {cite}`butcher2016numerical`

    1311:   Options Database Keys:
    1312: +  -ts_gl_type (irks)                               - the class of general linear method
    1313: .  -ts_gl_rtol tol                                  - relative error
    1314: .  -ts_gl_atol tol                                  - absolute error
    1315: .  -ts_gl_min_order p                               - minimum order method to consider (default=1)
    1316: .  -ts_gl_max_order p                               - maximum order method to consider (default=3)
    1317: .  -ts_gl_start_order p                             - order of starting method (default=1)
    1318: .  -ts_gl_complete (rescale|rescale-and-modify)     - method to use for completing the step (rescale-and-modify or rescale)
    1319: -  -ts_adapt_type (basic|dsp|none|cfl|glee|history) - adaptive controller to use (none step both)

    1321:   Level: beginner

    1323:   Notes:
    1324:   These methods contain Runge-Kutta and multistep schemes as special cases. These special cases
    1325:   have some fundamental limitations. For example, diagonally implicit Runge-Kutta cannot have
    1326:   stage order greater than 1 which limits their applicability to very stiff systems.
    1327:   Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF are not
    1328:   0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high
    1329:   stage order and reliable error estimates for both 1 and 2 orders higher to facilitate
    1330:   adaptive step sizes and adaptive order schemes. All this is possible while preserving a
    1331:   singly diagonally implicit structure.

    1333:   This integrator can be applied to DAE.

    1335:   Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit
    1336:   Runge-Kutta (DIRK). They are represented by the tableau

    1338: .vb
    1339:   A  |  U
    1340:   -------
    1341:   B  |  V
    1342: .ve

    1344:   combined with a vector c of abscissa. "Diagonally implicit" means that $A$ is lower
    1345:   triangular. A step of the general method reads

    1347:   $$
    1348:   \begin{align*}
    1349:   [ Y ] = [A  U] [  Y'   ] \\
    1350:   [X^k] = [B  V] [X^{k-1}]
    1351:   \end{align*}
    1352:   $$

    1354:   where Y is the multivector of stage values, $Y'$ is the multivector of stage derivatives, $X^k$
    1355:   is the Nordsieck vector of the solution at step $k$. The Nordsieck vector consists of the first
    1356:   $r$ moments of the solution, given by

    1358:   $$
    1359:   X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
    1360:   $$

    1362:   If $A$ is lower triangular, we can solve the stages $(Y, Y')$ sequentially

    1364:   $$
    1365:   y_i = h \sum_{j=0}^{s-1} (a_{ij} y'_j) + \sum_{j=0}^{r-1} u_{ij} x_j, \, \,   i=0,...,{s-1}
    1366:   $$

    1368:   and then construct the pieces to carry to the next step

    1370:   $$
    1371:   xx_i = h \sum_{j=0}^{s-1} b_{ij} y'_j  + \sum_{j=0}^{r-1} v_{ij} x_j,  \, \,  i=0,...,{r-1}
    1372:   $$

    1374:   Note that when the equations are cast in implicit form, we are using the stage equation to
    1375:   define $y'_i$ in terms of $y_i$ and known stuff ($y_j$ for $j<i$ and $x_j$ for all $j$).

    1377:   Error estimation

    1379:   At present, the most attractive GL methods for stiff problems are singly diagonally implicit
    1380:   schemes which posses Inherent Runge-Kutta Stability (`TSIRKS`).  These methods have $r=s$, the
    1381:   number of items passed between steps is equal to the number of stages.  The order and
    1382:   stage-order are one less than the number of stages.  We use the error estimates in the 2007
    1383:   paper which provide the following estimates

    1385:   $$
    1386:   \begin{align*}
    1387:   h^{p+1} X^{(p+1)}          = \phi_0^T Y' + [0 \psi_0^T] Xold \\
    1388:   h^{p+2} X^{(p+2)}          = \phi_1^T Y' + [0 \psi_1^T] Xold \\
    1389:   h^{p+2} (dx'/dx) X^{(p+1)} = \phi_2^T Y' + [0 \psi_2^T] Xold
    1390:   \end{align*}
    1391:   $$

    1393:   These estimates are accurate to $ O(h^{p+3})$.

    1395:   Changing the step size

    1397:   Uses the generalized "rescale and modify" scheme, see equation (4.5) of {cite}`butcher_2007`.

    1399: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType`, `TSGLLEType`, `TSGLLESetType()`, `TSGLLESetAcceptType()`, `TSGLLEGetAdapt()`
    1400: M*/
    1401: PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
    1402: {
    1403:   TS_GLLE *gl;

    1405:   PetscFunctionBegin;
    1406:   PetscCall(TSGLLEInitializePackage());

    1408:   PetscCall(PetscNew(&gl));
    1409:   ts->data = (void *)gl;

    1411:   ts->ops->reset          = TSReset_GLLE;
    1412:   ts->ops->destroy        = TSDestroy_GLLE;
    1413:   ts->ops->view           = TSView_GLLE;
    1414:   ts->ops->setup          = TSSetUp_GLLE;
    1415:   ts->ops->solve          = TSSolve_GLLE;
    1416:   ts->ops->setfromoptions = TSSetFromOptions_GLLE;
    1417:   ts->ops->snesfunction   = SNESTSFormFunction_GLLE;
    1418:   ts->ops->snesjacobian   = SNESTSFormJacobian_GLLE;

    1420:   ts->usessnes = PETSC_TRUE;

    1422:   gl->max_step_rejections = 1;
    1423:   gl->min_order           = 1;
    1424:   gl->max_order           = 3;
    1425:   gl->start_order         = 1;
    1426:   gl->current_scheme      = -1;
    1427:   gl->extrapolate         = PETSC_FALSE;

    1429:   gl->wrms_atol = 1e-8;
    1430:   gl->wrms_rtol = 1e-5;

    1432:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE));
    1433:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE));
    1434:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE));
    1435:   PetscFunctionReturn(PETSC_SUCCESS);
    1436: }