Actual source code: glle.c
2: #include <../src/ts/impls/implicit/glle/glle.h>
3: #include <petscdm.h>
4: #include <petscblaslapack.h>
6: static const char *TSGLLEErrorDirections[] = {"FORWARD", "BACKWARD", "TSGLLEErrorDirection", "TSGLLEERROR_", NULL};
7: static PetscFunctionList TSGLLEList;
8: static PetscFunctionList TSGLLEAcceptList;
9: static PetscBool TSGLLEPackageInitialized;
10: static PetscBool TSGLLERegisterAllCalled;
12: /* This function is pure */
13: static PetscScalar Factorial(PetscInt n)
14: {
15: PetscInt i;
16: if (n < 12) { /* Can compute with 32-bit integers */
17: PetscInt f = 1;
18: for (i = 2; i <= n; i++) f *= i;
19: return (PetscScalar)f;
20: } else {
21: PetscScalar f = 1.;
22: for (i = 2; i <= n; i++) f *= (PetscScalar)i;
23: return f;
24: }
25: }
27: /* This function is pure */
28: static PetscScalar CPowF(PetscScalar c, PetscInt p)
29: {
30: return PetscPowRealInt(PetscRealPart(c), p) / Factorial(p);
31: }
33: static PetscErrorCode TSGLLEGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
34: {
35: TS_GLLE *gl = (TS_GLLE *)ts->data;
37: if (Z) {
38: if (dm && dm != ts->dm) {
39: DMGetNamedGlobalVector(dm, "TSGLLE_Z", Z);
40: } else *Z = gl->Z;
41: }
42: if (Ydotstage) {
43: if (dm && dm != ts->dm) {
44: DMGetNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage);
45: } else *Ydotstage = gl->Ydot[gl->stage];
46: }
47: return 0;
48: }
50: static PetscErrorCode TSGLLERestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
51: {
52: if (Z) {
53: if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSGLLE_Z", Z);
54: }
55: if (Ydotstage) {
56: if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage);
57: }
58: return 0;
59: }
61: static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, void *ctx)
62: {
63: return 0;
64: }
66: static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
67: {
68: TS ts = (TS)ctx;
69: Vec Ydot, Ydot_c;
71: TSGLLEGetVecs(ts, fine, NULL, &Ydot);
72: TSGLLEGetVecs(ts, coarse, NULL, &Ydot_c);
73: MatRestrict(restrct, Ydot, Ydot_c);
74: VecPointwiseMult(Ydot_c, rscale, Ydot_c);
75: TSGLLERestoreVecs(ts, fine, NULL, &Ydot);
76: TSGLLERestoreVecs(ts, coarse, NULL, &Ydot_c);
77: return 0;
78: }
80: static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, void *ctx)
81: {
82: return 0;
83: }
85: static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
86: {
87: TS ts = (TS)ctx;
88: Vec Ydot, Ydot_s;
90: TSGLLEGetVecs(ts, dm, NULL, &Ydot);
91: TSGLLEGetVecs(ts, subdm, NULL, &Ydot_s);
93: VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD);
94: VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD);
96: TSGLLERestoreVecs(ts, dm, NULL, &Ydot);
97: TSGLLERestoreVecs(ts, subdm, NULL, &Ydot_s);
98: return 0;
99: }
101: 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)
102: {
103: TSGLLEScheme scheme;
104: PetscInt j;
110: *inscheme = NULL;
111: PetscNew(&scheme);
112: scheme->p = p;
113: scheme->q = q;
114: scheme->r = r;
115: scheme->s = s;
117: PetscMalloc5(s, &scheme->c, s * s, &scheme->a, r * s, &scheme->b, r * s, &scheme->u, r * r, &scheme->v);
118: PetscArraycpy(scheme->c, c, s);
119: for (j = 0; j < s * s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j];
120: for (j = 0; j < r * s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j];
121: for (j = 0; j < s * r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j];
122: for (j = 0; j < r * r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j];
124: PetscMalloc6(r, &scheme->alpha, r, &scheme->beta, r, &scheme->gamma, 3 * s, &scheme->phi, 3 * r, &scheme->psi, r, &scheme->stage_error);
125: {
126: PetscInt i, j, k, ss = s + 2;
127: PetscBLASInt m, n, one = 1, *ipiv, lwork = 4 * ((s + 3) * 3 + 3), info, ldb;
128: PetscReal rcond, *sing, *workreal;
129: PetscScalar *ImV, *H, *bmat, *workscalar, *c = scheme->c, *a = scheme->a, *b = scheme->b, *u = scheme->u, *v = scheme->v;
130: PetscBLASInt rank;
131: PetscMalloc7(PetscSqr(r), &ImV, 3 * s, &H, 3 * ss, &bmat, lwork, &workscalar, 5 * (3 + r), &workreal, r + s, &sing, r + s, &ipiv);
133: /* column-major input */
134: for (i = 0; i < r - 1; i++) {
135: for (j = 0; j < r - 1; j++) ImV[i + j * r] = 1.0 * (i == j) - v[(i + 1) * r + j + 1];
136: }
137: /* Build right hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */
138: for (i = 1; i < r; i++) {
139: scheme->alpha[i] = 1. / Factorial(p + 1 - i);
140: for (j = 0; j < s; j++) scheme->alpha[i] -= b[i * s + j] * CPowF(c[j], p);
141: }
142: PetscBLASIntCast(r - 1, &m);
143: PetscBLASIntCast(r, &n);
144: PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&m, &one, ImV, &n, ipiv, scheme->alpha + 1, &n, &info));
148: /* Build right hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */
149: for (i = 1; i < r; i++) {
150: scheme->beta[i] = 1. / Factorial(p + 2 - i) - scheme->alpha[i];
151: for (j = 0; j < s; j++) scheme->beta[i] -= b[i * s + j] * CPowF(c[j], p + 1);
152: }
153: PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->beta + 1, &n, &info));
157: /* Build stage_error vector
158: xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha;
159: */
160: for (i = 0; i < s; i++) {
161: scheme->stage_error[i] = CPowF(c[i], p + 1);
162: for (j = 0; j < s; j++) scheme->stage_error[i] -= a[i * s + j] * CPowF(c[j], p);
163: for (j = 1; j < r; j++) scheme->stage_error[i] += u[i * r + j] * scheme->alpha[j];
164: }
166: /* alpha[0] (epsilon in B,J,W 2007)
167: epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha;
168: */
169: scheme->alpha[0] = 1. / Factorial(p + 1);
170: for (j = 0; j < s; j++) scheme->alpha[0] -= b[0 * s + j] * CPowF(c[j], p);
171: for (j = 1; j < r; j++) scheme->alpha[0] += v[0 * r + j] * scheme->alpha[j];
173: /* right hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */
174: for (i = 1; i < r; i++) {
175: scheme->gamma[i] = (i == 1 ? -1. : 0) * scheme->alpha[0];
176: for (j = 0; j < s; j++) scheme->gamma[i] += b[i * s + j] * scheme->stage_error[j];
177: }
178: PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->gamma + 1, &n, &info));
182: /* beta[0] (rho in B,J,W 2007)
183: e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ...
184: + glm.V(1,2:end)*e.beta;% - e.epsilon;
185: % Note: The paper (B,J,W 2007) includes the last term in their definition
186: * */
187: scheme->beta[0] = 1. / Factorial(p + 2);
188: for (j = 0; j < s; j++) scheme->beta[0] -= b[0 * s + j] * CPowF(c[j], p + 1);
189: for (j = 1; j < r; j++) scheme->beta[0] += v[0 * r + j] * scheme->beta[j];
191: /* gamma[0] (sigma in B,J,W 2007)
192: * e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma;
193: * */
194: scheme->gamma[0] = 0.0;
195: for (j = 0; j < s; j++) scheme->gamma[0] += b[0 * s + j] * scheme->stage_error[j];
196: for (j = 1; j < r; j++) scheme->gamma[0] += v[0 * s + j] * scheme->gamma[j];
198: /* Assemble H
199: * % " PetscInt_FMT "etermine the error estimators phi
200: H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ...
201: [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]';
202: % Paper has formula above without the 0, but that term must be left
203: % out to satisfy the conditions they propose and to make the
204: % example schemes work
205: e.H = H;
206: e.phi = (H \ [1 0 0;1 1 0;0 0 -1])';
207: e.psi = -e.phi*C;
208: * */
209: for (j = 0; j < s; j++) {
210: H[0 + j * 3] = CPowF(c[j], p);
211: H[1 + j * 3] = CPowF(c[j], p + 1);
212: H[2 + j * 3] = scheme->stage_error[j];
213: for (k = 1; k < r; k++) {
214: H[0 + j * 3] += CPowF(c[j], k - 1) * scheme->alpha[k];
215: H[1 + j * 3] += CPowF(c[j], k - 1) * scheme->beta[k];
216: H[2 + j * 3] -= CPowF(c[j], k - 1) * scheme->gamma[k];
217: }
218: }
219: bmat[0 + 0 * ss] = 1.;
220: bmat[0 + 1 * ss] = 0.;
221: bmat[0 + 2 * ss] = 0.;
222: bmat[1 + 0 * ss] = 1.;
223: bmat[1 + 1 * ss] = 1.;
224: bmat[1 + 2 * ss] = 0.;
225: bmat[2 + 0 * ss] = 0.;
226: bmat[2 + 1 * ss] = 0.;
227: bmat[2 + 2 * ss] = -1.;
228: m = 3;
229: PetscBLASIntCast(s, &n);
230: PetscBLASIntCast(ss, &ldb);
231: rcond = 1e-12;
232: #if defined(PETSC_USE_COMPLEX)
233: /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
234: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, workreal, &info));
235: #else
236: /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
237: PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, &info));
238: #endif
242: for (j = 0; j < 3; j++) {
243: for (k = 0; k < s; k++) scheme->phi[k + j * s] = bmat[k + j * ss];
244: }
246: /* the other part of the error estimator, psi in B,J,W 2007 */
247: scheme->psi[0 * r + 0] = 0.;
248: scheme->psi[1 * r + 0] = 0.;
249: scheme->psi[2 * r + 0] = 0.;
250: for (j = 1; j < r; j++) {
251: scheme->psi[0 * r + j] = 0.;
252: scheme->psi[1 * r + j] = 0.;
253: scheme->psi[2 * r + j] = 0.;
254: for (k = 0; k < s; k++) {
255: scheme->psi[0 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[0 * s + k];
256: scheme->psi[1 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[1 * s + k];
257: scheme->psi[2 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[2 * s + k];
258: }
259: }
260: PetscFree7(ImV, H, bmat, workscalar, workreal, sing, ipiv);
261: }
262: /* Check which properties are satisfied */
263: scheme->stiffly_accurate = PETSC_TRUE;
264: if (scheme->c[s - 1] != 1.) scheme->stiffly_accurate = PETSC_FALSE;
265: for (j = 0; j < s; j++)
266: if (a[(s - 1) * s + j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE;
267: for (j = 0; j < r; j++)
268: if (u[(s - 1) * r + j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE;
269: scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */
270: for (j = 0; j < s - 1; j++)
271: if (r > 1 && b[1 * s + j] != 0.) scheme->fsal = PETSC_FALSE;
272: if (b[1 * s + r - 1] != 1.) scheme->fsal = PETSC_FALSE;
273: for (j = 0; j < r; j++)
274: if (r > 1 && v[1 * r + j] != 0.) scheme->fsal = PETSC_FALSE;
276: *inscheme = scheme;
277: return 0;
278: }
280: static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc)
281: {
282: PetscFree5(sc->c, sc->a, sc->b, sc->u, sc->v);
283: PetscFree6(sc->alpha, sc->beta, sc->gamma, sc->phi, sc->psi, sc->stage_error);
284: PetscFree(sc);
285: return 0;
286: }
288: static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
289: {
290: PetscInt i;
292: for (i = 0; i < gl->nschemes; i++) {
293: if (gl->schemes[i]) TSGLLESchemeDestroy(gl->schemes[i]);
294: }
295: PetscFree(gl->schemes);
296: gl->nschemes = 0;
297: PetscMemzero(gl->type_name, sizeof(gl->type_name));
298: return 0;
299: }
301: static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer, PetscInt m, PetscInt n, const PetscScalar a[], const char name[])
302: {
303: PetscBool iascii;
304: PetscInt i, j;
306: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
307: if (iascii) {
308: PetscViewerASCIIPrintf(viewer, "%30s = [", name);
309: for (i = 0; i < m; i++) {
310: if (i) PetscViewerASCIIPrintf(viewer, "%30s [", "");
311: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
312: for (j = 0; j < n; j++) PetscViewerASCIIPrintf(viewer, " %12.8g", (double)PetscRealPart(a[i * n + j]));
313: PetscViewerASCIIPrintf(viewer, "]\n");
314: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
315: }
316: }
317: return 0;
318: }
320: static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer)
321: {
322: PetscBool iascii;
324: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
325: if (iascii) {
326: 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);
327: PetscViewerASCIIPushTab(viewer);
328: PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s, FSAL: %s\n", sc->stiffly_accurate ? "yes" : "no", sc->fsal ? "yes" : "no");
329: 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]));
330: TSGLLEViewTable_Private(viewer, 1, sc->s, sc->c, "Abscissas c");
331: if (view_details) {
332: TSGLLEViewTable_Private(viewer, sc->s, sc->s, sc->a, "A");
333: TSGLLEViewTable_Private(viewer, sc->r, sc->s, sc->b, "B");
334: TSGLLEViewTable_Private(viewer, sc->s, sc->r, sc->u, "U");
335: TSGLLEViewTable_Private(viewer, sc->r, sc->r, sc->v, "V");
337: TSGLLEViewTable_Private(viewer, 3, sc->s, sc->phi, "Error estimate phi");
338: TSGLLEViewTable_Private(viewer, 3, sc->r, sc->psi, "Error estimate psi");
339: TSGLLEViewTable_Private(viewer, 1, sc->r, sc->alpha, "Modify alpha");
340: TSGLLEViewTable_Private(viewer, 1, sc->r, sc->beta, "Modify beta");
341: TSGLLEViewTable_Private(viewer, 1, sc->r, sc->gamma, "Modify gamma");
342: TSGLLEViewTable_Private(viewer, 1, sc->s, sc->stage_error, "Stage error xi");
343: }
344: PetscViewerASCIIPopTab(viewer);
345: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported", ((PetscObject)viewer)->type_name);
346: return 0;
347: }
349: static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[])
350: {
351: PetscInt i;
354: /* build error vectors*/
355: for (i = 0; i < 3; i++) {
356: PetscScalar phih[64];
357: PetscInt j;
358: for (j = 0; j < sc->s; j++) phih[j] = sc->phi[i * sc->s + j] * h;
359: VecZeroEntries(hm[i]);
360: VecMAXPY(hm[i], sc->s, phih, Ydot);
361: VecMAXPY(hm[i], sc->r, &sc->psi[i * sc->r], Xold);
362: }
363: return 0;
364: }
366: static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
367: {
368: PetscScalar brow[32], vrow[32];
369: PetscInt i, j, r, s;
371: /* Build the new solution from (X,Ydot) */
372: r = sc->r;
373: s = sc->s;
374: for (i = 0; i < r; i++) {
375: VecZeroEntries(X[i]);
376: for (j = 0; j < s; j++) brow[j] = h * sc->b[i * s + j];
377: VecMAXPY(X[i], s, brow, Ydot);
378: for (j = 0; j < r; j++) vrow[j] = sc->v[i * r + j];
379: VecMAXPY(X[i], r, vrow, Xold);
380: }
381: return 0;
382: }
384: static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
385: {
386: PetscScalar brow[32], vrow[32];
387: PetscReal ratio;
388: PetscInt i, j, p, r, s;
390: /* Build the new solution from (X,Ydot) */
391: p = sc->p;
392: r = sc->r;
393: s = sc->s;
394: ratio = next_h / h;
395: for (i = 0; i < r; i++) {
396: VecZeroEntries(X[i]);
397: for (j = 0; j < s; j++) {
398: 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]));
399: }
400: VecMAXPY(X[i], s, brow, Ydot);
401: for (j = 0; j < r; j++) {
402: 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]));
403: }
404: VecMAXPY(X[i], r, vrow, Xold);
405: }
406: if (r < next_sc->r) {
408: VecZeroEntries(X[r]);
409: for (j = 0; j < s; j++) brow[j] = h * PetscPowRealInt(ratio, p + 1) * sc->phi[0 * s + j];
410: VecMAXPY(X[r], s, brow, Ydot);
411: for (j = 0; j < r; j++) vrow[j] = PetscPowRealInt(ratio, p + 1) * sc->psi[0 * r + j];
412: VecMAXPY(X[r], r, vrow, Xold);
413: }
414: return 0;
415: }
417: static PetscErrorCode TSGLLECreate_IRKS(TS ts)
418: {
419: TS_GLLE *gl = (TS_GLLE *)ts->data;
421: gl->Destroy = TSGLLEDestroy_Default;
422: gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default;
423: gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
424: PetscMalloc1(10, &gl->schemes);
425: gl->nschemes = 0;
427: {
428: /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3
429: * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE.
430: * irks(0.3,0,[.3,1],[1],1)
431: * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2)
432: * but doing so would sacrifice the error estimator.
433: */
434: const PetscScalar c[2] = {3. / 10., 1.};
435: const PetscScalar a[2][2] = {
436: {3. / 10., 0 },
437: {7. / 10., 3. / 10.}
438: };
439: const PetscScalar b[2][2] = {
440: {7. / 10., 3. / 10.},
441: {0, 1 }
442: };
443: const PetscScalar u[2][2] = {
444: {1, 0},
445: {1, 0}
446: };
447: const PetscScalar v[2][2] = {
448: {1, 0},
449: {0, 0}
450: };
451: TSGLLESchemeCreate(1, 1, 2, 2, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
452: }
454: {
455: /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */
456: /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */
457: const PetscScalar c[3] = {1. / 3., 2. / 3., 1};
458: const PetscScalar a[3][3] = {
459: {4. / 9., 0, 0 },
460: {1.03750643704090e+00, 4. / 9., 0 },
461: {7.67024779410304e-01, -3.81140216918943e-01, 4. / 9.}
462: };
463: const PetscScalar b[3][3] = {
464: {0.767024779410304, -0.381140216918943, 4. / 9. },
465: {0.000000000000000, 0.000000000000000, 1.000000000000000},
466: {-2.075048385225385, 0.621728385225383, 1.277197204924873}
467: };
468: const PetscScalar u[3][3] = {
469: {1.0000000000000000, -0.1111111111111109, -0.0925925925925922},
470: {1.0000000000000000, -0.8152842148186744, -0.4199095530877056},
471: {1.0000000000000000, 0.1696709930641948, 0.0539741070314165 }
472: };
473: const PetscScalar v[3][3] = {
474: {1.0000000000000000, 0.1696709930641948, 0.0539741070314165},
475: {0.000000000000000, 0.000000000000000, 0.000000000000000 },
476: {0.000000000000000, 0.176122795075129, 0.000000000000000 }
477: };
478: TSGLLESchemeCreate(2, 2, 3, 3, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
479: }
480: {
481: /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */
482: const PetscScalar c[4] = {0.25, 0.5, 0.75, 1.0};
483: const PetscScalar a[4][4] = {
484: {9. / 40., 0, 0, 0 },
485: {2.11286958887701e-01, 9. / 40., 0, 0 },
486: {9.46338294287584e-01, -3.42942861246094e-01, 9. / 40., 0 },
487: {0.521490453970721, -0.662474225622980, 0.490476425459734, 9. / 40.}
488: };
489: const PetscScalar b[4][4] = {
490: {0.521490453970721, -0.662474225622980, 0.490476425459734, 9. / 40. },
491: {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000},
492: {-0.084677029310348, 1.390757514776085, -1.568157386206001, 2.023192696767826},
493: {0.465383797936408, 1.478273530625148, -1.930836081010182, 1.644872111193354}
494: };
495: const PetscScalar u[4][4] = {
496: {1.00000000000000000, 0.02500000000001035, -0.02499999999999053, -0.00442708333332865},
497: {1.00000000000000000, 0.06371304111232945, -0.04032173972189845, -0.01389438413189452},
498: {1.00000000000000000, -0.07839543304147778, 0.04738685705116663, 0.02032603595928376 },
499: {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034}
500: };
501: const PetscScalar v[4][4] = {
502: {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034},
503: {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000 },
504: {0.000000000000000, -1.761115796027561, -0.521284157173780, 0.258249384305463 },
505: {0.000000000000000, -1.657693358744728, -1.052227765232394, 0.521284157173780 }
506: };
507: TSGLLESchemeCreate(3, 3, 4, 4, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
508: }
509: {
510: /* p=q=4, r=s=5:
511: irks(3/11,0,[1:5]/5, [0.1715 -0.1238 0.6617],...
512: [ -0.0812 0.4079 1.0000
513: 1.0000 0 0
514: 0.8270 1.0000 0])
515: */
516: const PetscScalar c[5] = {0.2, 0.4, 0.6, 0.8, 1.0};
517: const PetscScalar a[5][5] = {
518: {2.72727272727352e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00},
519: {-1.03980153733431e-01, 2.72727272727405e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00},
520: {-1.58615400341492e+00, 7.44168951881122e-01, 2.72727272727309e-01, 0.00000000000000e+00, 0.00000000000000e+00},
521: {-8.73658042865628e-01, 5.37884671894595e-01, -1.63298538799523e-01, 2.72727272726996e-01, 0.00000000000000e+00},
522: {2.95489397443992e-01, -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01}
523: };
524: const PetscScalar b[5][5] = {
525: {2.95489397443992e-01, -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01},
526: {0.00000000000000e+00, 1.11022302462516e-16, -2.22044604925031e-16, 0.00000000000000e+00, 1.00000000000000e+00},
527: {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00, 6.32331093108427e-01},
528: {8.35690179937017e+00, -2.26640927349732e+00, 6.86647884973826e+00, -5.22595158025740e+00, 4.50893068837431e+00},
529: {1.27656267027479e+01, 2.80882153840821e+00, 8.91173096522890e+00, -1.07936444078906e+01, 4.82534148988854e+00}
530: };
531: const PetscScalar u[5][5] = {
532: {1.00000000000000e+00, -7.27272727273551e-02, -3.45454545454419e-02, -4.12121212119565e-03, -2.96969696964014e-04},
533: {1.00000000000000e+00, 2.31252881006154e-01, -8.29487834416481e-03, -9.07191207681020e-03, -1.70378403743473e-03},
534: {1.00000000000000e+00, 1.16925777880663e+00, 3.59268562942635e-02, -4.09013451730615e-02, -1.02411119670164e-02},
535: {1.00000000000000e+00, 1.02634463704356e+00, 1.59375044913405e-01, 1.89673015035370e-03, -4.89987231897569e-03},
536: {1.00000000000000e+00, 1.27746320298021e+00, 2.37186008132728e-01, -8.28694373940065e-02, -5.34396510196430e-02}
537: };
538: const PetscScalar v[5][5] = {
539: {1.00000000000000e+00, 1.27746320298021e+00, 2.37186008132728e-01, -8.28694373940065e-02, -5.34396510196430e-02},
540: {0.00000000000000e+00, -1.77635683940025e-15, -1.99840144432528e-15, -9.99200722162641e-16, -3.33066907387547e-16},
541: {0.00000000000000e+00, 4.37280081906924e+00, 5.49221645016377e-02, -8.88913177394943e-02, 1.12879077989154e-01 },
542: {0.00000000000000e+00, -1.22399504837280e+01, -5.21287338448645e+00, -8.03952325565291e-01, 4.60298678047147e-01 },
543: {0.00000000000000e+00, -1.85178762883829e+01, -5.21411849862624e+00, -1.04283436528809e+00, 7.49030161063651e-01 }
544: };
545: TSGLLESchemeCreate(4, 4, 5, 5, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
546: }
547: {
548: /* p=q=5, r=s=6;
549: irks(1/3,0,[1:6]/6,...
550: [-0.0489 0.4228 -0.8814 0.9021],...
551: [-0.3474 -0.6617 0.6294 0.2129
552: 0.0044 -0.4256 -0.1427 -0.8936
553: -0.8267 0.4821 0.1371 -0.2557
554: -0.4426 -0.3855 -0.7514 0.3014])
555: */
556: const PetscScalar c[6] = {1. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1.};
557: const PetscScalar a[6][6] = {
558: {3.33333333333940e-01, 0, 0, 0, 0, 0 },
559: {-8.64423857333350e-02, 3.33333333332888e-01, 0, 0, 0, 0 },
560: {-2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01, 0, 0, 0 },
561: {-4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01, 0, 0 },
562: {-6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01, -4.48352364517632e-01, 3.33333333328483e-01, 0 },
563: {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}
564: };
565: const PetscScalar b[6][6] = {
566: {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01 },
567: {-8.88178419700125e-16, 4.44089209850063e-16, -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00, 1.00000000000001e+00 },
568: {-2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01, 2.56943874812797e+01, -3.06702268304488e+01, 6.68067773510103e+00 },
569: {5.47971245256474e+01, 6.80366875868284e+01, -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01, -1.17819043489036e+01},
570: {-2.33332114788869e+02, 6.12942539462634e+01, -4.91850135865944e+01, 1.82716844135480e+02, -1.29788173979395e+02, 3.09968095651099e+01 },
571: {-1.72049132343751e+02, 8.60194713593999e+00, 7.98154219170200e-01, 1.50371386053218e+02, -1.18515423962066e+02, 2.50898277784663e+01 }
572: };
573: const PetscScalar u[6][6] = {
574: {1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06},
575: {1.00000000000000e+00, 8.64423857327162e-02, -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04},
576: {1.00000000000000e+00, 4.57135912953434e+00, 1.06514719719137e+00, 1.33517564218007e-01, 1.11365952968659e-02, 6.12382756769504e-04 },
577: {1.00000000000000e+00, 9.23391519753404e+00, 2.22431212392095e+00, 2.91823807741891e-01, 2.52058456411084e-02, 1.22800542949647e-03 },
578: {1.00000000000000e+00, 1.48175480533865e+01, 3.73439117461835e+00, 5.14648336541804e-01, 4.76430038853402e-02, 2.56798515502156e-03 },
579: {1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03}
580: };
581: const PetscScalar v[6][6] = {
582: {1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03},
583: {0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16},
584: {0.00000000000000e+00, 1.22250171233141e+01, -1.77150760606169e+00, 3.54516769879390e-01, 6.22298845883398e-01, 2.31647447450276e-01 },
585: {0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01, 6.55727990241799e-02, 1.63175368287079e-01 },
586: {0.00000000000000e+00, 1.37297394708005e+02, -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01, 9.16629423682464e-01 },
587: {0.00000000000000e+00, 1.05703241119022e+02, -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00 }
588: };
589: TSGLLESchemeCreate(5, 5, 6, 6, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
590: }
591: return 0;
592: }
594: /*@C
595: TSGLLESetType - sets the class of general linear method, `TSGLLE` to use for time-stepping
597: Collective
599: Input Parameters:
600: + ts - the `TS` context
601: - type - a method
603: Options Database Key:
604: . -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks)
606: Level: intermediate
608: Notes:
609: See "petsc/include/petscts.h" for available methods (for instance)
610: . TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems)
612: Normally, it is best to use the `TSSetFromOptions()` command and
613: then set the `TSGLLE` type from the options database rather than by using
614: this routine. Using the options database provides the user with
615: maximum flexibility in evaluating the many different solvers.
616: The `TSGLLESetType()` routine is provided for those situations where it
617: is necessary to set the timestepping solver independently of the
618: command line or options database. This might be the case, for example,
619: when the choice of solver changes during the execution of the
620: program, and the user's application is taking responsibility for
621: choosing the appropriate method.
623: .seealso: [](chapter_ts), `TS`, `TSGLLEType`, `TSGLLE`
624: @*/
625: PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type)
626: {
629: PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type));
630: return 0;
631: }
633: /*@C
634: TSGLLESetAcceptType - sets the acceptance test for `TSGLLE`
636: Time integrators that need to control error must have the option to reject a time step based on local error
637: estimates. This function allows different schemes to be set.
639: Logically Collective
641: Input Parameters:
642: + ts - the `TS` context
643: - type - the type
645: Options Database Key:
646: . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step
648: Level: intermediate
650: .seealso: [](chapter_ts), `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt`
651: @*/
652: PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type)
653: {
656: PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type));
657: return 0;
658: }
660: /*@C
661: TSGLLEGetAdapt - gets the `TSGLLEAdapt` object from the `TS`
663: Not Collective
665: Input Parameter:
666: . ts - the `TS` context
668: Output Parameter:
669: . adapt - the `TSGLLEAdapt` context
671: Level: advanced
673: Note:
674: This allows the user set options on the `TSGLLEAdapt` object. Usually it is better to do this using the options
675: database, so this function is rarely needed.
677: .seealso: [](chapter_ts), `TS`, `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegister()`
678: @*/
679: PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt)
680: {
683: PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt));
684: return 0;
685: }
687: static PetscErrorCode TSGLLEAccept_Always(TS ts, PetscReal tleft, PetscReal h, const PetscReal enorms[], PetscBool *accept)
688: {
689: *accept = PETSC_TRUE;
690: return 0;
691: }
693: static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
694: {
695: TS_GLLE *gl = (TS_GLLE *)ts->data;
696: PetscScalar *x, *w;
697: PetscInt n, i;
699: VecGetArray(gl->X[0], &x);
700: VecGetArray(gl->W, &w);
701: VecGetLocalSize(gl->W, &n);
702: for (i = 0; i < n; i++) w[i] = 1. / (gl->wrms_atol + gl->wrms_rtol * PetscAbsScalar(x[i]));
703: VecRestoreArray(gl->X[0], &x);
704: VecRestoreArray(gl->W, &w);
705: return 0;
706: }
708: static PetscErrorCode TSGLLEVecNormWRMS(TS ts, Vec X, PetscReal *nrm)
709: {
710: TS_GLLE *gl = (TS_GLLE *)ts->data;
711: PetscScalar *x, *w;
712: PetscReal sum = 0.0, gsum;
713: PetscInt n, N, i;
715: VecGetArray(X, &x);
716: VecGetArray(gl->W, &w);
717: VecGetLocalSize(gl->W, &n);
718: for (i = 0; i < n; i++) sum += PetscAbsScalar(PetscSqr(x[i] * w[i]));
719: VecRestoreArray(X, &x);
720: VecRestoreArray(gl->W, &w);
721: MPIU_Allreduce(&sum, &gsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));
722: VecGetSize(gl->W, &N);
723: *nrm = PetscSqrtReal(gsum / (1. * N));
724: return 0;
725: }
727: static PetscErrorCode TSGLLESetType_GLLE(TS ts, TSGLLEType type)
728: {
729: PetscBool same;
730: TS_GLLE *gl = (TS_GLLE *)ts->data;
731: PetscErrorCode (*r)(TS);
733: if (gl->type_name[0]) {
734: PetscStrcmp(gl->type_name, type, &same);
735: if (same) return 0;
736: (*gl->Destroy)(gl);
737: }
739: PetscFunctionListFind(TSGLLEList, type, &r);
741: (*r)(ts);
742: PetscStrcpy(gl->type_name, type);
743: return 0;
744: }
746: static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type)
747: {
748: TSGLLEAcceptFunction r;
749: TS_GLLE *gl = (TS_GLLE *)ts->data;
751: PetscFunctionListFind(TSGLLEAcceptList, type, &r);
753: gl->Accept = r;
754: PetscStrncpy(gl->accept_name, type, sizeof(gl->accept_name));
755: return 0;
756: }
758: static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt)
759: {
760: TS_GLLE *gl = (TS_GLLE *)ts->data;
762: if (!gl->adapt) {
763: TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts), &gl->adapt);
764: PetscObjectIncrementTabLevel((PetscObject)gl->adapt, (PetscObject)ts, 1);
765: }
766: *adapt = gl->adapt;
767: return 0;
768: }
770: static PetscErrorCode TSGLLEChooseNextScheme(TS ts, PetscReal h, const PetscReal hmnorm[], PetscInt *next_scheme, PetscReal *next_h, PetscBool *finish)
771: {
772: TS_GLLE *gl = (TS_GLLE *)ts->data;
773: PetscInt i, n, cur_p, cur, next_sc, candidates[64], orders[64];
774: PetscReal errors[64], costs[64], tleft;
776: cur = -1;
777: cur_p = gl->schemes[gl->current_scheme]->p;
778: tleft = ts->max_time - (ts->ptime + ts->time_step);
779: for (i = 0, n = 0; i < gl->nschemes; i++) {
780: TSGLLEScheme sc = gl->schemes[i];
781: if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
782: if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[0];
783: else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[1];
784: else if (sc->p == cur_p + 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * (hmnorm[2] + hmnorm[3]);
785: else continue;
786: candidates[n] = i;
787: orders[n] = PetscMin(sc->p, sc->q); /* order of global truncation error */
788: costs[n] = sc->s; /* estimate the cost as the number of stages */
789: if (i == gl->current_scheme) cur = n;
790: n++;
791: }
793: TSGLLEAdaptChoose(gl->adapt, n, orders, errors, costs, cur, h, tleft, &next_sc, next_h, finish);
794: *next_scheme = candidates[next_sc];
795: 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,
796: gl->schemes[*next_scheme]->r, gl->schemes[*next_scheme]->s, (double)*next_h, PetscBools[*finish]));
797: return 0;
798: }
800: static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s)
801: {
802: TS_GLLE *gl = (TS_GLLE *)ts->data;
804: *max_r = gl->schemes[gl->nschemes - 1]->r;
805: *max_s = gl->schemes[gl->nschemes - 1]->s;
806: return 0;
807: }
809: static PetscErrorCode TSSolve_GLLE(TS ts)
810: {
811: TS_GLLE *gl = (TS_GLLE *)ts->data;
812: PetscInt i, k, its, lits, max_r, max_s;
813: PetscBool final_step, finish;
814: SNESConvergedReason snesreason;
816: TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);
818: TSGLLEGetMaxSizes(ts, &max_r, &max_s);
819: VecCopy(ts->vec_sol, gl->X[0]);
820: for (i = 1; i < max_r; i++) VecZeroEntries(gl->X[i]);
821: TSGLLEUpdateWRMS(ts);
823: if (0) {
824: /* Find consistent initial data for DAE */
825: gl->stage_time = ts->ptime + ts->time_step;
826: gl->scoeff = 1.;
827: gl->stage = 0;
829: VecCopy(ts->vec_sol, gl->Z);
830: VecCopy(ts->vec_sol, gl->Y);
831: SNESSolve(ts->snes, NULL, gl->Y);
832: SNESGetIterationNumber(ts->snes, &its);
833: SNESGetLinearSolveIterations(ts->snes, &lits);
834: SNESGetConvergedReason(ts->snes, &snesreason);
836: ts->snes_its += its;
837: ts->ksp_its += lits;
838: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
839: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
840: PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures);
841: return 0;
842: }
843: }
847: for (k = 0, final_step = PETSC_FALSE, finish = PETSC_FALSE; k < ts->max_steps && !finish; k++) {
848: PetscInt j, r, s, next_scheme = 0, rejections;
849: PetscReal h, hmnorm[4], enorm[3], next_h;
850: PetscBool accept;
851: const PetscScalar *c, *a, *u;
852: Vec *X, *Ydot, Y;
853: TSGLLEScheme scheme = gl->schemes[gl->current_scheme];
855: r = scheme->r;
856: s = scheme->s;
857: c = scheme->c;
858: a = scheme->a;
859: u = scheme->u;
860: h = ts->time_step;
861: X = gl->X;
862: Ydot = gl->Ydot;
863: Y = gl->Y;
865: if (ts->ptime > ts->max_time) break;
867: /*
868: We only call PreStep at the start of each STEP, not each STAGE. This is because it is
869: possible to fail (have to restart a step) after multiple stages.
870: */
871: TSPreStep(ts);
873: rejections = 0;
874: while (1) {
875: for (i = 0; i < s; i++) {
876: PetscScalar shift;
877: gl->scoeff = 1. / PetscRealPart(a[i * s + i]);
878: shift = gl->scoeff / ts->time_step;
879: gl->stage = i;
880: gl->stage_time = ts->ptime + PetscRealPart(c[i]) * h;
882: /*
883: * Stage equation: Y = h A Y' + U X
884: * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
885: * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
886: * Then y'_i = z + 1/(h a_ii) y_i
887: */
888: VecZeroEntries(gl->Z);
889: for (j = 0; j < r; j++) VecAXPY(gl->Z, -shift * u[i * r + j], X[j]);
890: for (j = 0; j < i; j++) VecAXPY(gl->Z, -shift * h * a[i * s + j], Ydot[j]);
891: /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */
893: /* Compute an estimate of Y to start Newton iteration */
894: if (gl->extrapolate) {
895: if (i == 0) {
896: /* Linear extrapolation on the first stage */
897: VecWAXPY(Y, c[i] * h, X[1], X[0]);
898: } else {
899: /* Linear extrapolation from the last stage */
900: VecAXPY(Y, (c[i] - c[i - 1]) * h, Ydot[i - 1]);
901: }
902: } else if (i == 0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
903: VecCopy(X[0], Y);
904: }
906: /* Solve this stage (Ydot[i] is computed during function evaluation) */
907: SNESSolve(ts->snes, NULL, Y);
908: SNESGetIterationNumber(ts->snes, &its);
909: SNESGetLinearSolveIterations(ts->snes, &lits);
910: SNESGetConvergedReason(ts->snes, &snesreason);
911: ts->snes_its += its;
912: ts->ksp_its += lits;
913: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
914: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
915: PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures);
916: return 0;
917: }
918: }
920: gl->stage_time = ts->ptime + ts->time_step;
922: (*gl->EstimateHigherMoments)(scheme, h, Ydot, gl->X, gl->himom);
923: /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */
924: for (i = 0; i < 3; i++) TSGLLEVecNormWRMS(ts, gl->himom[i], &hmnorm[i + 1]);
925: enorm[0] = PetscRealPart(scheme->alpha[0]) * hmnorm[1];
926: enorm[1] = PetscRealPart(scheme->beta[0]) * hmnorm[2];
927: enorm[2] = PetscRealPart(scheme->gamma[0]) * hmnorm[3];
928: (*gl->Accept)(ts, ts->max_time - gl->stage_time, h, enorm, &accept);
929: if (accept) goto accepted;
930: rejections++;
931: PetscInfo(ts, "Step %" PetscInt_FMT " (t=%g) not accepted, rejections=%" PetscInt_FMT "\n", k, (double)gl->stage_time, rejections);
932: if (rejections > gl->max_step_rejections) break;
933: /*
934: There are lots of reasons why a step might be rejected, including solvers not converging and other factors that
935: TSGLLEChooseNextScheme does not support. Additionally, the error estimates may be very screwed up, so I'm not
936: convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor
937: (the adaptor interface probably has to change). Here we make an arbitrary and naive choice. This assumes that
938: steps were written in Nordsieck form. The "correct" method would be to re-complete the previous time step with
939: the correct "next" step size. It is unclear to me whether the present ad-hoc method of rescaling X is stable.
940: */
941: h *= 0.5;
942: for (i = 1; i < scheme->r; i++) VecScale(X[i], PetscPowRealInt(0.5, i));
943: }
944: 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);
946: accepted:
947: /* This term is not error, but it *would* be the leading term for a lower order method */
948: TSGLLEVecNormWRMS(ts, gl->X[scheme->r - 1], &hmnorm[0]);
949: /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */
951: 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]);
952: if (!final_step) {
953: TSGLLEChooseNextScheme(ts, h, hmnorm, &next_scheme, &next_h, &final_step);
954: } else {
955: /* Dummy values to complete the current step in a consistent manner */
956: next_scheme = gl->current_scheme;
957: next_h = h;
958: finish = PETSC_TRUE;
959: }
961: X = gl->Xold;
962: gl->Xold = gl->X;
963: gl->X = X;
964: (*gl->CompleteStep)(scheme, h, gl->schemes[next_scheme], next_h, Ydot, gl->Xold, gl->X);
966: TSGLLEUpdateWRMS(ts);
968: /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */
969: VecCopy(gl->X[0], ts->vec_sol);
970: ts->ptime += h;
971: ts->steps++;
973: TSPostEvaluate(ts);
974: TSPostStep(ts);
975: TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);
977: gl->current_scheme = next_scheme;
978: ts->time_step = next_h;
979: }
980: return 0;
981: }
983: /*------------------------------------------------------------*/
985: static PetscErrorCode TSReset_GLLE(TS ts)
986: {
987: TS_GLLE *gl = (TS_GLLE *)ts->data;
988: PetscInt max_r, max_s;
990: if (gl->setupcalled) {
991: TSGLLEGetMaxSizes(ts, &max_r, &max_s);
992: VecDestroyVecs(max_r, &gl->Xold);
993: VecDestroyVecs(max_r, &gl->X);
994: VecDestroyVecs(max_s, &gl->Ydot);
995: VecDestroyVecs(3, &gl->himom);
996: VecDestroy(&gl->W);
997: VecDestroy(&gl->Y);
998: VecDestroy(&gl->Z);
999: }
1000: gl->setupcalled = PETSC_FALSE;
1001: return 0;
1002: }
1004: static PetscErrorCode TSDestroy_GLLE(TS ts)
1005: {
1006: TS_GLLE *gl = (TS_GLLE *)ts->data;
1008: TSReset_GLLE(ts);
1009: if (ts->dm) {
1010: DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts);
1011: DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts);
1012: }
1013: if (gl->adapt) TSGLLEAdaptDestroy(&gl->adapt);
1014: if (gl->Destroy) (*gl->Destroy)(gl);
1015: PetscFree(ts->data);
1016: PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL);
1017: PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL);
1018: PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL);
1019: return 0;
1020: }
1022: /*
1023: This defines the nonlinear equation that is to be solved with SNES
1024: g(x) = f(t,x,z+shift*x) = 0
1025: */
1026: static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts)
1027: {
1028: TS_GLLE *gl = (TS_GLLE *)ts->data;
1029: Vec Z, Ydot;
1030: DM dm, dmsave;
1032: SNESGetDM(snes, &dm);
1033: TSGLLEGetVecs(ts, dm, &Z, &Ydot);
1034: VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z);
1035: dmsave = ts->dm;
1036: ts->dm = dm;
1037: TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE);
1038: ts->dm = dmsave;
1039: TSGLLERestoreVecs(ts, dm, &Z, &Ydot);
1040: return 0;
1041: }
1043: static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts)
1044: {
1045: TS_GLLE *gl = (TS_GLLE *)ts->data;
1046: Vec Z, Ydot;
1047: DM dm, dmsave;
1049: SNESGetDM(snes, &dm);
1050: TSGLLEGetVecs(ts, dm, &Z, &Ydot);
1051: dmsave = ts->dm;
1052: ts->dm = dm;
1053: /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
1054: TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE);
1055: ts->dm = dmsave;
1056: TSGLLERestoreVecs(ts, dm, &Z, &Ydot);
1057: return 0;
1058: }
1060: static PetscErrorCode TSSetUp_GLLE(TS ts)
1061: {
1062: TS_GLLE *gl = (TS_GLLE *)ts->data;
1063: PetscInt max_r, max_s;
1064: DM dm;
1066: gl->setupcalled = PETSC_TRUE;
1067: TSGLLEGetMaxSizes(ts, &max_r, &max_s);
1068: VecDuplicateVecs(ts->vec_sol, max_r, &gl->X);
1069: VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold);
1070: VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot);
1071: VecDuplicateVecs(ts->vec_sol, 3, &gl->himom);
1072: VecDuplicate(ts->vec_sol, &gl->W);
1073: VecDuplicate(ts->vec_sol, &gl->Y);
1074: VecDuplicate(ts->vec_sol, &gl->Z);
1076: /* Default acceptance tests and adaptivity */
1077: if (!gl->Accept) TSGLLESetAcceptType(ts, TSGLLEACCEPT_ALWAYS);
1078: if (!gl->adapt) TSGLLEGetAdapt(ts, &gl->adapt);
1080: if (gl->current_scheme < 0) {
1081: PetscInt i;
1082: for (i = 0;; i++) {
1083: if (gl->schemes[i]->p == gl->start_order) break;
1085: }
1086: gl->current_scheme = i;
1087: }
1088: TSGetDM(ts, &dm);
1089: DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts);
1090: DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts);
1091: return 0;
1092: }
1093: /*------------------------------------------------------------*/
1095: static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems *PetscOptionsObject)
1096: {
1097: TS_GLLE *gl = (TS_GLLE *)ts->data;
1098: char tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify";
1100: PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options");
1101: {
1102: PetscBool flg;
1103: PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg);
1104: if (flg || !gl->type_name[0]) TSGLLESetType(ts, tname);
1105: PetscOptionsInt("-ts_gl_max_step_rejections", "Maximum number of times to attempt a step", "None", gl->max_step_rejections, &gl->max_step_rejections, NULL);
1106: PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL);
1107: PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL);
1108: PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL);
1109: PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL);
1110: PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL);
1111: PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL);
1112: PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL);
1113: PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg);
1114: if (flg) {
1115: PetscBool match1, match2;
1116: PetscStrcmp(completef, "rescale", &match1);
1117: PetscStrcmp(completef, "rescale-and-modify", &match2);
1118: if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale;
1119: else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
1120: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef);
1121: }
1122: {
1123: char type[256] = TSGLLEACCEPT_ALWAYS;
1124: 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);
1125: if (flg || !gl->accept_name[0]) TSGLLESetAcceptType(ts, type);
1126: }
1127: {
1128: TSGLLEAdapt adapt;
1129: TSGLLEGetAdapt(ts, &adapt);
1130: TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject);
1131: }
1132: }
1133: PetscOptionsHeadEnd();
1134: return 0;
1135: }
1137: static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer)
1138: {
1139: TS_GLLE *gl = (TS_GLLE *)ts->data;
1140: PetscInt i;
1141: PetscBool iascii, details;
1143: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1144: if (iascii) {
1145: 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);
1146: PetscViewerASCIIPrintf(viewer, " Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction]);
1147: PetscViewerASCIIPrintf(viewer, " Extrapolation: %s\n", gl->extrapolate ? "yes" : "no");
1148: PetscViewerASCIIPrintf(viewer, " Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)");
1149: PetscViewerASCIIPushTab(viewer);
1150: TSGLLEAdaptView(gl->adapt, viewer);
1151: PetscViewerASCIIPopTab(viewer);
1152: PetscViewerASCIIPrintf(viewer, " type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)");
1153: PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes);
1154: details = PETSC_FALSE;
1155: PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL);
1156: PetscViewerASCIIPushTab(viewer);
1157: for (i = 0; i < gl->nschemes; i++) TSGLLESchemeView(gl->schemes[i], details, viewer);
1158: if (gl->View) (*gl->View)(gl, viewer);
1159: PetscViewerASCIIPopTab(viewer);
1160: }
1161: return 0;
1162: }
1164: /*@C
1165: TSGLLERegister - adds a `TSGLLE` implementation
1167: Not Collective
1169: Input Parameters:
1170: + name_scheme - name of user-defined general linear scheme
1171: - routine_create - routine to create method context
1173: Level: advanced
1175: Note:
1176: `TSGLLERegister()` may be called multiple times to add several user-defined families.
1178: Sample usage:
1179: .vb
1180: TSGLLERegister("my_scheme",MySchemeCreate);
1181: .ve
1183: Then, your scheme can be chosen with the procedural interface via
1184: $ TSGLLESetType(ts,"my_scheme")
1185: or at runtime via the option
1186: $ -ts_gl_type my_scheme
1188: .seealso: [](chapter_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`
1189: @*/
1190: PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS))
1191: {
1192: TSGLLEInitializePackage();
1193: PetscFunctionListAdd(&TSGLLEList, sname, function);
1194: return 0;
1195: }
1197: /*@C
1198: TSGLLEAcceptRegister - adds a `TSGLLE` acceptance scheme
1200: Not Collective
1202: Input Parameters:
1203: + name_scheme - name of user-defined acceptance scheme
1204: - routine_create - routine to create method context
1206: Level: advanced
1208: Note:
1209: `TSGLLEAcceptRegister()` may be called multiple times to add several user-defined families.
1211: Sample usage:
1212: .vb
1213: TSGLLEAcceptRegister("my_scheme",MySchemeCreate);
1214: .ve
1216: Then, your scheme can be chosen with the procedural interface via
1217: $ TSGLLESetAcceptType(ts,"my_scheme")
1218: or at runtime via the option
1219: $ -ts_gl_accept_type my_scheme
1221: .seealso: [](chapter_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLEAcceptFunction`
1222: @*/
1223: PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFunction function)
1224: {
1225: PetscFunctionListAdd(&TSGLLEAcceptList, sname, function);
1226: return 0;
1227: }
1229: /*@C
1230: TSGLLERegisterAll - Registers all of the general linear methods in `TSGLLE`
1232: Not Collective
1234: Level: advanced
1236: .seealso: [](chapter_ts), `TSGLLE`, `TSGLLERegisterDestroy()`
1237: @*/
1238: PetscErrorCode TSGLLERegisterAll(void)
1239: {
1240: if (TSGLLERegisterAllCalled) return 0;
1241: TSGLLERegisterAllCalled = PETSC_TRUE;
1243: TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS);
1244: TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always);
1245: return 0;
1246: }
1248: /*@C
1249: TSGLLEInitializePackage - This function initializes everything in the `TSGLLE` package. It is called
1250: from `TSInitializePackage()`.
1252: Level: developer
1254: .seealso: [](chapter_ts), `PetscInitialize()`, `TSInitializePackage()`, `TSGLLEFinalizePackage()`
1255: @*/
1256: PetscErrorCode TSGLLEInitializePackage(void)
1257: {
1258: if (TSGLLEPackageInitialized) return 0;
1259: TSGLLEPackageInitialized = PETSC_TRUE;
1260: TSGLLERegisterAll();
1261: PetscRegisterFinalize(TSGLLEFinalizePackage);
1262: return 0;
1263: }
1265: /*@C
1266: TSGLLEFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
1267: called from `PetscFinalize()`.
1269: Level: developer
1271: .seealso: [](chapter_ts), `PetscFinalize()`, `TSGLLEInitializePackage()`, `TSInitializePackage()`
1272: @*/
1273: PetscErrorCode TSGLLEFinalizePackage(void)
1274: {
1275: PetscFunctionListDestroy(&TSGLLEList);
1276: PetscFunctionListDestroy(&TSGLLEAcceptList);
1277: TSGLLEPackageInitialized = PETSC_FALSE;
1278: TSGLLERegisterAllCalled = PETSC_FALSE;
1279: return 0;
1280: }
1282: /* ------------------------------------------------------------ */
1283: /*MC
1284: TSGLLE - DAE solver using implicit General Linear methods
1286: These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental
1287: limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their
1288: applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF
1289: are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and
1290: reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes.
1291: All this is possible while preserving a singly diagonally implicit structure.
1293: Options Database Keys:
1294: + -ts_gl_type <type> - the class of general linear method (irks)
1295: . -ts_gl_rtol <tol> - relative error
1296: . -ts_gl_atol <tol> - absolute error
1297: . -ts_gl_min_order <p> - minimum order method to consider (default=1)
1298: . -ts_gl_max_order <p> - maximum order method to consider (default=3)
1299: . -ts_gl_start_order <p> - order of starting method (default=1)
1300: . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
1301: - -ts_adapt_type <method> - adaptive controller to use (none step both)
1303: Level: beginner
1305: Notes:
1306: This integrator can be applied to DAE.
1308: Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK).
1309: They are represented by the tableau
1311: .vb
1312: A | U
1313: -------
1314: B | V
1315: .ve
1317: combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular.
1318: A step of the general method reads
1320: .vb
1321: [ Y ] = [A U] [ Y' ]
1322: [X^k] = [B V] [X^{k-1}]
1323: .ve
1325: where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of
1326: the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by
1328: .vb
1329: X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
1330: .ve
1332: If A is lower triangular, we can solve the stages (Y,Y') sequentially
1334: .vb
1335: 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}
1336: .ve
1338: and then construct the pieces to carry to the next step
1340: .vb
1341: 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}
1342: .ve
1344: Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i
1345: in terms of y_i and known stuff (y_j for j<i and x_j for all j).
1347: Error estimation
1349: At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses
1350: Inherent Runge-Kutta Stability (`TSIRKS`). These methods have r=s, the number of items passed between steps is equal to
1351: the number of stages. The order and stage-order are one less than the number of stages. We use the error estimates
1352: in the 2007 paper which provide the following estimates
1354: .vb
1355: h^{p+1} X^{(p+1)} = phi_0^T Y' + [0 psi_0^T] Xold
1356: h^{p+2} X^{(p+2)} = phi_1^T Y' + [0 psi_1^T] Xold
1357: h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold
1358: .ve
1360: These estimates are accurate to O(h^{p+3}).
1362: Changing the step size
1364: We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper.
1366: References:
1367: + * - John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for
1368: ordinary differential equations, Journal of Complexity, Vol 23, 2007.
1369: - * - John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009.
1371: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType`
1372: M*/
1373: PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
1374: {
1375: TS_GLLE *gl;
1377: TSGLLEInitializePackage();
1379: PetscNew(&gl);
1380: ts->data = (void *)gl;
1382: ts->ops->reset = TSReset_GLLE;
1383: ts->ops->destroy = TSDestroy_GLLE;
1384: ts->ops->view = TSView_GLLE;
1385: ts->ops->setup = TSSetUp_GLLE;
1386: ts->ops->solve = TSSolve_GLLE;
1387: ts->ops->setfromoptions = TSSetFromOptions_GLLE;
1388: ts->ops->snesfunction = SNESTSFormFunction_GLLE;
1389: ts->ops->snesjacobian = SNESTSFormJacobian_GLLE;
1391: ts->usessnes = PETSC_TRUE;
1393: gl->max_step_rejections = 1;
1394: gl->min_order = 1;
1395: gl->max_order = 3;
1396: gl->start_order = 1;
1397: gl->current_scheme = -1;
1398: gl->extrapolate = PETSC_FALSE;
1400: gl->wrms_atol = 1e-8;
1401: gl->wrms_rtol = 1e-5;
1403: PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE);
1404: PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE);
1405: PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE);
1406: return 0;
1407: }