1211a84d6SLisandro Dalcin /* 2211a84d6SLisandro Dalcin Code for timestepping with BDF methods 3211a84d6SLisandro Dalcin */ 4211a84d6SLisandro Dalcin #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 51117961dSLisandro Dalcin #include <petscdm.h> 6211a84d6SLisandro Dalcin 7211a84d6SLisandro Dalcin static PetscBool cited = PETSC_FALSE; 89371c9d4SSatish Balay static const char citation[] = "@book{Brenan1995,\n" 9211a84d6SLisandro Dalcin " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n" 10211a84d6SLisandro Dalcin " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n" 11211a84d6SLisandro Dalcin " publisher = {Society for Industrial and Applied Mathematics},\n" 12211a84d6SLisandro Dalcin " year = {1995},\n" 13211a84d6SLisandro Dalcin " doi = {10.1137/1.9781611971224},\n}\n"; 14211a84d6SLisandro Dalcin 15211a84d6SLisandro Dalcin typedef struct { 16211a84d6SLisandro Dalcin PetscInt k, n; 17211a84d6SLisandro Dalcin PetscReal time[6 + 2]; 18211a84d6SLisandro Dalcin Vec work[6 + 2]; 19e3c11fc1SJed Brown Vec tvwork[6 + 2]; 20211a84d6SLisandro Dalcin PetscReal shift; 21e3c11fc1SJed Brown Vec vec_dot; /* Xdot when !transientvar, else Cdot where C(X) is the transient variable. */ 221117961dSLisandro Dalcin Vec vec_wrk; 23211a84d6SLisandro Dalcin Vec vec_lte; 24211a84d6SLisandro Dalcin 25e3c11fc1SJed Brown PetscBool transientvar; 26211a84d6SLisandro Dalcin PetscInt order; 27211a84d6SLisandro Dalcin TSStepStatus status; 28211a84d6SLisandro Dalcin } TS_BDF; 29211a84d6SLisandro Dalcin 30e3c11fc1SJed Brown /* Compute Lagrange polynomials on T[:n] evaluated at t. 31e3c11fc1SJed Brown * If one has data (T[i], Y[i]), then the interpolation/extrapolation is f(t) = \sum_i L[i]*Y[i]. 32e3c11fc1SJed Brown */ 339371c9d4SSatish Balay static inline void LagrangeBasisVals(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar L[]) { 34211a84d6SLisandro Dalcin PetscInt k, j; 35211a84d6SLisandro Dalcin for (k = 0; k < n; k++) 36211a84d6SLisandro Dalcin for (L[k] = 1, j = 0; j < n; j++) 379371c9d4SSatish Balay if (j != k) L[k] *= (t - T[j]) / (T[k] - T[j]); 38211a84d6SLisandro Dalcin } 39211a84d6SLisandro Dalcin 409371c9d4SSatish Balay static inline void LagrangeBasisDers(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar dL[]) { 41211a84d6SLisandro Dalcin PetscInt k, j, i; 42211a84d6SLisandro Dalcin for (k = 0; k < n; k++) 43211a84d6SLisandro Dalcin for (dL[k] = 0, j = 0; j < n; j++) 44211a84d6SLisandro Dalcin if (j != k) { 45211a84d6SLisandro Dalcin PetscReal L = 1 / (T[k] - T[j]); 46211a84d6SLisandro Dalcin for (i = 0; i < n; i++) 479371c9d4SSatish Balay if (i != j && i != k) L *= (t - T[i]) / (T[k] - T[i]); 48211a84d6SLisandro Dalcin dL[k] += L; 49211a84d6SLisandro Dalcin } 50211a84d6SLisandro Dalcin } 51211a84d6SLisandro Dalcin 529371c9d4SSatish Balay static PetscErrorCode TSBDF_GetVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot) { 531117961dSLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 541117961dSLisandro Dalcin 551117961dSLisandro Dalcin PetscFunctionBegin; 561117961dSLisandro Dalcin if (dm && dm != ts->dm) { 579566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot)); 589566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot)); 591117961dSLisandro Dalcin } else { 601117961dSLisandro Dalcin *Xdot = bdf->vec_dot; 611117961dSLisandro Dalcin *Ydot = bdf->vec_wrk; 621117961dSLisandro Dalcin } 631117961dSLisandro Dalcin PetscFunctionReturn(0); 641117961dSLisandro Dalcin } 651117961dSLisandro Dalcin 669371c9d4SSatish Balay static PetscErrorCode TSBDF_RestoreVecs(TS ts, DM dm, Vec *Xdot, Vec *Ydot) { 671117961dSLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 681117961dSLisandro Dalcin 691117961dSLisandro Dalcin PetscFunctionBegin; 701117961dSLisandro Dalcin if (dm && dm != ts->dm) { 719566063dSJacob Faibussowitsch PetscCall(DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Xdot", Xdot)); 729566063dSJacob Faibussowitsch PetscCall(DMRestoreNamedGlobalVector(dm, "TSBDF_Vec_Ydot", Ydot)); 731117961dSLisandro Dalcin } else { 743c633725SBarry Smith PetscCheck(*Xdot == bdf->vec_dot, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Vec does not match the cache"); 753c633725SBarry Smith PetscCheck(*Ydot == bdf->vec_wrk, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "Vec does not match the cache"); 761117961dSLisandro Dalcin *Xdot = NULL; 771117961dSLisandro Dalcin *Ydot = NULL; 781117961dSLisandro Dalcin } 791117961dSLisandro Dalcin PetscFunctionReturn(0); 801117961dSLisandro Dalcin } 811117961dSLisandro Dalcin 829371c9d4SSatish Balay static PetscErrorCode DMCoarsenHook_TSBDF(DM fine, DM coarse, void *ctx) { 831117961dSLisandro Dalcin PetscFunctionBegin; 841117961dSLisandro Dalcin PetscFunctionReturn(0); 851117961dSLisandro Dalcin } 861117961dSLisandro Dalcin 879371c9d4SSatish Balay static PetscErrorCode DMRestrictHook_TSBDF(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) { 881117961dSLisandro Dalcin TS ts = (TS)ctx; 891117961dSLisandro Dalcin Vec Ydot, Ydot_c; 901117961dSLisandro Dalcin Vec Xdot, Xdot_c; 911117961dSLisandro Dalcin 921117961dSLisandro Dalcin PetscFunctionBegin; 939566063dSJacob Faibussowitsch PetscCall(TSBDF_GetVecs(ts, fine, &Xdot, &Ydot)); 949566063dSJacob Faibussowitsch PetscCall(TSBDF_GetVecs(ts, coarse, &Xdot_c, &Ydot_c)); 951117961dSLisandro Dalcin 969566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Ydot, Ydot_c)); 979566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c)); 981117961dSLisandro Dalcin 999566063dSJacob Faibussowitsch PetscCall(TSBDF_RestoreVecs(ts, fine, &Xdot, &Ydot)); 1009566063dSJacob Faibussowitsch PetscCall(TSBDF_RestoreVecs(ts, coarse, &Xdot_c, &Ydot_c)); 1011117961dSLisandro Dalcin PetscFunctionReturn(0); 1021117961dSLisandro Dalcin } 1031117961dSLisandro Dalcin 1049371c9d4SSatish Balay static PetscErrorCode TSBDF_Advance(TS ts, PetscReal t, Vec X) { 105211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 106211a84d6SLisandro Dalcin PetscInt i, n = (PetscInt)(sizeof(bdf->work) / sizeof(Vec)); 107e3c11fc1SJed Brown Vec tail = bdf->work[n - 1], tvtail = bdf->tvwork[n - 1]; 108211a84d6SLisandro Dalcin 109211a84d6SLisandro Dalcin PetscFunctionBegin; 110211a84d6SLisandro Dalcin for (i = n - 1; i >= 2; i--) { 111211a84d6SLisandro Dalcin bdf->time[i] = bdf->time[i - 1]; 112211a84d6SLisandro Dalcin bdf->work[i] = bdf->work[i - 1]; 113e3c11fc1SJed Brown bdf->tvwork[i] = bdf->tvwork[i - 1]; 114211a84d6SLisandro Dalcin } 115211a84d6SLisandro Dalcin bdf->n = PetscMin(bdf->n + 1, n - 1); 116211a84d6SLisandro Dalcin bdf->time[1] = t; 117211a84d6SLisandro Dalcin bdf->work[1] = tail; 118e3c11fc1SJed Brown bdf->tvwork[1] = tvtail; 1199566063dSJacob Faibussowitsch PetscCall(VecCopy(X, tail)); 1209566063dSJacob Faibussowitsch PetscCall(TSComputeTransientVariable(ts, tail, tvtail)); 121211a84d6SLisandro Dalcin PetscFunctionReturn(0); 122211a84d6SLisandro Dalcin } 123211a84d6SLisandro Dalcin 1249371c9d4SSatish Balay static PetscErrorCode TSBDF_VecLTE(TS ts, PetscInt order, Vec lte) { 125211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 126211a84d6SLisandro Dalcin PetscInt i, n = order + 1; 127211a84d6SLisandro Dalcin PetscReal *time = bdf->time; 128211a84d6SLisandro Dalcin Vec *vecs = bdf->work; 129211a84d6SLisandro Dalcin PetscScalar a[8], b[8], alpha[8]; 130211a84d6SLisandro Dalcin 131211a84d6SLisandro Dalcin PetscFunctionBegin; 1329371c9d4SSatish Balay LagrangeBasisDers(n + 0, time[0], time, a); 1339371c9d4SSatish Balay a[n] = 0; 134211a84d6SLisandro Dalcin LagrangeBasisDers(n + 1, time[0], time, b); 135211a84d6SLisandro Dalcin for (i = 0; i < n + 1; i++) alpha[i] = (a[i] - b[i]) / a[0]; 1369566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(lte)); 1379566063dSJacob Faibussowitsch PetscCall(VecMAXPY(lte, n + 1, alpha, vecs)); 138211a84d6SLisandro Dalcin PetscFunctionReturn(0); 139211a84d6SLisandro Dalcin } 140211a84d6SLisandro Dalcin 1419371c9d4SSatish Balay static PetscErrorCode TSBDF_Extrapolate(TS ts, PetscInt order, PetscReal t, Vec X) { 142211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 143211a84d6SLisandro Dalcin PetscInt n = order + 1; 144211a84d6SLisandro Dalcin PetscReal *time = bdf->time + 1; 145211a84d6SLisandro Dalcin Vec *vecs = bdf->work + 1; 146211a84d6SLisandro Dalcin PetscScalar alpha[7]; 147211a84d6SLisandro Dalcin 148211a84d6SLisandro Dalcin PetscFunctionBegin; 149211a84d6SLisandro Dalcin n = PetscMin(n, bdf->n); 150211a84d6SLisandro Dalcin LagrangeBasisVals(n, t, time, alpha); 1519566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(X)); 1529566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, n, alpha, vecs)); 153211a84d6SLisandro Dalcin PetscFunctionReturn(0); 154211a84d6SLisandro Dalcin } 155211a84d6SLisandro Dalcin 1569371c9d4SSatish Balay static PetscErrorCode TSBDF_Interpolate(TS ts, PetscInt order, PetscReal t, Vec X) { 157211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 158211a84d6SLisandro Dalcin PetscInt n = order + 1; 159211a84d6SLisandro Dalcin PetscReal *time = bdf->time; 160211a84d6SLisandro Dalcin Vec *vecs = bdf->work; 161211a84d6SLisandro Dalcin PetscScalar alpha[7]; 162211a84d6SLisandro Dalcin 163211a84d6SLisandro Dalcin PetscFunctionBegin; 164211a84d6SLisandro Dalcin LagrangeBasisVals(n, t, time, alpha); 1659566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(X)); 1669566063dSJacob Faibussowitsch PetscCall(VecMAXPY(X, n, alpha, vecs)); 167211a84d6SLisandro Dalcin PetscFunctionReturn(0); 168211a84d6SLisandro Dalcin } 169211a84d6SLisandro Dalcin 170e3c11fc1SJed Brown /* Compute the affine term V0 such that Xdot = shift*X + V0. 171e3c11fc1SJed Brown * 172e3c11fc1SJed Brown * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork. 173e3c11fc1SJed Brown */ 1749371c9d4SSatish Balay static PetscErrorCode TSBDF_PreSolve(TS ts) { 1751117961dSLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 1761117961dSLisandro Dalcin PetscInt i, n = PetscMax(bdf->k, 1) + 1; 1771117961dSLisandro Dalcin Vec V, V0; 1781117961dSLisandro Dalcin Vec vecs[7]; 1791117961dSLisandro Dalcin PetscScalar alpha[7]; 1801117961dSLisandro Dalcin 1811117961dSLisandro Dalcin PetscFunctionBegin; 1829566063dSJacob Faibussowitsch PetscCall(TSBDF_GetVecs(ts, NULL, &V, &V0)); 1831117961dSLisandro Dalcin LagrangeBasisDers(n, bdf->time[0], bdf->time, alpha); 1849371c9d4SSatish Balay for (i = 1; i < n; i++) { vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i]; } 1859566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(V0)); 1869566063dSJacob Faibussowitsch PetscCall(VecMAXPY(V0, n - 1, alpha + 1, vecs + 1)); 1871117961dSLisandro Dalcin bdf->shift = PetscRealPart(alpha[0]); 1889566063dSJacob Faibussowitsch PetscCall(TSBDF_RestoreVecs(ts, NULL, &V, &V0)); 1891117961dSLisandro Dalcin PetscFunctionReturn(0); 1901117961dSLisandro Dalcin } 1911117961dSLisandro Dalcin 1929371c9d4SSatish Balay static PetscErrorCode TSBDF_SNESSolve(TS ts, Vec b, Vec x) { 193211a84d6SLisandro Dalcin PetscInt nits, lits; 194211a84d6SLisandro Dalcin 195211a84d6SLisandro Dalcin PetscFunctionBegin; 1969566063dSJacob Faibussowitsch PetscCall(TSBDF_PreSolve(ts)); 1979566063dSJacob Faibussowitsch PetscCall(SNESSolve(ts->snes, b, x)); 1989566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(ts->snes, &nits)); 1999566063dSJacob Faibussowitsch PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits)); 2009371c9d4SSatish Balay ts->snes_its += nits; 2019371c9d4SSatish Balay ts->ksp_its += lits; 202211a84d6SLisandro Dalcin PetscFunctionReturn(0); 203211a84d6SLisandro Dalcin } 204211a84d6SLisandro Dalcin 2059371c9d4SSatish Balay static PetscErrorCode TSBDF_Restart(TS ts, PetscBool *accept) { 206211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 207211a84d6SLisandro Dalcin 208211a84d6SLisandro Dalcin PetscFunctionBegin; 2099371c9d4SSatish Balay bdf->k = 1; 2109371c9d4SSatish Balay bdf->n = 0; 2119566063dSJacob Faibussowitsch PetscCall(TSBDF_Advance(ts, ts->ptime, ts->vec_sol)); 212211a84d6SLisandro Dalcin 213211a84d6SLisandro Dalcin bdf->time[0] = ts->ptime + ts->time_step / 2; 2149566063dSJacob Faibussowitsch PetscCall(VecCopy(bdf->work[1], bdf->work[0])); 2159566063dSJacob Faibussowitsch PetscCall(TSPreStage(ts, bdf->time[0])); 2169566063dSJacob Faibussowitsch PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0])); 2179566063dSJacob Faibussowitsch PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0])); 2189566063dSJacob Faibussowitsch PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], accept)); 219211a84d6SLisandro Dalcin if (!*accept) PetscFunctionReturn(0); 220211a84d6SLisandro Dalcin 2219371c9d4SSatish Balay bdf->k = PetscMin(2, bdf->order); 2229371c9d4SSatish Balay bdf->n++; 2239566063dSJacob Faibussowitsch PetscCall(VecCopy(bdf->work[0], bdf->work[2])); 224211a84d6SLisandro Dalcin bdf->time[2] = bdf->time[0]; 2259566063dSJacob Faibussowitsch PetscCall(TSComputeTransientVariable(ts, bdf->work[2], bdf->tvwork[2])); 226211a84d6SLisandro Dalcin PetscFunctionReturn(0); 227211a84d6SLisandro Dalcin } 228211a84d6SLisandro Dalcin 229211a84d6SLisandro Dalcin static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"}; 230211a84d6SLisandro Dalcin 2319371c9d4SSatish Balay static PetscErrorCode TSStep_BDF(TS ts) { 232211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 233211a84d6SLisandro Dalcin PetscInt rejections = 0; 234211a84d6SLisandro Dalcin PetscBool stageok, accept = PETSC_TRUE; 235211a84d6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 236211a84d6SLisandro Dalcin 237211a84d6SLisandro Dalcin PetscFunctionBegin; 2389566063dSJacob Faibussowitsch PetscCall(PetscCitationsRegister(citation, &cited)); 239211a84d6SLisandro Dalcin 240211a84d6SLisandro Dalcin if (!ts->steprollback && !ts->steprestart) { 241211a84d6SLisandro Dalcin bdf->k = PetscMin(bdf->k + 1, bdf->order); 2429566063dSJacob Faibussowitsch PetscCall(TSBDF_Advance(ts, ts->ptime, ts->vec_sol)); 243211a84d6SLisandro Dalcin } 244211a84d6SLisandro Dalcin 245211a84d6SLisandro Dalcin bdf->status = TS_STEP_INCOMPLETE; 246211a84d6SLisandro Dalcin while (!ts->reason && bdf->status != TS_STEP_COMPLETE) { 247211a84d6SLisandro Dalcin if (ts->steprestart) { 2489566063dSJacob Faibussowitsch PetscCall(TSBDF_Restart(ts, &stageok)); 249211a84d6SLisandro Dalcin if (!stageok) goto reject_step; 250211a84d6SLisandro Dalcin } 251211a84d6SLisandro Dalcin 252211a84d6SLisandro Dalcin bdf->time[0] = ts->ptime + ts->time_step; 2539566063dSJacob Faibussowitsch PetscCall(TSBDF_Extrapolate(ts, bdf->k - (accept ? 0 : 1), bdf->time[0], bdf->work[0])); 2549566063dSJacob Faibussowitsch PetscCall(TSPreStage(ts, bdf->time[0])); 2559566063dSJacob Faibussowitsch PetscCall(TSBDF_SNESSolve(ts, NULL, bdf->work[0])); 2569566063dSJacob Faibussowitsch PetscCall(TSPostStage(ts, bdf->time[0], 0, &bdf->work[0])); 2579566063dSJacob Faibussowitsch PetscCall(TSAdaptCheckStage(ts->adapt, ts, bdf->time[0], bdf->work[0], &stageok)); 258211a84d6SLisandro Dalcin if (!stageok) goto reject_step; 259211a84d6SLisandro Dalcin 260211a84d6SLisandro Dalcin bdf->status = TS_STEP_PENDING; 2619566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidatesClear(ts->adapt)); 2629566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidateAdd(ts->adapt, BDF_SchemeName[bdf->k], bdf->k, 1, 1.0, 1.0, PETSC_TRUE)); 2639566063dSJacob Faibussowitsch PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 264211a84d6SLisandro Dalcin bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 2659371c9d4SSatish Balay if (!accept) { 2669371c9d4SSatish Balay ts->time_step = next_time_step; 2679371c9d4SSatish Balay goto reject_step; 2689371c9d4SSatish Balay } 269211a84d6SLisandro Dalcin 2709566063dSJacob Faibussowitsch PetscCall(VecCopy(bdf->work[0], ts->vec_sol)); 271211a84d6SLisandro Dalcin ts->ptime += ts->time_step; 272211a84d6SLisandro Dalcin ts->time_step = next_time_step; 273211a84d6SLisandro Dalcin break; 274211a84d6SLisandro Dalcin 275211a84d6SLisandro Dalcin reject_step: 2769371c9d4SSatish Balay ts->reject++; 2779371c9d4SSatish Balay accept = PETSC_FALSE; 278211a84d6SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 27963a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 280211a84d6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 281211a84d6SLisandro Dalcin } 282211a84d6SLisandro Dalcin } 283211a84d6SLisandro Dalcin PetscFunctionReturn(0); 284211a84d6SLisandro Dalcin } 285211a84d6SLisandro Dalcin 2869371c9d4SSatish Balay static PetscErrorCode TSInterpolate_BDF(TS ts, PetscReal t, Vec X) { 287211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 288211a84d6SLisandro Dalcin 289211a84d6SLisandro Dalcin PetscFunctionBegin; 2909566063dSJacob Faibussowitsch PetscCall(TSBDF_Interpolate(ts, bdf->k, t, X)); 291211a84d6SLisandro Dalcin PetscFunctionReturn(0); 292211a84d6SLisandro Dalcin } 293211a84d6SLisandro Dalcin 2949371c9d4SSatish Balay static PetscErrorCode TSEvaluateWLTE_BDF(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte) { 295211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 296211a84d6SLisandro Dalcin PetscInt k = bdf->k; 2977453f775SEmil Constantinescu PetscReal wltea, wlter; 298211a84d6SLisandro Dalcin Vec X = bdf->work[0], Y = bdf->vec_lte; 299211a84d6SLisandro Dalcin 300211a84d6SLisandro Dalcin PetscFunctionBegin; 301211a84d6SLisandro Dalcin k = PetscMin(k, bdf->n - 1); 3029566063dSJacob Faibussowitsch PetscCall(TSBDF_VecLTE(ts, k, Y)); 3039566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y, 1, X)); 3049566063dSJacob Faibussowitsch PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter)); 305211a84d6SLisandro Dalcin if (order) *order = k + 1; 306211a84d6SLisandro Dalcin PetscFunctionReturn(0); 307211a84d6SLisandro Dalcin } 308211a84d6SLisandro Dalcin 3099371c9d4SSatish Balay static PetscErrorCode TSRollBack_BDF(TS ts) { 310211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 311211a84d6SLisandro Dalcin 312211a84d6SLisandro Dalcin PetscFunctionBegin; 3139566063dSJacob Faibussowitsch PetscCall(VecCopy(bdf->work[1], ts->vec_sol)); 314211a84d6SLisandro Dalcin PetscFunctionReturn(0); 315211a84d6SLisandro Dalcin } 316211a84d6SLisandro Dalcin 3179371c9d4SSatish Balay static PetscErrorCode SNESTSFormFunction_BDF(SNES snes, Vec X, Vec F, TS ts) { 318211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 3191117961dSLisandro Dalcin DM dm, dmsave = ts->dm; 320211a84d6SLisandro Dalcin PetscReal t = bdf->time[0]; 3211117961dSLisandro Dalcin PetscReal shift = bdf->shift; 3221117961dSLisandro Dalcin Vec V, V0; 323211a84d6SLisandro Dalcin 324211a84d6SLisandro Dalcin PetscFunctionBegin; 3259566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 3269566063dSJacob Faibussowitsch PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0)); 327e3c11fc1SJed Brown if (bdf->transientvar) { /* shift*C(X) + V0 */ 3289566063dSJacob Faibussowitsch PetscCall(TSComputeTransientVariable(ts, X, V)); 3299566063dSJacob Faibussowitsch PetscCall(VecAYPX(V, shift, V0)); 330e3c11fc1SJed Brown } else { /* shift*X + V0 */ 3319566063dSJacob Faibussowitsch PetscCall(VecWAXPY(V, shift, X, V0)); 332e3c11fc1SJed Brown } 3331117961dSLisandro Dalcin 334211a84d6SLisandro Dalcin /* F = Function(t,X,V) */ 3351117961dSLisandro Dalcin ts->dm = dm; 3369566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, t, X, V, F, PETSC_FALSE)); 3371117961dSLisandro Dalcin ts->dm = dmsave; 3381117961dSLisandro Dalcin 3399566063dSJacob Faibussowitsch PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0)); 340211a84d6SLisandro Dalcin PetscFunctionReturn(0); 341211a84d6SLisandro Dalcin } 342211a84d6SLisandro Dalcin 3439371c9d4SSatish Balay static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes, Vec X, Mat J, Mat P, TS ts) { 344211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 3451117961dSLisandro Dalcin DM dm, dmsave = ts->dm; 346211a84d6SLisandro Dalcin PetscReal t = bdf->time[0]; 3471117961dSLisandro Dalcin PetscReal shift = bdf->shift; 3481117961dSLisandro Dalcin Vec V, V0; 349211a84d6SLisandro Dalcin 350211a84d6SLisandro Dalcin PetscFunctionBegin; 3519566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 3529566063dSJacob Faibussowitsch PetscCall(TSBDF_GetVecs(ts, dm, &V, &V0)); 3531117961dSLisandro Dalcin 354211a84d6SLisandro Dalcin /* J,P = Jacobian(t,X,V) */ 3551117961dSLisandro Dalcin ts->dm = dm; 3569566063dSJacob Faibussowitsch PetscCall(TSComputeIJacobian(ts, t, X, V, shift, J, P, PETSC_FALSE)); 3571117961dSLisandro Dalcin ts->dm = dmsave; 3581117961dSLisandro Dalcin 3599566063dSJacob Faibussowitsch PetscCall(TSBDF_RestoreVecs(ts, dm, &V, &V0)); 360211a84d6SLisandro Dalcin PetscFunctionReturn(0); 361211a84d6SLisandro Dalcin } 362211a84d6SLisandro Dalcin 3639371c9d4SSatish Balay static PetscErrorCode TSReset_BDF(TS ts) { 364211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 365211a84d6SLisandro Dalcin size_t i, n = sizeof(bdf->work) / sizeof(Vec); 366211a84d6SLisandro Dalcin 367211a84d6SLisandro Dalcin PetscFunctionBegin; 3681117961dSLisandro Dalcin bdf->k = bdf->n = 0; 369e3c11fc1SJed Brown for (i = 0; i < n; i++) { 3709566063dSJacob Faibussowitsch PetscCall(VecDestroy(&bdf->work[i])); 3719566063dSJacob Faibussowitsch PetscCall(VecDestroy(&bdf->tvwork[i])); 372e3c11fc1SJed Brown } 3739566063dSJacob Faibussowitsch PetscCall(VecDestroy(&bdf->vec_dot)); 3749566063dSJacob Faibussowitsch PetscCall(VecDestroy(&bdf->vec_wrk)); 3759566063dSJacob Faibussowitsch PetscCall(VecDestroy(&bdf->vec_lte)); 3769566063dSJacob Faibussowitsch if (ts->dm) PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts)); 377211a84d6SLisandro Dalcin PetscFunctionReturn(0); 378211a84d6SLisandro Dalcin } 379211a84d6SLisandro Dalcin 3809371c9d4SSatish Balay static PetscErrorCode TSDestroy_BDF(TS ts) { 381211a84d6SLisandro Dalcin PetscFunctionBegin; 3829566063dSJacob Faibussowitsch PetscCall(TSReset_BDF(ts)); 3839566063dSJacob Faibussowitsch PetscCall(PetscFree(ts->data)); 3849566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", NULL)); 3859566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", NULL)); 386211a84d6SLisandro Dalcin PetscFunctionReturn(0); 387211a84d6SLisandro Dalcin } 388211a84d6SLisandro Dalcin 3899371c9d4SSatish Balay static PetscErrorCode TSSetUp_BDF(TS ts) { 390211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 391211a84d6SLisandro Dalcin size_t i, n = sizeof(bdf->work) / sizeof(Vec); 3922ffb9264SLisandro Dalcin PetscReal low, high, two = 2; 393211a84d6SLisandro Dalcin 394211a84d6SLisandro Dalcin PetscFunctionBegin; 3959566063dSJacob Faibussowitsch PetscCall(TSHasTransientVariable(ts, &bdf->transientvar)); 396211a84d6SLisandro Dalcin bdf->k = bdf->n = 0; 397e3c11fc1SJed Brown for (i = 0; i < n; i++) { 3989566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &bdf->work[i])); 399*48a46eb9SPierre Jolivet if (i && bdf->transientvar) PetscCall(VecDuplicate(ts->vec_sol, &bdf->tvwork[i])); 400e3c11fc1SJed Brown } 4019566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_dot)); 4029566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_wrk)); 4039566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &bdf->vec_lte)); 4049566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &ts->dm)); 4059566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookAdd(ts->dm, DMCoarsenHook_TSBDF, DMRestrictHook_TSBDF, ts)); 406211a84d6SLisandro Dalcin 4079566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &ts->adapt)); 4089566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidatesClear(ts->adapt)); 4099566063dSJacob Faibussowitsch PetscCall(TSAdaptGetClip(ts->adapt, &low, &high)); 4109566063dSJacob Faibussowitsch PetscCall(TSAdaptSetClip(ts->adapt, low, PetscMin(high, two))); 411211a84d6SLisandro Dalcin 4129566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &ts->snes)); 413211a84d6SLisandro Dalcin PetscFunctionReturn(0); 414211a84d6SLisandro Dalcin } 415211a84d6SLisandro Dalcin 4169371c9d4SSatish Balay static PetscErrorCode TSSetFromOptions_BDF(TS ts, PetscOptionItems *PetscOptionsObject) { 417211a84d6SLisandro Dalcin PetscFunctionBegin; 418d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "BDF ODE solver options"); 419211a84d6SLisandro Dalcin { 420211a84d6SLisandro Dalcin PetscBool flg; 421e5b8ffdfSLisandro Dalcin PetscInt order; 4229566063dSJacob Faibussowitsch PetscCall(TSBDFGetOrder(ts, &order)); 4239566063dSJacob Faibussowitsch PetscCall(PetscOptionsInt("-ts_bdf_order", "Order of the BDF method", "TSBDFSetOrder", order, &order, &flg)); 4249566063dSJacob Faibussowitsch if (flg) PetscCall(TSBDFSetOrder(ts, order)); 425211a84d6SLisandro Dalcin } 426d0609cedSBarry Smith PetscOptionsHeadEnd(); 427211a84d6SLisandro Dalcin PetscFunctionReturn(0); 428211a84d6SLisandro Dalcin } 429211a84d6SLisandro Dalcin 4309371c9d4SSatish Balay static PetscErrorCode TSView_BDF(TS ts, PetscViewer viewer) { 431211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 432211a84d6SLisandro Dalcin PetscBool iascii; 433211a84d6SLisandro Dalcin 434211a84d6SLisandro Dalcin PetscFunctionBegin; 4359566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 436*48a46eb9SPierre Jolivet if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Order=%" PetscInt_FMT "\n", bdf->order)); 437211a84d6SLisandro Dalcin PetscFunctionReturn(0); 438211a84d6SLisandro Dalcin } 439211a84d6SLisandro Dalcin 440211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */ 441211a84d6SLisandro Dalcin 4429371c9d4SSatish Balay static PetscErrorCode TSBDFSetOrder_BDF(TS ts, PetscInt order) { 443211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 444211a84d6SLisandro Dalcin 445211a84d6SLisandro Dalcin PetscFunctionBegin; 446211a84d6SLisandro Dalcin if (order == bdf->order) PetscFunctionReturn(0); 44763a3b9bcSJacob Faibussowitsch PetscCheck(order >= 1 && order <= 6, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "BDF Order %" PetscInt_FMT " not implemented", order); 448211a84d6SLisandro Dalcin bdf->order = order; 449211a84d6SLisandro Dalcin PetscFunctionReturn(0); 450211a84d6SLisandro Dalcin } 451211a84d6SLisandro Dalcin 4529371c9d4SSatish Balay static PetscErrorCode TSBDFGetOrder_BDF(TS ts, PetscInt *order) { 453211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF *)ts->data; 454211a84d6SLisandro Dalcin 455211a84d6SLisandro Dalcin PetscFunctionBegin; 456211a84d6SLisandro Dalcin *order = bdf->order; 457211a84d6SLisandro Dalcin PetscFunctionReturn(0); 458211a84d6SLisandro Dalcin } 459211a84d6SLisandro Dalcin 460211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */ 461211a84d6SLisandro Dalcin 462211a84d6SLisandro Dalcin /*MC 463211a84d6SLisandro Dalcin TSBDF - DAE solver using BDF methods 464211a84d6SLisandro Dalcin 465211a84d6SLisandro Dalcin Level: beginner 466211a84d6SLisandro Dalcin 467db781477SPatrick Sanan .seealso: `TS`, `TSCreate()`, `TSSetType()` 468211a84d6SLisandro Dalcin M*/ 4699371c9d4SSatish Balay PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts) { 470211a84d6SLisandro Dalcin TS_BDF *bdf; 471211a84d6SLisandro Dalcin 472211a84d6SLisandro Dalcin PetscFunctionBegin; 473211a84d6SLisandro Dalcin ts->ops->reset = TSReset_BDF; 474211a84d6SLisandro Dalcin ts->ops->destroy = TSDestroy_BDF; 475211a84d6SLisandro Dalcin ts->ops->view = TSView_BDF; 476211a84d6SLisandro Dalcin ts->ops->setup = TSSetUp_BDF; 477211a84d6SLisandro Dalcin ts->ops->setfromoptions = TSSetFromOptions_BDF; 478211a84d6SLisandro Dalcin ts->ops->step = TSStep_BDF; 479211a84d6SLisandro Dalcin ts->ops->evaluatewlte = TSEvaluateWLTE_BDF; 480211a84d6SLisandro Dalcin ts->ops->rollback = TSRollBack_BDF; 481211a84d6SLisandro Dalcin ts->ops->interpolate = TSInterpolate_BDF; 482211a84d6SLisandro Dalcin ts->ops->snesfunction = SNESTSFormFunction_BDF; 483211a84d6SLisandro Dalcin ts->ops->snesjacobian = SNESTSFormJacobian_BDF; 4842ffb9264SLisandro Dalcin ts->default_adapt_type = TSADAPTBASIC; 485211a84d6SLisandro Dalcin 486efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 487efd4aadfSBarry Smith 4889566063dSJacob Faibussowitsch PetscCall(PetscNewLog(ts, &bdf)); 489211a84d6SLisandro Dalcin ts->data = (void *)bdf; 490211a84d6SLisandro Dalcin 491211a84d6SLisandro Dalcin bdf->status = TS_STEP_COMPLETE; 492211a84d6SLisandro Dalcin 4939566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFSetOrder_C", TSBDFSetOrder_BDF)); 4949566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSBDFGetOrder_C", TSBDFGetOrder_BDF)); 4959566063dSJacob Faibussowitsch PetscCall(TSBDFSetOrder(ts, 2)); 496211a84d6SLisandro Dalcin PetscFunctionReturn(0); 497211a84d6SLisandro Dalcin } 498211a84d6SLisandro Dalcin 499211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */ 500211a84d6SLisandro Dalcin 501211a84d6SLisandro Dalcin /*@ 502211a84d6SLisandro Dalcin TSBDFSetOrder - Set the order of the BDF method 503211a84d6SLisandro Dalcin 504211a84d6SLisandro Dalcin Logically Collective on TS 505211a84d6SLisandro Dalcin 506d8d19677SJose E. Roman Input Parameters: 507211a84d6SLisandro Dalcin + ts - timestepping context 508211a84d6SLisandro Dalcin - order - order of the method 509211a84d6SLisandro Dalcin 510211a84d6SLisandro Dalcin Options Database: 511147403d9SBarry Smith . -ts_bdf_order <order> - select the order 512211a84d6SLisandro Dalcin 513211a84d6SLisandro Dalcin Level: intermediate 514211a84d6SLisandro Dalcin 515211a84d6SLisandro Dalcin @*/ 5169371c9d4SSatish Balay PetscErrorCode TSBDFSetOrder(TS ts, PetscInt order) { 517211a84d6SLisandro Dalcin PetscFunctionBegin; 518211a84d6SLisandro Dalcin PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 519211a84d6SLisandro Dalcin PetscValidLogicalCollectiveInt(ts, order, 2); 520cac4c232SBarry Smith PetscTryMethod(ts, "TSBDFSetOrder_C", (TS, PetscInt), (ts, order)); 521211a84d6SLisandro Dalcin PetscFunctionReturn(0); 522211a84d6SLisandro Dalcin } 523211a84d6SLisandro Dalcin 524211a84d6SLisandro Dalcin /*@ 525211a84d6SLisandro Dalcin TSBDFGetOrder - Get the order of the BDF method 526211a84d6SLisandro Dalcin 527211a84d6SLisandro Dalcin Not Collective 528211a84d6SLisandro Dalcin 529211a84d6SLisandro Dalcin Input Parameter: 530211a84d6SLisandro Dalcin . ts - timestepping context 531211a84d6SLisandro Dalcin 532211a84d6SLisandro Dalcin Output Parameter: 533211a84d6SLisandro Dalcin . order - order of the method 534211a84d6SLisandro Dalcin 535211a84d6SLisandro Dalcin Level: intermediate 536211a84d6SLisandro Dalcin 537211a84d6SLisandro Dalcin @*/ 5389371c9d4SSatish Balay PetscErrorCode TSBDFGetOrder(TS ts, PetscInt *order) { 539211a84d6SLisandro Dalcin PetscFunctionBegin; 540211a84d6SLisandro Dalcin PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 541211a84d6SLisandro Dalcin PetscValidIntPointer(order, 2); 542cac4c232SBarry Smith PetscUseMethod(ts, "TSBDFGetOrder_C", (TS, PetscInt *), (ts, order)); 543211a84d6SLisandro Dalcin PetscFunctionReturn(0); 544211a84d6SLisandro Dalcin } 545