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) {
     39:       PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Z", Z));
     40:     } else *Z = gl->Z;
     41:   }
     42:   if (Ydotstage) {
     43:     if (dm && dm != ts->dm) {
     44:       PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage));
     45:     } else *Ydotstage = gl->Ydot[gl->stage];
     46:   }
     47:   PetscFunctionReturn(PETSC_SUCCESS);
     48: }

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

     62: static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, void *ctx)
     63: {
     64:   PetscFunctionBegin;
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
     69: {
     70:   TS  ts = (TS)ctx;
     71:   Vec Ydot, Ydot_c;

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

     83: static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, void *ctx)
     84: {
     85:   PetscFunctionBegin;
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
     90: {
     91:   TS  ts = (TS)ctx;
     92:   Vec Ydot, Ydot_s;

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

     98:   PetscCall(VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));
     99:   PetscCall(VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));

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

    106: 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)
    107: {
    108:   TSGLLEScheme scheme;
    109:   PetscInt     j;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    284:   *inscheme = scheme;
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

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

    297: static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
    298: {
    299:   PetscInt i;

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

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

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

    331: static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer)
    332: {
    333:   PetscBool iascii;

    335:   PetscFunctionBegin;
    336:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    337:   if (iascii) {
    338:     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));
    339:     PetscCall(PetscViewerASCIIPushTab(viewer));
    340:     PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s,  FSAL: %s\n", sc->stiffly_accurate ? "yes" : "no", sc->fsal ? "yes" : "no"));
    341:     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])));
    342:     PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->c, "Abscissas c"));
    343:     if (view_details) {
    344:       PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->s, sc->a, "A"));
    345:       PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->s, sc->b, "B"));
    346:       PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->r, sc->u, "U"));
    347:       PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->r, sc->v, "V"));

    349:       PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->s, sc->phi, "Error estimate phi"));
    350:       PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->r, sc->psi, "Error estimate psi"));
    351:       PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->alpha, "Modify alpha"));
    352:       PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->beta, "Modify beta"));
    353:       PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->gamma, "Modify gamma"));
    354:       PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->stage_error, "Stage error xi"));
    355:     }
    356:     PetscCall(PetscViewerASCIIPopTab(viewer));
    357:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported", ((PetscObject)viewer)->type_name);
    358:   PetscFunctionReturn(PETSC_SUCCESS);
    359: }

    361: static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[])
    362: {
    363:   PetscInt i;

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

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

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

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

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

    432: static PetscErrorCode TSGLLECreate_IRKS(TS ts)
    433: {
    434:   TS_GLLE *gl = (TS_GLLE *)ts->data;

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

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

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

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

    613:   Collective

    615:   Input Parameters:
    616: + ts   - the `TS` context
    617: - type - a method

    619:   Options Database Key:
    620: . -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks)

    622:   Level: intermediate

    624:   Notes:
    625:   See "petsc/include/petscts.h" for available methods (for instance)
    626: .    TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems)

    628:   Normally, it is best to use the `TSSetFromOptions()` command and then set the `TSGLLE` type
    629:   from the options database rather than by using this routine.  Using the options database
    630:   provides the user with maximum flexibility in evaluating the many different solvers.  The
    631:   `TSGLLESetType()` routine is provided for those situations where it is necessary to set the
    632:   timestepping solver independently of the command line or options database.  This might be the
    633:   case, for example, when the choice of solver changes during the execution of the program, and
    634:   the user's application is taking responsibility for choosing the appropriate method.

    636: .seealso: [](ch_ts), `TS`, `TSGLLEType`, `TSGLLE`
    637: @*/
    638: PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type)
    639: {
    640:   PetscFunctionBegin;
    642:   PetscAssertPointer(type, 2);
    643:   PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type));
    644:   PetscFunctionReturn(PETSC_SUCCESS);
    645: }

    647: /*@C
    648:   TSGLLESetAcceptType - sets the acceptance test for `TSGLLE`

    650:   Logically Collective

    652:   Input Parameters:
    653: + ts   - the `TS` context
    654: - type - the type

    656:   Options Database Key:
    657: . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step

    659:   Level: intermediate

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

    665: .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt`
    666: @*/
    667: PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type)
    668: {
    669:   PetscFunctionBegin;
    671:   PetscAssertPointer(type, 2);
    672:   PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type));
    673:   PetscFunctionReturn(PETSC_SUCCESS);
    674: }

    676: /*@
    677:   TSGLLEGetAdapt - gets the `TSGLLEAdapt` object from the `TS`

    679:   Not Collective

    681:   Input Parameter:
    682: . ts - the `TS` context

    684:   Output Parameter:
    685: . adapt - the `TSGLLEAdapt` context

    687:   Level: advanced

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

    693: .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegister()`
    694: @*/
    695: PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt)
    696: {
    697:   PetscFunctionBegin;
    699:   PetscAssertPointer(adapt, 2);
    700:   PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt));
    701:   PetscFunctionReturn(PETSC_SUCCESS);
    702: }

    704: static PetscErrorCode TSGLLEAccept_Always(TS ts, PetscReal tleft, PetscReal h, const PetscReal enorms[], PetscBool *accept)
    705: {
    706:   PetscFunctionBegin;
    707:   *accept = PETSC_TRUE;
    708:   PetscFunctionReturn(PETSC_SUCCESS);
    709: }

    711: static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
    712: {
    713:   TS_GLLE     *gl = (TS_GLLE *)ts->data;
    714:   PetscScalar *x, *w;
    715:   PetscInt     n, i;

    717:   PetscFunctionBegin;
    718:   PetscCall(VecGetArray(gl->X[0], &x));
    719:   PetscCall(VecGetArray(gl->W, &w));
    720:   PetscCall(VecGetLocalSize(gl->W, &n));
    721:   for (i = 0; i < n; i++) w[i] = 1. / (gl->wrms_atol + gl->wrms_rtol * PetscAbsScalar(x[i]));
    722:   PetscCall(VecRestoreArray(gl->X[0], &x));
    723:   PetscCall(VecRestoreArray(gl->W, &w));
    724:   PetscFunctionReturn(PETSC_SUCCESS);
    725: }

    727: static PetscErrorCode TSGLLEVecNormWRMS(TS ts, Vec X, PetscReal *nrm)
    728: {
    729:   TS_GLLE     *gl = (TS_GLLE *)ts->data;
    730:   PetscScalar *x, *w;
    731:   PetscReal    sum = 0.0, gsum;
    732:   PetscInt     n, N, i;

    734:   PetscFunctionBegin;
    735:   PetscCall(VecGetArray(X, &x));
    736:   PetscCall(VecGetArray(gl->W, &w));
    737:   PetscCall(VecGetLocalSize(gl->W, &n));
    738:   for (i = 0; i < n; i++) sum += PetscAbsScalar(PetscSqr(x[i] * w[i]));
    739:   PetscCall(VecRestoreArray(X, &x));
    740:   PetscCall(VecRestoreArray(gl->W, &w));
    741:   PetscCallMPI(MPIU_Allreduce(&sum, &gsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
    742:   PetscCall(VecGetSize(gl->W, &N));
    743:   *nrm = PetscSqrtReal(gsum / (1. * N));
    744:   PetscFunctionReturn(PETSC_SUCCESS);
    745: }

    747: static PetscErrorCode TSGLLESetType_GLLE(TS ts, TSGLLEType type)
    748: {
    749:   PetscBool same;
    750:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
    751:   PetscErrorCode (*r)(TS);

    753:   PetscFunctionBegin;
    754:   if (gl->type_name[0]) {
    755:     PetscCall(PetscStrcmp(gl->type_name, type, &same));
    756:     if (same) PetscFunctionReturn(PETSC_SUCCESS);
    757:     PetscCall((*gl->Destroy)(gl));
    758:   }

    760:   PetscCall(PetscFunctionListFind(TSGLLEList, type, &r));
    761:   PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLE type \"%s\" given", type);
    762:   PetscCall((*r)(ts));
    763:   PetscCall(PetscStrncpy(gl->type_name, type, sizeof(gl->type_name)));
    764:   PetscFunctionReturn(PETSC_SUCCESS);
    765: }

    767: static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type)
    768: {
    769:   TSGLLEAcceptFn *r;
    770:   TS_GLLE        *gl = (TS_GLLE *)ts->data;

    772:   PetscFunctionBegin;
    773:   PetscCall(PetscFunctionListFind(TSGLLEAcceptList, type, &r));
    774:   PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAccept type \"%s\" given", type);
    775:   gl->Accept = r;
    776:   PetscCall(PetscStrncpy(gl->accept_name, type, sizeof(gl->accept_name)));
    777:   PetscFunctionReturn(PETSC_SUCCESS);
    778: }

    780: static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt)
    781: {
    782:   TS_GLLE *gl = (TS_GLLE *)ts->data;

    784:   PetscFunctionBegin;
    785:   if (!gl->adapt) {
    786:     PetscCall(TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts), &gl->adapt));
    787:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)gl->adapt, (PetscObject)ts, 1));
    788:   }
    789:   *adapt = gl->adapt;
    790:   PetscFunctionReturn(PETSC_SUCCESS);
    791: }

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

    799:   PetscFunctionBegin;
    800:   cur   = -1;
    801:   cur_p = gl->schemes[gl->current_scheme]->p;
    802:   tleft = ts->max_time - (ts->ptime + ts->time_step);
    803:   for (i = 0, n = 0; i < gl->nschemes; i++) {
    804:     TSGLLEScheme sc = gl->schemes[i];
    805:     if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
    806:     if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[0];
    807:     else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[1];
    808:     else if (sc->p == cur_p + 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * (hmnorm[2] + hmnorm[3]);
    809:     else continue;
    810:     candidates[n] = i;
    811:     orders[n]     = PetscMin(sc->p, sc->q); /* order of global truncation error */
    812:     costs[n]      = sc->s;                  /* estimate the cost as the number of stages */
    813:     if (i == gl->current_scheme) cur = n;
    814:     n++;
    815:   }
    816:   PetscCheck(cur >= 0 && gl->nschemes > cur, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Current scheme not found in scheme list");
    817:   PetscCall(TSGLLEAdaptChoose(gl->adapt, n, orders, errors, costs, cur, h, tleft, &next_sc, next_h, finish));
    818:   *next_scheme = candidates[next_sc];
    819:   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,
    820:                       gl->schemes[*next_scheme]->r, gl->schemes[*next_scheme]->s, (double)*next_h, PetscBools[*finish]));
    821:   PetscFunctionReturn(PETSC_SUCCESS);
    822: }

    824: static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s)
    825: {
    826:   TS_GLLE *gl = (TS_GLLE *)ts->data;

    828:   PetscFunctionBegin;
    829:   *max_r = gl->schemes[gl->nschemes - 1]->r;
    830:   *max_s = gl->schemes[gl->nschemes - 1]->s;
    831:   PetscFunctionReturn(PETSC_SUCCESS);
    832: }

    834: static PetscErrorCode TSSolve_GLLE(TS ts)
    835: {
    836:   TS_GLLE            *gl = (TS_GLLE *)ts->data;
    837:   PetscInt            i, k, its, lits, max_r, max_s;
    838:   PetscBool           final_step, finish;
    839:   SNESConvergedReason snesreason;

    841:   PetscFunctionBegin;
    842:   PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

    844:   PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
    845:   PetscCall(VecCopy(ts->vec_sol, gl->X[0]));
    846:   for (i = 1; i < max_r; i++) PetscCall(VecZeroEntries(gl->X[i]));
    847:   PetscCall(TSGLLEUpdateWRMS(ts));

    849:   if (0) {
    850:     /* Find consistent initial data for DAE */
    851:     gl->stage_time = ts->ptime + ts->time_step;
    852:     gl->scoeff     = 1.;
    853:     gl->stage      = 0;

    855:     PetscCall(VecCopy(ts->vec_sol, gl->Z));
    856:     PetscCall(VecCopy(ts->vec_sol, gl->Y));
    857:     PetscCall(SNESSolve(ts->snes, NULL, gl->Y));
    858:     PetscCall(SNESGetIterationNumber(ts->snes, &its));
    859:     PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
    860:     PetscCall(SNESGetConvergedReason(ts->snes, &snesreason));

    862:     ts->snes_its += its;
    863:     ts->ksp_its += lits;
    864:     if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) {
    865:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
    866:       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));
    867:       PetscFunctionReturn(PETSC_SUCCESS);
    868:     }
    869:   }

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

    873:   for (k = 0, final_step = PETSC_FALSE, finish = PETSC_FALSE; k < ts->max_steps && !finish; k++) {
    874:     PetscInt           j, r, s, next_scheme = 0, rejections;
    875:     PetscReal          h, hmnorm[4], enorm[3], next_h;
    876:     PetscBool          accept;
    877:     const PetscScalar *c, *a, *u;
    878:     Vec               *X, *Ydot, Y;
    879:     TSGLLEScheme       scheme = gl->schemes[gl->current_scheme];

    881:     r    = scheme->r;
    882:     s    = scheme->s;
    883:     c    = scheme->c;
    884:     a    = scheme->a;
    885:     u    = scheme->u;
    886:     h    = ts->time_step;
    887:     X    = gl->X;
    888:     Ydot = gl->Ydot;
    889:     Y    = gl->Y;

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

    893:     /*
    894:       We only call PreStep at the start of each STEP, not each STAGE.  This is because it is
    895:       possible to fail (have to restart a step) after multiple stages.
    896:     */
    897:     PetscCall(TSPreStep(ts));

    899:     rejections = 0;
    900:     while (1) {
    901:       for (i = 0; i < s; i++) {
    902:         PetscScalar shift;
    903:         gl->scoeff     = 1. / PetscRealPart(a[i * s + i]);
    904:         shift          = gl->scoeff / ts->time_step;
    905:         gl->stage      = i;
    906:         gl->stage_time = ts->ptime + PetscRealPart(c[i]) * h;

    908:         /*
    909:         * Stage equation: Y = h A Y' + U X
    910:         * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
    911:         * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
    912:         * Then y'_i = z + 1/(h a_ii) y_i
    913:         */
    914:         PetscCall(VecZeroEntries(gl->Z));
    915:         for (j = 0; j < r; j++) PetscCall(VecAXPY(gl->Z, -shift * u[i * r + j], X[j]));
    916:         for (j = 0; j < i; j++) PetscCall(VecAXPY(gl->Z, -shift * h * a[i * s + j], Ydot[j]));
    917:         /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */

    919:         /* Compute an estimate of Y to start Newton iteration */
    920:         if (gl->extrapolate) {
    921:           if (i == 0) {
    922:             /* Linear extrapolation on the first stage */
    923:             PetscCall(VecWAXPY(Y, c[i] * h, X[1], X[0]));
    924:           } else {
    925:             /* Linear extrapolation from the last stage */
    926:             PetscCall(VecAXPY(Y, (c[i] - c[i - 1]) * h, Ydot[i - 1]));
    927:           }
    928:         } else if (i == 0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
    929:           PetscCall(VecCopy(X[0], Y));
    930:         }

    932:         /* Solve this stage (Ydot[i] is computed during function evaluation) */
    933:         PetscCall(SNESSolve(ts->snes, NULL, Y));
    934:         PetscCall(SNESGetIterationNumber(ts->snes, &its));
    935:         PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
    936:         PetscCall(SNESGetConvergedReason(ts->snes, &snesreason));
    937:         ts->snes_its += its;
    938:         ts->ksp_its += lits;
    939:         if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) {
    940:           ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
    941:           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));
    942:           PetscFunctionReturn(PETSC_SUCCESS);
    943:         }
    944:       }

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

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

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

    977:     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]));
    978:     if (!final_step) {
    979:       PetscCall(TSGLLEChooseNextScheme(ts, h, hmnorm, &next_scheme, &next_h, &final_step));
    980:     } else {
    981:       /* Dummy values to complete the current step in a consistent manner */
    982:       next_scheme = gl->current_scheme;
    983:       next_h      = h;
    984:       finish      = PETSC_TRUE;
    985:     }

    987:     X        = gl->Xold;
    988:     gl->Xold = gl->X;
    989:     gl->X    = X;
    990:     PetscCall((*gl->CompleteStep)(scheme, h, gl->schemes[next_scheme], next_h, Ydot, gl->Xold, gl->X));

    992:     PetscCall(TSGLLEUpdateWRMS(ts));

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

    999:     PetscCall(TSPostEvaluate(ts));
    1000:     PetscCall(TSPostStep(ts));
    1001:     PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

    1003:     gl->current_scheme = next_scheme;
    1004:     ts->time_step      = next_h;
    1005:   }
    1006:   PetscFunctionReturn(PETSC_SUCCESS);
    1007: }

    1009: /*------------------------------------------------------------*/

    1011: static PetscErrorCode TSReset_GLLE(TS ts)
    1012: {
    1013:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    1014:   PetscInt max_r, max_s;

    1016:   PetscFunctionBegin;
    1017:   if (gl->setupcalled) {
    1018:     PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
    1019:     PetscCall(VecDestroyVecs(max_r, &gl->Xold));
    1020:     PetscCall(VecDestroyVecs(max_r, &gl->X));
    1021:     PetscCall(VecDestroyVecs(max_s, &gl->Ydot));
    1022:     PetscCall(VecDestroyVecs(3, &gl->himom));
    1023:     PetscCall(VecDestroy(&gl->W));
    1024:     PetscCall(VecDestroy(&gl->Y));
    1025:     PetscCall(VecDestroy(&gl->Z));
    1026:   }
    1027:   gl->setupcalled = PETSC_FALSE;
    1028:   PetscFunctionReturn(PETSC_SUCCESS);
    1029: }

    1031: static PetscErrorCode TSDestroy_GLLE(TS ts)
    1032: {
    1033:   TS_GLLE *gl = (TS_GLLE *)ts->data;

    1035:   PetscFunctionBegin;
    1036:   PetscCall(TSReset_GLLE(ts));
    1037:   if (ts->dm) {
    1038:     PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
    1039:     PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
    1040:   }
    1041:   if (gl->adapt) PetscCall(TSGLLEAdaptDestroy(&gl->adapt));
    1042:   if (gl->Destroy) PetscCall((*gl->Destroy)(gl));
    1043:   PetscCall(PetscFree(ts->data));
    1044:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL));
    1045:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL));
    1046:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL));
    1047:   PetscFunctionReturn(PETSC_SUCCESS);
    1048: }

    1050: /*
    1051:     This defines the nonlinear equation that is to be solved with SNES
    1052:     g(x) = f(t,x,z+shift*x) = 0
    1053: */
    1054: static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts)
    1055: {
    1056:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    1057:   Vec      Z, Ydot;
    1058:   DM       dm, dmsave;

    1060:   PetscFunctionBegin;
    1061:   PetscCall(SNESGetDM(snes, &dm));
    1062:   PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
    1063:   PetscCall(VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z));
    1064:   dmsave = ts->dm;
    1065:   ts->dm = dm;
    1066:   PetscCall(TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE));
    1067:   ts->dm = dmsave;
    1068:   PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
    1069:   PetscFunctionReturn(PETSC_SUCCESS);
    1070: }

    1072: static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts)
    1073: {
    1074:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    1075:   Vec      Z, Ydot;
    1076:   DM       dm, dmsave;

    1078:   PetscFunctionBegin;
    1079:   PetscCall(SNESGetDM(snes, &dm));
    1080:   PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
    1081:   dmsave = ts->dm;
    1082:   ts->dm = dm;
    1083:   /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
    1084:   PetscCall(TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE));
    1085:   ts->dm = dmsave;
    1086:   PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
    1087:   PetscFunctionReturn(PETSC_SUCCESS);
    1088: }

    1090: static PetscErrorCode TSSetUp_GLLE(TS ts)
    1091: {
    1092:   TS_GLLE *gl = (TS_GLLE *)ts->data;
    1093:   PetscInt max_r, max_s;
    1094:   DM       dm;

    1096:   PetscFunctionBegin;
    1097:   gl->setupcalled = PETSC_TRUE;
    1098:   PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
    1099:   PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->X));
    1100:   PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold));
    1101:   PetscCall(VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot));
    1102:   PetscCall(VecDuplicateVecs(ts->vec_sol, 3, &gl->himom));
    1103:   PetscCall(VecDuplicate(ts->vec_sol, &gl->W));
    1104:   PetscCall(VecDuplicate(ts->vec_sol, &gl->Y));
    1105:   PetscCall(VecDuplicate(ts->vec_sol, &gl->Z));

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

    1111:   if (gl->current_scheme < 0) {
    1112:     PetscInt i;
    1113:     for (i = 0;; i++) {
    1114:       if (gl->schemes[i]->p == gl->start_order) break;
    1115:       PetscCheck(i + 1 != gl->nschemes, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No schemes available with requested start order %" PetscInt_FMT, i);
    1116:     }
    1117:     gl->current_scheme = i;
    1118:   }
    1119:   PetscCall(TSGetDM(ts, &dm));
    1120:   PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
    1121:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
    1122:   PetscFunctionReturn(PETSC_SUCCESS);
    1123: }
    1124: /*------------------------------------------------------------*/

    1126: static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems *PetscOptionsObject)
    1127: {
    1128:   TS_GLLE *gl         = (TS_GLLE *)ts->data;
    1129:   char     tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify";

    1131:   PetscFunctionBegin;
    1132:   PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options");
    1133:   {
    1134:     PetscBool flg;
    1135:     PetscCall(PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg));
    1136:     if (flg || !gl->type_name[0]) PetscCall(TSGLLESetType(ts, tname));
    1137:     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));
    1138:     PetscCall(PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL));
    1139:     PetscCall(PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL));
    1140:     PetscCall(PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL));
    1141:     PetscCall(PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL));
    1142:     PetscCall(PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL));
    1143:     PetscCall(PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL));
    1144:     PetscCall(PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL));
    1145:     PetscCall(PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg));
    1146:     if (flg) {
    1147:       PetscBool match1, match2;
    1148:       PetscCall(PetscStrcmp(completef, "rescale", &match1));
    1149:       PetscCall(PetscStrcmp(completef, "rescale-and-modify", &match2));
    1150:       if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale;
    1151:       else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
    1152:       else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef);
    1153:     }
    1154:     {
    1155:       char type[256] = TSGLLEACCEPT_ALWAYS;
    1156:       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));
    1157:       if (flg || !gl->accept_name[0]) PetscCall(TSGLLESetAcceptType(ts, type));
    1158:     }
    1159:     {
    1160:       TSGLLEAdapt adapt;
    1161:       PetscCall(TSGLLEGetAdapt(ts, &adapt));
    1162:       PetscCall(TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject));
    1163:     }
    1164:   }
    1165:   PetscOptionsHeadEnd();
    1166:   PetscFunctionReturn(PETSC_SUCCESS);
    1167: }

    1169: static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer)
    1170: {
    1171:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
    1172:   PetscInt  i;
    1173:   PetscBool iascii, details;

    1175:   PetscFunctionBegin;
    1176:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    1177:   if (iascii) {
    1178:     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));
    1179:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction]));
    1180:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Extrapolation: %s\n", gl->extrapolate ? "yes" : "no"));
    1181:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)"));
    1182:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1183:     PetscCall(TSGLLEAdaptView(gl->adapt, viewer));
    1184:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1185:     PetscCall(PetscViewerASCIIPrintf(viewer, "  type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)"));
    1186:     PetscCall(PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes));
    1187:     details = PETSC_FALSE;
    1188:     PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL));
    1189:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1190:     for (i = 0; i < gl->nschemes; i++) PetscCall(TSGLLESchemeView(gl->schemes[i], details, viewer));
    1191:     if (gl->View) PetscCall((*gl->View)(gl, viewer));
    1192:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1193:   }
    1194:   PetscFunctionReturn(PETSC_SUCCESS);
    1195: }

    1197: /*@C
    1198:   TSGLLERegister -  adds a `TSGLLE` implementation

    1200:   Not Collective, No Fortran Support

    1202:   Input Parameters:
    1203: + sname    - name of user-defined general linear scheme
    1204: - function - routine to create method context

    1206:   Level: advanced

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

    1211:   Example Usage:
    1212: .vb
    1213:   TSGLLERegister("my_scheme", MySchemeCreate);
    1214: .ve

    1216:   Then, your scheme can be chosen with the procedural interface via
    1217: $ TSGLLESetType(ts, "my_scheme")
    1218:   or at runtime via the option
    1219: $ -ts_gl_type my_scheme

    1221: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`
    1222: @*/
    1223: PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS))
    1224: {
    1225:   PetscFunctionBegin;
    1226:   PetscCall(TSGLLEInitializePackage());
    1227:   PetscCall(PetscFunctionListAdd(&TSGLLEList, sname, function));
    1228:   PetscFunctionReturn(PETSC_SUCCESS);
    1229: }

    1231: /*@C
    1232:   TSGLLEAcceptRegister -  adds a `TSGLLE` acceptance scheme

    1234:   Not Collective

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

    1240:   Level: advanced

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

    1245:   Example Usage:
    1246: .vb
    1247:   TSGLLEAcceptRegister("my_scheme", MySchemeCreate);
    1248: .ve

    1250:   Then, your scheme can be chosen with the procedural interface via
    1251: .vb
    1252:   TSGLLESetAcceptType(ts, "my_scheme")
    1253: .ve
    1254:   or at runtime via the option `-ts_gl_accept_type my_scheme`

    1256: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLEAcceptFn`
    1257: @*/
    1258: PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFn *function)
    1259: {
    1260:   PetscFunctionBegin;
    1261:   PetscCall(PetscFunctionListAdd(&TSGLLEAcceptList, sname, function));
    1262:   PetscFunctionReturn(PETSC_SUCCESS);
    1263: }

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

    1268:   Not Collective

    1270:   Level: advanced

    1272: .seealso: [](ch_ts), `TSGLLE`, `TSGLLERegisterDestroy()`
    1273: @*/
    1274: PetscErrorCode TSGLLERegisterAll(void)
    1275: {
    1276:   PetscFunctionBegin;
    1277:   if (TSGLLERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
    1278:   TSGLLERegisterAllCalled = PETSC_TRUE;

    1280:   PetscCall(TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS));
    1281:   PetscCall(TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always));
    1282:   PetscFunctionReturn(PETSC_SUCCESS);
    1283: }

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

    1289:   Level: developer

    1291: .seealso: [](ch_ts), `PetscInitialize()`, `TSInitializePackage()`, `TSGLLEFinalizePackage()`
    1292: @*/
    1293: PetscErrorCode TSGLLEInitializePackage(void)
    1294: {
    1295:   PetscFunctionBegin;
    1296:   if (TSGLLEPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
    1297:   TSGLLEPackageInitialized = PETSC_TRUE;
    1298:   PetscCall(TSGLLERegisterAll());
    1299:   PetscCall(PetscRegisterFinalize(TSGLLEFinalizePackage));
    1300:   PetscFunctionReturn(PETSC_SUCCESS);
    1301: }

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

    1307:   Level: developer

    1309: .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEInitializePackage()`, `TSInitializePackage()`
    1310: @*/
    1311: PetscErrorCode TSGLLEFinalizePackage(void)
    1312: {
    1313:   PetscFunctionBegin;
    1314:   PetscCall(PetscFunctionListDestroy(&TSGLLEList));
    1315:   PetscCall(PetscFunctionListDestroy(&TSGLLEAcceptList));
    1316:   TSGLLEPackageInitialized = PETSC_FALSE;
    1317:   TSGLLERegisterAllCalled  = PETSC_FALSE;
    1318:   PetscFunctionReturn(PETSC_SUCCESS);
    1319: }

    1321: /* ------------------------------------------------------------ */
    1322: /*MC
    1323:   TSGLLE - DAE solver using implicit General Linear methods {cite}`butcher_2007` {cite}`butcher2016numerical`

    1325:   Options Database Keys:
    1326: +  -ts_gl_type <type> - the class of general linear method (irks)
    1327: .  -ts_gl_rtol <tol>  - relative error
    1328: .  -ts_gl_atol <tol>  - absolute error
    1329: .  -ts_gl_min_order <p> - minimum order method to consider (default=1)
    1330: .  -ts_gl_max_order <p> - maximum order method to consider (default=3)
    1331: .  -ts_gl_start_order <p> - order of starting method (default=1)
    1332: .  -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
    1333: -  -ts_adapt_type <method> - adaptive controller to use (none step both)

    1335:   Level: beginner

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

    1347:   This integrator can be applied to DAE.

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

    1352: .vb
    1353:   A  |  U
    1354:   -------
    1355:   B  |  V
    1356: .ve

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

    1361:   $$
    1362:   \begin{align*}
    1363:   [ Y ] = [A  U] [  Y'   ] \\
    1364:   [X^k] = [B  V] [X^{k-1}]
    1365:   \end{align*}
    1366:   $$

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

    1372:   $$
    1373:   X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
    1374:   $$

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

    1378:   $$
    1379:   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}
    1380:   $$

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

    1384:   $$
    1385:   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}
    1386:   $$

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

    1391:   Error estimation

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

    1399:   $$
    1400:   \begin{align*}
    1401:   h^{p+1} X^{(p+1)}          = \phi_0^T Y' + [0 \psi_0^T] Xold \\
    1402:   h^{p+2} X^{(p+2)}          = \phi_1^T Y' + [0 \psi_1^T] Xold \\
    1403:   h^{p+2} (dx'/dx) X^{(p+1)} = \phi_2^T Y' + [0 \psi_2^T] Xold
    1404:   \end{align*}
    1405:   $$

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

    1409:   Changing the step size

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

    1413: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType`
    1414: M*/
    1415: PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
    1416: {
    1417:   TS_GLLE *gl;

    1419:   PetscFunctionBegin;
    1420:   PetscCall(TSGLLEInitializePackage());

    1422:   PetscCall(PetscNew(&gl));
    1423:   ts->data = (void *)gl;

    1425:   ts->ops->reset          = TSReset_GLLE;
    1426:   ts->ops->destroy        = TSDestroy_GLLE;
    1427:   ts->ops->view           = TSView_GLLE;
    1428:   ts->ops->setup          = TSSetUp_GLLE;
    1429:   ts->ops->solve          = TSSolve_GLLE;
    1430:   ts->ops->setfromoptions = TSSetFromOptions_GLLE;
    1431:   ts->ops->snesfunction   = SNESTSFormFunction_GLLE;
    1432:   ts->ops->snesjacobian   = SNESTSFormJacobian_GLLE;

    1434:   ts->usessnes = PETSC_TRUE;

    1436:   gl->max_step_rejections = 1;
    1437:   gl->min_order           = 1;
    1438:   gl->max_order           = 3;
    1439:   gl->start_order         = 1;
    1440:   gl->current_scheme      = -1;
    1441:   gl->extrapolate         = PETSC_FALSE;

    1443:   gl->wrms_atol = 1e-8;
    1444:   gl->wrms_rtol = 1e-5;

    1446:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE));
    1447:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE));
    1448:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE));
    1449:   PetscFunctionReturn(PETSC_SUCCESS);
    1450: }