/* Code for time stepping with the General Linear with Error Estimation method Notes: The general system is written as Udot = F(t,U) */ #include /*I "petscts.h" I*/ #include static PetscBool cited = PETSC_FALSE; static const char citation[] = "@ARTICLE{Constantinescu_TR2016b,\n" " author = {Constantinescu, E.M.},\n" " title = {Estimating Global Errors in Time Stepping},\n" " journal = {ArXiv e-prints},\n" " year = 2016,\n" " adsurl = {http://adsabs.harvard.edu/abs/2015arXiv150305166C}\n}\n"; static TSGLEEType TSGLEEDefaultType = TSGLEE35; static PetscBool TSGLEERegisterAllCalled; static PetscBool TSGLEEPackageInitialized; static PetscInt explicit_stage_time_id; typedef struct _GLEETableau *GLEETableau; struct _GLEETableau { char *name; PetscInt order; /* Classical approximation order of the method i*/ PetscInt s; /* Number of stages */ PetscInt r; /* Number of steps */ PetscReal gamma; /* LTE ratio */ PetscReal *A, *B, *U, *V, *S, *F, *c; /* Tableau */ PetscReal *Fembed; /* Embedded final method coefficients */ PetscReal *Ferror; /* Coefficients for computing error */ PetscReal *Serror; /* Coefficients for initializing the error */ PetscInt pinterp; /* Interpolation order */ PetscReal *binterp; /* Interpolation coefficients */ PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ }; typedef struct _GLEETableauLink *GLEETableauLink; struct _GLEETableauLink { struct _GLEETableau tab; GLEETableauLink next; }; static GLEETableauLink GLEETableauList; typedef struct { GLEETableau tableau; Vec *Y; /* Solution vector (along with auxiliary solution y~ or eps) */ Vec *X; /* Temporary solution vector */ Vec *YStage; /* Stage values */ Vec *YdotStage; /* Stage right-hand side */ Vec W; /* Right-hand-side for implicit stage solve */ Vec Ydot; /* Work vector holding Ydot during residual evaluation */ Vec yGErr; /* Vector holding the global error after a step is completed */ PetscScalar *swork; /* Scalar work (size of the number of stages)*/ PetscScalar *rwork; /* Scalar work (size of the number of steps)*/ PetscReal scoeff; /* shift = scoeff/dt */ PetscReal stage_time; TSStepStatus status; } TS_GLEE; /*MC TSGLEEi1 - Second order three stage implicit GLEE method This method has two stages. s = 3, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*MC TSGLEE23 - Second order three stage explicit GLEE method This method has three stages. s = 3, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*MC TSGLEE24 - Second order four stage explicit GLEE method This method has four stages. s = 4, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*MC TSGLEE25i - Second order five stage explicit GLEE method This method has five stages. s = 5, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*MC TSGLEE35 - Third order five stage explicit GLEE method This method has five stages. s = 5, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*MC TSGLEEEXRK2A - Second order six stage explicit GLEE method This method has six stages. s = 6, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*MC TSGLEERK32G1 - Third order eight stage explicit GLEE method This method has eight stages. s = 8, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*MC TSGLEERK285EX - Second order nine stage explicit GLEE method This method has nine stages. s = 9, r = 2 Level: advanced .seealso: [](ch_ts), `TSGLEE` M*/ /*@C TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in `TSGLEE` Not Collective, but should be called by all processes which will need the schemes to be registered Level: advanced .seealso: [](ch_ts), `TSGLEERegisterDestroy()` @*/ PetscErrorCode TSGLEERegisterAll(void) { PetscFunctionBegin; if (TSGLEERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); TSGLEERegisterAllCalled = PETSC_TRUE; { #define GAMMA 0.5 /* y-eps form */ const PetscInt p = 1, s = 3, r = 2; const PetscReal A[3][3] = { {1.0, 0, 0 }, {0, 0.5, 0 }, {0, 0.5, 0.5} }, B[2][3] = {{1.0, 0, 0}, {-2.0, 1.0, 1.0}}, U[3][2] = {{1.0, 0}, {1.0, 0.5}, {1.0, 0.5}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0}; PetscCall(TSGLEERegister(TSGLEEi1, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } { #undef GAMMA #define GAMMA 0.0 /* y-eps form */ const PetscInt p = 2, s = 3, r = 2; const PetscReal A[3][3] = { {0, 0, 0}, {1, 0, 0}, {0.25, 0.25, 0} }, B[2][3] = {{1.0 / 12.0, 1.0 / 12.0, 5.0 / 6.0}, {1.0 / 12.0, 1.0 / 12.0, -1.0 / 6.0}}, U[3][2] = {{1, 0}, {1, 10}, {1, -1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0}; PetscCall(TSGLEERegister(TSGLEE23, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } { #undef GAMMA #define GAMMA 0.0 /* y-y~ form */ const PetscInt p = 2, s = 4, r = 2; const PetscReal A[4][4] = { {0, 0, 0, 0}, {0.75, 0, 0, 0}, {0.25, 29.0 / 60.0, 0, 0}, {-21.0 / 44.0, 145.0 / 44.0, -20.0 / 11.0, 0} }, B[2][4] = {{109.0 / 275.0, 58.0 / 75.0, -37.0 / 110.0, 1.0 / 6.0}, {3.0 / 11.0, 0, 75.0 / 88.0, -1.0 / 8.0}}, U[4][2] = {{0, 1}, {75.0 / 58.0, -17.0 / 58.0}, {0, 1}, {0, 1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0}; PetscCall(TSGLEERegister(TSGLEE24, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } { #undef GAMMA #define GAMMA 0.0 /* y-y~ form */ const PetscInt p = 2, s = 5, r = 2; const PetscReal A[5][5] = { {0, 0, 0, 0, 0}, {-0.94079244066783383269, 0, 0, 0, 0}, {0.64228187778301907108, 0.10915356933958500042, 0, 0, 0}, {-0.51764297742287450812, 0.74414270351096040738, -0.71404164927824538121, 0, 0}, {-0.44696561556825969206, -0.76768425657590196518, 0.20111608138142987881, 0.93828186737840469796, 0} }, B[2][5] = {{-0.029309178948150356153, -0.49671981884013874923, 0.34275801517650053274, 0.32941112623949194988, 0.85385985637229662276}, {0.78133219686062535272, 0.074238691892675897635, 0.57957363498384957966, -0.24638502829674959968, -0.18875949544040123033}}, U[5][2] = {{0.16911424754448327735, 0.83088575245551672265}, {0.53638465733199574340, 0.46361534266800425660}, {0.39901579167169582526, 0.60098420832830417474}, {0.87689005530618575480, 0.12310994469381424520}, {0.99056100455550913009, 0.0094389954444908699092}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0}; PetscCall(TSGLEERegister(TSGLEE25I, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } { #undef GAMMA #define GAMMA 0.0 /* y-y~ form */ const PetscInt p = 3, s = 5, r = 2; const PetscReal A[5][5] = { {0, 0, 0, 0, 0}, {-2169604947363702313.0 / 24313474998937147335.0, 0, 0, 0, 0}, {46526746497697123895.0 / 94116917485856474137.0, -10297879244026594958.0 / 49199457603717988219.0, 0, 0, 0}, {23364788935845982499.0 / 87425311444725389446.0, -79205144337496116638.0 / 148994349441340815519.0, 40051189859317443782.0 / 36487615018004984309.0, 0, 0}, {42089522664062539205.0 / 124911313006412840286.0, -15074384760342762939.0 / 137927286865289746282.0, -62274678522253371016.0 / 125918573676298591413.0, 13755475729852471739.0 / 79257927066651693390.0, 0} }, B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0, -55810892792806293355.0 / 206957624151308356511.0, 24061048952676379087.0 / 158739347956038723465.0, 3577972206874351339.0 / 7599733370677197135.0, -59449832954780563947.0 / 137360038685338563670.0}, {-9738262186984159168.0 / 99299082461487742983.0, -32797097931948613195.0 / 61521565616362163366.0, 42895514606418420631.0 / 71714201188501437336.0, 22608567633166065068.0 / 55371917805607957003.0, 94655809487476459565.0 / 151517167160302729021.0}}, U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0, 10043614439674808267.0 / 80863923579509469826.0}, {161694774978034105510.0 / 106187653640211060371.0, -55507121337823045139.0 / 106187653640211060371.0}, {78486094644566264568.0 / 88171030896733822981.0, 9684936252167558413.0 / 88171030896733822981.0}, {65394922146334854435.0 / 84570853840405479554.0, 19175931694070625119.0 / 84570853840405479554.0}, {8607282770183754108.0 / 108658046436496925911.0, 100050763666313171803.0 / 108658046436496925911.0}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0}; PetscCall(TSGLEERegister(TSGLEE35, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } { #undef GAMMA #define GAMMA 0.25 /* y-eps form */ const PetscInt p = 2, s = 6, r = 2; const PetscReal A[6][6] = { {0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0.5, 0, 0, 0}, {0, 0, 0.25, 0.25, 0, 0}, {0, 0, 0.25, 0.25, 0.5, 0} }, B[2][6] = {{0.5, 0.5, 0, 0, 0, 0}, {-2.0 / 3.0, -2.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}}, U[6][2] = {{1, 0}, {1, 0}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0}; PetscCall(TSGLEERegister(TSGLEEEXRK2A, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } { #undef GAMMA #define GAMMA 0.0 /* y-eps form */ const PetscInt p = 3, s = 8, r = 2; const PetscReal A[8][8] = { {0, 0, 0, 0, 0, 0, 0, 0}, {0.5, 0, 0, 0, 0, 0, 0, 0}, {-1, 2, 0, 0, 0, 0, 0, 0}, {1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0}, {-7.0 / 24.0, 1.0 / 3.0, 1.0 / 12.0, -1.0 / 8.0, 0.5, 0, 0, 0}, {7.0 / 6.0, -4.0 / 3.0, -1.0 / 3.0, 0.5, -1.0, 2.0, 0, 0}, {0, 0, 0, 0, 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0} }, B[2][8] = {{1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0, 0, 0, 0, 0}, {-1.0 / 6.0, -2.0 / 3.0, -1.0 / 6.0, 0, 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0}}, U[8][2] = {{1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 1}, {1, 1}, {1, 1}, {1, 1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0}; PetscCall(TSGLEERegister(TSGLEERK32G1, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } { #undef GAMMA #define GAMMA 0.25 /* y-eps form */ const PetscInt p = 2, s = 9, r = 2; const PetscReal A[9][9] = { {0, 0, 0, 0, 0, 0, 0, 0, 0}, {0.585786437626904966, 0, 0, 0, 0, 0, 0, 0, 0}, {0.149999999999999994, 0.849999999999999978, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0.292893218813452483, 0, 0, 0, 0, 0}, {0, 0, 0, 0.0749999999999999972, 0.424999999999999989, 0, 0, 0, 0}, {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0, 0, 0}, {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0.292893218813452483, 0, 0}, {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0.0749999999999999972, 0.424999999999999989, 0} }, B[2][9] = {{0.353553390593273786, 0.353553390593273786, 0.292893218813452483, 0, 0, 0, 0, 0, 0}, {-0.471404520791031678, -0.471404520791031678, -0.390524291751269959, 0.235702260395515839, 0.235702260395515839, 0.195262145875634979, 0.235702260395515839, 0.235702260395515839, 0.195262145875634979}}, U[9][2] = {{1, 0}, {1, 0}, {1, 0}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0}; PetscCall(TSGLEERegister(TSGLEERK285EX, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL)); } PetscFunctionReturn(PETSC_SUCCESS); } /*@C TSGLEERegisterDestroy - Frees the list of schemes that were registered by `TSGLEERegister()`. Not Collective Level: advanced .seealso: [](ch_ts), `TSGLEERegister()`, `TSGLEERegisterAll()` @*/ PetscErrorCode TSGLEERegisterDestroy(void) { GLEETableauLink link; PetscFunctionBegin; while ((link = GLEETableauList)) { GLEETableau t = &link->tab; GLEETableauList = link->next; PetscCall(PetscFree5(t->A, t->B, t->U, t->V, t->c)); PetscCall(PetscFree2(t->S, t->F)); PetscCall(PetscFree(t->Fembed)); PetscCall(PetscFree(t->Ferror)); PetscCall(PetscFree(t->Serror)); PetscCall(PetscFree(t->binterp)); PetscCall(PetscFree(t->name)); PetscCall(PetscFree(link)); } TSGLEERegisterAllCalled = PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } /*@C TSGLEEInitializePackage - This function initializes everything in the `TSGLEE` package. It is called from `TSInitializePackage()`. Level: developer .seealso: [](ch_ts), `PetscInitialize()` @*/ PetscErrorCode TSGLEEInitializePackage(void) { PetscFunctionBegin; if (TSGLEEPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); TSGLEEPackageInitialized = PETSC_TRUE; PetscCall(TSGLEERegisterAll()); PetscCall(PetscObjectComposedDataRegister(&explicit_stage_time_id)); PetscCall(PetscRegisterFinalize(TSGLEEFinalizePackage)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C TSGLEEFinalizePackage - This function destroys everything in the `TSGLEE` package. It is called from `PetscFinalize()`. Level: developer .seealso: [](ch_ts), `PetscFinalize()` @*/ PetscErrorCode TSGLEEFinalizePackage(void) { PetscFunctionBegin; TSGLEEPackageInitialized = PETSC_FALSE; PetscCall(TSGLEERegisterDestroy()); PetscFunctionReturn(PETSC_SUCCESS); } /*@C TSGLEERegister - register a new `TSGLEE` scheme by providing the entries in the Butcher tableau Not Collective, but the same schemes should be registered on all processes on which they will be used, No Fortran Support Input Parameters: + name - identifier for method . order - order of method . s - number of stages . r - number of steps . gamma - LTE ratio . A - stage coefficients (dimension s*s, row-major) . B - step completion coefficients (dimension r*s, row-major) . U - method coefficients (dimension s*r, row-major) . V - method coefficients (dimension r*r, row-major) . S - starting coefficients . F - finishing coefficients . c - abscissa (dimension s; NULL to use row sums of A) . Fembed - step completion coefficients for embedded method . Ferror - error computation coefficients . Serror - error initialization coefficients . pinterp - order of interpolation (0 if unavailable) - binterp - array of interpolation coefficients (NULL if unavailable) Level: advanced Note: Several `TSGLEE` methods are provided, this function is only needed to create new methods. .seealso: [](ch_ts), `TSGLEE` @*/ PetscErrorCode TSGLEERegister(TSGLEEType name, PetscInt order, PetscInt s, PetscInt r, PetscReal gamma, const PetscReal A[], const PetscReal B[], const PetscReal U[], const PetscReal V[], const PetscReal S[], const PetscReal F[], const PetscReal c[], const PetscReal Fembed[], const PetscReal Ferror[], const PetscReal Serror[], PetscInt pinterp, const PetscReal binterp[]) { GLEETableauLink link; GLEETableau t; PetscInt i, j; PetscFunctionBegin; PetscCall(TSGLEEInitializePackage()); PetscCall(PetscNew(&link)); t = &link->tab; PetscCall(PetscStrallocpy(name, &t->name)); t->order = order; t->s = s; t->r = r; t->gamma = gamma; PetscCall(PetscMalloc5(s * s, &t->A, r * r, &t->V, s, &t->c, r * s, &t->B, s * r, &t->U)); PetscCall(PetscMalloc2(r, &t->S, r, &t->F)); PetscCall(PetscArraycpy(t->A, A, s * s)); PetscCall(PetscArraycpy(t->B, B, r * s)); PetscCall(PetscArraycpy(t->U, U, s * r)); PetscCall(PetscArraycpy(t->V, V, r * r)); PetscCall(PetscArraycpy(t->S, S, r)); PetscCall(PetscArraycpy(t->F, F, r)); if (c) { PetscCall(PetscArraycpy(t->c, c, s)); } else { for (i = 0; i < s; i++) for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j]; } PetscCall(PetscMalloc1(r, &t->Fembed)); PetscCall(PetscMalloc1(r, &t->Ferror)); PetscCall(PetscMalloc1(r, &t->Serror)); PetscCall(PetscArraycpy(t->Fembed, Fembed, r)); PetscCall(PetscArraycpy(t->Ferror, Ferror, r)); PetscCall(PetscArraycpy(t->Serror, Serror, r)); t->pinterp = pinterp; PetscCall(PetscMalloc1(s * pinterp, &t->binterp)); PetscCall(PetscArraycpy(t->binterp, binterp, s * pinterp)); link->next = GLEETableauList; GLEETableauList = link; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSEvaluateStep_GLEE(TS ts, PetscInt order, Vec X, PetscBool *done) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; PetscReal h, *B = tab->B, *V = tab->V, *F = tab->F, *Fembed = tab->Fembed; PetscInt s = tab->s, r = tab->r, i, j; Vec *Y = glee->Y, *YdotStage = glee->YdotStage; PetscScalar *ws = glee->swork, *wr = glee->rwork; PetscFunctionBegin; switch (glee->status) { case TS_STEP_INCOMPLETE: case TS_STEP_PENDING: h = ts->time_step; break; case TS_STEP_COMPLETE: h = ts->ptime - ts->ptime_prev; break; default: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); } if (order == tab->order) { /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE or TS_STEP_COMPLETE, glee->X has the solution at the beginning of the time step. So no need to roll-back. */ if (glee->status == TS_STEP_INCOMPLETE) { for (i = 0; i < r; i++) { PetscCall(VecZeroEntries(Y[i])); for (j = 0; j < r; j++) wr[j] = V[i * r + j]; PetscCall(VecMAXPY(Y[i], r, wr, glee->X)); for (j = 0; j < s; j++) ws[j] = h * B[i * s + j]; PetscCall(VecMAXPY(Y[i], s, ws, YdotStage)); } PetscCall(VecZeroEntries(X)); for (j = 0; j < r; j++) wr[j] = F[j]; PetscCall(VecMAXPY(X, r, wr, Y)); } else PetscCall(VecCopy(ts->vec_sol, X)); PetscFunctionReturn(PETSC_SUCCESS); } else if (order == tab->order - 1) { /* Complete with the embedded method (Fembed) */ for (i = 0; i < r; i++) { PetscCall(VecZeroEntries(Y[i])); for (j = 0; j < r; j++) wr[j] = V[i * r + j]; PetscCall(VecMAXPY(Y[i], r, wr, glee->X)); for (j = 0; j < s; j++) ws[j] = h * B[i * s + j]; PetscCall(VecMAXPY(Y[i], s, ws, YdotStage)); } PetscCall(VecZeroEntries(X)); for (j = 0; j < r; j++) wr[j] = Fembed[j]; PetscCall(VecMAXPY(X, r, wr, Y)); if (done) *done = PETSC_TRUE; PetscFunctionReturn(PETSC_SUCCESS); } PetscCheck(done, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "GLEE '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT, tab->name, tab->order, order); *done = PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSStep_GLEE(TS ts) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; const PetscInt s = tab->s, r = tab->r; PetscReal *A = tab->A, *U = tab->U, *F = tab->F, *c = tab->c; Vec *Y = glee->Y, *X = glee->X, *YStage = glee->YStage, *YdotStage = glee->YdotStage, W = glee->W; SNES snes; PetscScalar *ws = glee->swork, *wr = glee->rwork; TSAdapt adapt; PetscInt i, j, reject, next_scheme, its, lits; PetscReal next_time_step; PetscReal t; PetscBool accept; PetscFunctionBegin; PetscCall(PetscCitationsRegister(citation, &cited)); for (i = 0; i < r; i++) PetscCall(VecCopy(Y[i], X[i])); PetscCall(TSGetSNES(ts, &snes)); next_time_step = ts->time_step; t = ts->ptime; accept = PETSC_TRUE; glee->status = TS_STEP_INCOMPLETE; for (reject = 0; reject < ts->max_reject && !ts->reason; reject++, ts->reject++) { PetscReal h = ts->time_step; PetscCall(TSPreStep(ts)); for (i = 0; i < s; i++) { glee->stage_time = t + h * c[i]; PetscCall(TSPreStage(ts, glee->stage_time)); if (A[i * s + i] == 0) { /* Explicit stage */ PetscCall(VecZeroEntries(YStage[i])); for (j = 0; j < r; j++) wr[j] = U[i * r + j]; PetscCall(VecMAXPY(YStage[i], r, wr, X)); for (j = 0; j < i; j++) ws[j] = h * A[i * s + j]; PetscCall(VecMAXPY(YStage[i], i, ws, YdotStage)); } else { /* Implicit stage */ glee->scoeff = 1.0 / A[i * s + i]; /* compute right-hand side */ PetscCall(VecZeroEntries(W)); for (j = 0; j < r; j++) wr[j] = U[i * r + j]; PetscCall(VecMAXPY(W, r, wr, X)); for (j = 0; j < i; j++) ws[j] = h * A[i * s + j]; PetscCall(VecMAXPY(W, i, ws, YdotStage)); PetscCall(VecScale(W, glee->scoeff / h)); /* set initial guess */ PetscCall(VecCopy(i > 0 ? YStage[i - 1] : ts->vec_sol, YStage[i])); /* solve for this stage */ PetscCall(SNESSolve(snes, W, YStage[i])); PetscCall(SNESGetIterationNumber(snes, &its)); PetscCall(SNESGetLinearSolveIterations(snes, &lits)); ts->snes_its += its; ts->ksp_its += lits; } PetscCall(TSGetAdapt(ts, &adapt)); PetscCall(TSAdaptCheckStage(adapt, ts, glee->stage_time, YStage[i], &accept)); if (!accept) goto reject_step; PetscCall(TSPostStage(ts, glee->stage_time, i, YStage)); PetscCall(TSComputeRHSFunction(ts, t + h * c[i], YStage[i], YdotStage[i])); } PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL)); glee->status = TS_STEP_PENDING; /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ PetscCall(TSGetAdapt(ts, &adapt)); PetscCall(TSAdaptCandidatesClear(adapt)); PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, &next_scheme, &next_time_step, &accept)); if (accept) { /* ignore next_scheme for now */ ts->ptime += ts->time_step; ts->time_step = next_time_step; glee->status = TS_STEP_COMPLETE; /* compute and store the global error */ /* Note: this is not needed if TSAdaptGLEE is not used */ PetscCall(TSGetTimeError(ts, 0, &glee->yGErr)); PetscCall(PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol, explicit_stage_time_id, ts->ptime)); break; } else { /* Roll back the current step */ for (j = 0; j < r; j++) wr[j] = F[j]; PetscCall(VecMAXPY(ts->vec_sol, r, wr, X)); ts->time_step = next_time_step; glee->status = TS_STEP_INCOMPLETE; } reject_step: continue; } if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSInterpolate_GLEE(TS ts, PetscReal itime, Vec X) { TS_GLEE *glee = (TS_GLEE *)ts->data; PetscInt s = glee->tableau->s, pinterp = glee->tableau->pinterp, i, j; PetscReal h, tt, t; PetscScalar *b; const PetscReal *B = glee->tableau->binterp; PetscFunctionBegin; PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSGLEE %s does not have an interpolation formula", glee->tableau->name); switch (glee->status) { case TS_STEP_INCOMPLETE: case TS_STEP_PENDING: h = ts->time_step; t = (itime - ts->ptime) / h; break; case TS_STEP_COMPLETE: h = ts->ptime - ts->ptime_prev; t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ break; default: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); } PetscCall(PetscMalloc1(s, &b)); for (i = 0; i < s; i++) b[i] = 0; for (j = 0, tt = t; j < pinterp; j++, tt *= t) { for (i = 0; i < s; i++) b[i] += h * B[i * pinterp + j] * tt; } PetscCall(VecCopy(glee->YStage[0], X)); PetscCall(VecMAXPY(X, s, b, glee->YdotStage)); PetscCall(PetscFree(b)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSReset_GLEE(TS ts) { TS_GLEE *glee = (TS_GLEE *)ts->data; PetscInt s, r; PetscFunctionBegin; if (!glee->tableau) PetscFunctionReturn(PETSC_SUCCESS); s = glee->tableau->s; r = glee->tableau->r; PetscCall(VecDestroyVecs(r, &glee->Y)); PetscCall(VecDestroyVecs(r, &glee->X)); PetscCall(VecDestroyVecs(s, &glee->YStage)); PetscCall(VecDestroyVecs(s, &glee->YdotStage)); PetscCall(VecDestroy(&glee->Ydot)); PetscCall(VecDestroy(&glee->yGErr)); PetscCall(VecDestroy(&glee->W)); PetscCall(PetscFree2(glee->swork, glee->rwork)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSGLEEGetVecs(TS ts, DM dm, Vec *Ydot) { TS_GLEE *glee = (TS_GLEE *)ts->data; PetscFunctionBegin; if (Ydot) { if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot)); else *Ydot = glee->Ydot; } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSGLEERestoreVecs(TS ts, DM dm, Vec *Ydot) { PetscFunctionBegin; if (Ydot) { if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot)); } PetscFunctionReturn(PETSC_SUCCESS); } /* This defines the nonlinear equation that is to be solved with SNES */ static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes, Vec X, Vec F, TS ts) { TS_GLEE *glee = (TS_GLEE *)ts->data; DM dm, dmsave; Vec Ydot; PetscReal shift = glee->scoeff / ts->time_step; PetscFunctionBegin; PetscCall(SNESGetDM(snes, &dm)); PetscCall(TSGLEEGetVecs(ts, dm, &Ydot)); /* Set Ydot = shift*X */ PetscCall(VecCopy(X, Ydot)); PetscCall(VecScale(Ydot, shift)); dmsave = ts->dm; ts->dm = dm; PetscCall(TSComputeIFunction(ts, glee->stage_time, X, Ydot, F, PETSC_FALSE)); ts->dm = dmsave; PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes, Vec X, Mat A, Mat B, TS ts) { TS_GLEE *glee = (TS_GLEE *)ts->data; DM dm, dmsave; Vec Ydot; PetscReal shift = glee->scoeff / ts->time_step; PetscFunctionBegin; PetscCall(SNESGetDM(snes, &dm)); PetscCall(TSGLEEGetVecs(ts, dm, &Ydot)); /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ dmsave = ts->dm; ts->dm = dm; PetscCall(TSComputeIJacobian(ts, glee->stage_time, X, Ydot, shift, A, B, PETSC_FALSE)); ts->dm = dmsave; PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine, DM coarse, PetscCtx ctx) { PetscFunctionBegin; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMRestrictHook_TSGLEE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx) { PetscFunctionBegin; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm, DM subdm, PetscCtx ctx) { PetscFunctionBegin; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx) { PetscFunctionBegin; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSSetUp_GLEE(TS ts) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab; PetscInt s, r; DM dm; PetscFunctionBegin; if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType)); tab = glee->tableau; s = tab->s; r = tab->r; PetscCall(VecDuplicateVecs(ts->vec_sol, r, &glee->Y)); PetscCall(VecDuplicateVecs(ts->vec_sol, r, &glee->X)); PetscCall(VecDuplicateVecs(ts->vec_sol, s, &glee->YStage)); PetscCall(VecDuplicateVecs(ts->vec_sol, s, &glee->YdotStage)); PetscCall(VecDuplicate(ts->vec_sol, &glee->Ydot)); PetscCall(VecDuplicate(ts->vec_sol, &glee->yGErr)); PetscCall(VecZeroEntries(glee->yGErr)); PetscCall(VecDuplicate(ts->vec_sol, &glee->W)); PetscCall(PetscMalloc2(s, &glee->swork, r, &glee->rwork)); PetscCall(TSGetDM(ts, &dm)); PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts)); PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSStartingMethod_GLEE(TS ts) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; PetscInt r = tab->r, i; PetscReal *S = tab->S; PetscFunctionBegin; for (i = 0; i < r; i++) { PetscCall(VecZeroEntries(glee->Y[i])); PetscCall(VecAXPY(glee->Y[i], S[i], ts->vec_sol)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSSetFromOptions_GLEE(TS ts, PetscOptionItems PetscOptionsObject) { char gleetype[256]; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "GLEE ODE solver options"); { GLEETableauLink link; PetscInt count, choice; PetscBool flg; const char **namelist; PetscCall(PetscStrncpy(gleetype, TSGLEEDefaultType, sizeof(gleetype))); for (link = GLEETableauList, count = 0; link; link = link->next, count++); PetscCall(PetscMalloc1(count, (char ***)&namelist)); for (link = GLEETableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; PetscCall(PetscOptionsEList("-ts_glee_type", "Family of GLEE method", "TSGLEESetType", (const char *const *)namelist, count, gleetype, &choice, &flg)); PetscCall(TSGLEESetType(ts, flg ? namelist[choice] : gleetype)); PetscCall(PetscFree(namelist)); } PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSView_GLEE(TS ts, PetscViewer viewer) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; PetscBool isascii; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { TSGLEEType gleetype; char buf[512]; PetscCall(TSGLEEGetType(ts, &gleetype)); PetscCall(PetscViewerASCIIPrintf(viewer, " GLEE type %s\n", gleetype)); PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf)); /* Note: print out r as well */ } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSLoad_GLEE(TS ts, PetscViewer viewer) { SNES snes; TSAdapt tsadapt; PetscFunctionBegin; PetscCall(TSGetAdapt(ts, &tsadapt)); PetscCall(TSAdaptLoad(tsadapt, viewer)); PetscCall(TSGetSNES(ts, &snes)); PetscCall(SNESLoad(snes, viewer)); /* function and Jacobian context for SNES when used with TS is always ts object */ PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ TSGLEESetType - Set the type of `TSGLEE` scheme Logically Collective Input Parameters: + ts - timestepping context - gleetype - type of `TSGLEE` scheme Level: intermediate .seealso: [](ch_ts), `TSGLEEGetType()`, `TSGLEE` @*/ PetscErrorCode TSGLEESetType(TS ts, TSGLEEType gleetype) { PetscFunctionBegin; PetscValidHeaderSpecific(ts, TS_CLASSID, 1); PetscAssertPointer(gleetype, 2); PetscTryMethod(ts, "TSGLEESetType_C", (TS, TSGLEEType), (ts, gleetype)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ TSGLEEGetType - Get the type of `TSGLEE` scheme Logically Collective Input Parameter: . ts - timestepping context Output Parameter: . gleetype - type of `TSGLEE` scheme Level: intermediate .seealso: [](ch_ts), `TSGLEE`, `TSGLEESetType()` @*/ PetscErrorCode TSGLEEGetType(TS ts, TSGLEEType *gleetype) { PetscFunctionBegin; PetscValidHeaderSpecific(ts, TS_CLASSID, 1); PetscUseMethod(ts, "TSGLEEGetType_C", (TS, TSGLEEType *), (ts, gleetype)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSGLEEGetType_GLEE(TS ts, TSGLEEType *gleetype) { TS_GLEE *glee = (TS_GLEE *)ts->data; PetscFunctionBegin; if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType)); *gleetype = glee->tableau->name; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSGLEESetType_GLEE(TS ts, TSGLEEType gleetype) { TS_GLEE *glee = (TS_GLEE *)ts->data; PetscBool match; GLEETableauLink link; PetscFunctionBegin; if (glee->tableau) { PetscCall(PetscStrcmp(glee->tableau->name, gleetype, &match)); if (match) PetscFunctionReturn(PETSC_SUCCESS); } for (link = GLEETableauList; link; link = link->next) { PetscCall(PetscStrcmp(link->tab.name, gleetype, &match)); if (match) { PetscCall(TSReset_GLEE(ts)); glee->tableau = &link->tab; PetscFunctionReturn(PETSC_SUCCESS); } } SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", gleetype); } static PetscErrorCode TSGetStages_GLEE(TS ts, PetscInt *ns, Vec **Y) { TS_GLEE *glee = (TS_GLEE *)ts->data; PetscFunctionBegin; if (ns) *ns = glee->tableau->s; if (Y) *Y = glee->YStage; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSGetSolutionComponents_GLEE(TS ts, PetscInt *n, Vec *Y) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; PetscFunctionBegin; if (!Y) *n = tab->r; else { PetscCheck(*n >= 0 && *n < tab->r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Second argument (%" PetscInt_FMT ") out of range[0,%" PetscInt_FMT "].", *n, tab->r - 1); PetscCall(VecCopy(glee->Y[*n], *Y)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSGetAuxSolution_GLEE(TS ts, Vec *X) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; PetscReal *F = tab->Fembed; PetscInt r = tab->r; Vec *Y = glee->Y; PetscScalar *wr = glee->rwork; PetscInt i; PetscFunctionBegin; PetscCall(VecZeroEntries(*X)); for (i = 0; i < r; i++) wr[i] = F[i]; PetscCall(VecMAXPY(*X, r, wr, Y)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSGetTimeError_GLEE(TS ts, PetscInt n, Vec *X) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; PetscReal *F = tab->Ferror; PetscInt r = tab->r; Vec *Y = glee->Y; PetscScalar *wr = glee->rwork; PetscInt i; PetscFunctionBegin; PetscCall(VecZeroEntries(*X)); if (n == 0) { for (i = 0; i < r; i++) wr[i] = F[i]; PetscCall(VecMAXPY(*X, r, wr, Y)); } else if (n == -1) { *X = glee->yGErr; } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSSetTimeError_GLEE(TS ts, Vec X) { TS_GLEE *glee = (TS_GLEE *)ts->data; GLEETableau tab = glee->tableau; PetscReal *S = tab->Serror; PetscInt r = tab->r, i; Vec *Y = glee->Y; PetscFunctionBegin; PetscCheck(r == 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSSetTimeError_GLEE not supported for '%s' with r=%" PetscInt_FMT ".", tab->name, tab->r); for (i = 1; i < r; i++) { PetscCall(VecCopy(ts->vec_sol, Y[i])); PetscCall(VecAXPBY(Y[i], S[0], S[1], X)); PetscCall(VecCopy(X, glee->yGErr)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode TSDestroy_GLEE(TS ts) { PetscFunctionBegin; PetscCall(TSReset_GLEE(ts)); if (ts->dm) { PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts)); PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts)); } PetscCall(PetscFree(ts->data)); PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", NULL)); PetscFunctionReturn(PETSC_SUCCESS); } /*MC TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes The user should provide the right-hand side of the equation using `TSSetRHSFunction()` and for `TSGLEEi1` the Jacobian of the right-hand side using `TSSetRHSJacobian()` Level: beginner Notes: The default is `TSGLEE35`, it can be changed with `TSGLEESetType()` or `-ts_glee_type type` The only implicit scheme is `TSGLEEi1` .seealso: [](ch_ts), [](sec_ts_glee), `TSCreate()`, `TS`, `TSSetType()`, `TSGLEESetType()`, `TSGLEEGetType()`, `TSGLEE23`, `TSGLEE24`, `TSGLEE35`, `TSGLEE25I`, `TSGLEEEXRK2A`, `TSGLEERK32G1`, `TSGLEERK285EX`, `TSGLEEType`, `TSGLEERegister()`, `TSType` M*/ PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) { TS_GLEE *th; PetscFunctionBegin; PetscCall(TSGLEEInitializePackage()); ts->ops->reset = TSReset_GLEE; ts->ops->destroy = TSDestroy_GLEE; ts->ops->view = TSView_GLEE; ts->ops->load = TSLoad_GLEE; ts->ops->setup = TSSetUp_GLEE; ts->ops->step = TSStep_GLEE; ts->ops->interpolate = TSInterpolate_GLEE; ts->ops->evaluatestep = TSEvaluateStep_GLEE; ts->ops->setfromoptions = TSSetFromOptions_GLEE; ts->ops->getstages = TSGetStages_GLEE; ts->ops->snesfunction = SNESTSFormFunction_GLEE; ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; ts->ops->getauxsolution = TSGetAuxSolution_GLEE; ts->ops->gettimeerror = TSGetTimeError_GLEE; ts->ops->settimeerror = TSSetTimeError_GLEE; ts->ops->startingmethod = TSStartingMethod_GLEE; ts->default_adapt_type = TSADAPTGLEE; ts->usessnes = PETSC_TRUE; PetscCall(PetscNew(&th)); ts->data = (void *)th; PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", TSGLEEGetType_GLEE)); PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", TSGLEESetType_GLEE)); PetscFunctionReturn(PETSC_SUCCESS); }