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