1 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/
2 #include <petscdm.h>
3 #include <petscblaslapack.h>
4
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;
10
11 /* This function is pure */
Factorial(PetscInt n)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 }
25
26 /* This function is pure */
CPowF(PetscScalar c,PetscInt p)27 static PetscScalar CPowF(PetscScalar c, PetscInt p)
28 {
29 return PetscPowRealInt(PetscRealPart(c), p) / Factorial(p);
30 }
31
TSGLLEGetVecs(TS ts,DM dm,Vec * Z,Vec * Ydotstage)32 static PetscErrorCode TSGLLEGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
33 {
34 TS_GLLE *gl = (TS_GLLE *)ts->data;
35
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 }
47
TSGLLERestoreVecs(TS ts,DM dm,Vec * Z,Vec * Ydotstage)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 }
59
DMCoarsenHook_TSGLLE(DM fine,DM coarse,PetscCtx ctx)60 static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, PetscCtx ctx)
61 {
62 PetscFunctionBegin;
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
DMRestrictHook_TSGLLE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,PetscCtx ctx)66 static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
67 {
68 TS ts = (TS)ctx;
69 Vec Ydot, Ydot_c;
70
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 }
80
DMSubDomainHook_TSGLLE(DM dm,DM subdm,PetscCtx ctx)81 static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, PetscCtx ctx)
82 {
83 PetscFunctionBegin;
84 PetscFunctionReturn(PETSC_SUCCESS);
85 }
86
DMSubDomainRestrictHook_TSGLLE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,PetscCtx ctx)87 static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
88 {
89 TS ts = (TS)ctx;
90 Vec Ydot, Ydot_s;
91
92 PetscFunctionBegin;
93 PetscCall(TSGLLEGetVecs(ts, dm, NULL, &Ydot));
94 PetscCall(TSGLLEGetVecs(ts, subdm, NULL, &Ydot_s));
95
96 PetscCall(VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));
97 PetscCall(VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));
98
99 PetscCall(TSGLLERestoreVecs(ts, dm, NULL, &Ydot));
100 PetscCall(TSGLLERestoreVecs(ts, subdm, NULL, &Ydot_s));
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
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)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;
108
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;
120
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];
127
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;
135
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));
138
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");
153
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");
162
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 }
171
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];
178
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");
187
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];
196
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];
203
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");
247
248 for (j = 0; j < 3; j++) {
249 for (k = 0; k < s; k++) scheme->phi[k + j * s] = bmat[k + j * ss];
250 }
251
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;
281
282 *inscheme = scheme;
283 PetscFunctionReturn(PETSC_SUCCESS);
284 }
285
TSGLLESchemeDestroy(TSGLLEScheme sc)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 }
294
TSGLLEDestroy_Default(TS_GLLE * gl)295 static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
296 {
297 PetscInt i;
298
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 }
308
TSGLLEViewTable_Private(PetscViewer viewer,PetscInt m,PetscInt n,const PetscScalar a[],const char name[])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;
313
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 }
328
TSGLLESchemeView(TSGLLEScheme sc,PetscBool view_details,PetscViewer viewer)329 static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer)
330 {
331 PetscBool isascii;
332
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"));
346
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 }
357
TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc,PetscReal h,Vec Ydot[],Vec Xold[],Vec hm[])358 static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[])
359 {
360 PetscInt i;
361
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 }
375
TSGLLECompleteStep_Rescale(TSGLLEScheme sc,PetscReal h,TSGLLEScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[])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;
380
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 }
394
TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc,PetscReal h,TSGLLEScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[])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;
400
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 }
428
TSGLLECreate_IRKS(TS ts)429 static PetscErrorCode TSGLLECreate_IRKS(TS ts)
430 {
431 TS_GLLE *gl = (TS_GLLE *)ts->data;
432
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;
439
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 }
466
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 }
606
607 /*@
608 TSGLLESetType - sets the class of general linear method, `TSGLLE` to use for time-stepping
609
610 Collective
611
612 Input Parameters:
613 + ts - the `TS` context
614 - type - a method
615
616 Options Database Key:
617 . -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks)
618
619 Level: intermediate
620
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)
624
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.
632
633 .seealso: [](ch_ts), `TS`, `TSGLLEType`, `TSGLLE`
634 @*/
TSGLLESetType(TS ts,TSGLLEType type)635 PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type)
636 {
637 PetscFunctionBegin;
638 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
639 PetscAssertPointer(type, 2);
640 PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type));
641 PetscFunctionReturn(PETSC_SUCCESS);
642 }
643
644 /*@C
645 TSGLLESetAcceptType - sets the acceptance test for `TSGLLE`
646
647 Logically Collective
648
649 Input Parameters:
650 + ts - the `TS` context
651 - type - the type
652
653 Options Database Key:
654 . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step
655
656 Level: intermediate
657
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.
661
662 .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt`
663 @*/
TSGLLESetAcceptType(TS ts,TSGLLEAcceptType type)664 PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type)
665 {
666 PetscFunctionBegin;
667 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
668 PetscAssertPointer(type, 2);
669 PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type));
670 PetscFunctionReturn(PETSC_SUCCESS);
671 }
672
673 /*@
674 TSGLLEGetAdapt - gets the `TSGLLEAdapt` object from the `TS`
675
676 Not Collective
677
678 Input Parameter:
679 . ts - the `TS` context
680
681 Output Parameter:
682 . adapt - the `TSGLLEAdapt` context
683
684 Level: advanced
685
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.
689
690 .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegister()`
691 @*/
TSGLLEGetAdapt(TS ts,TSGLLEAdapt * adapt)692 PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt)
693 {
694 PetscFunctionBegin;
695 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
696 PetscAssertPointer(adapt, 2);
697 PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt));
698 PetscFunctionReturn(PETSC_SUCCESS);
699 }
700
TSGLLEAccept_Always(TS ts,PetscReal tleft,PetscReal h,const PetscReal enorms[],PetscBool * accept)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 }
707
TSGLLEUpdateWRMS(TS ts)708 static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
709 {
710 TS_GLLE *gl = (TS_GLLE *)ts->data;
711 PetscScalar *x, *w;
712 PetscInt n, i;
713
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 }
723
TSGLLEVecNormWRMS(TS ts,Vec X,PetscReal * nrm)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;
730
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 }
743
TSGLLESetType_GLLE(TS ts,TSGLLEType type)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);
749
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 }
756
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 }
763
TSGLLESetAcceptType_GLLE(TS ts,TSGLLEAcceptType type)764 static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type)
765 {
766 TSGLLEAcceptFn *r;
767 TS_GLLE *gl = (TS_GLLE *)ts->data;
768
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 }
776
TSGLLEGetAdapt_GLLE(TS ts,TSGLLEAdapt * adapt)777 static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt)
778 {
779 TS_GLLE *gl = (TS_GLLE *)ts->data;
780
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 }
789
TSGLLEChooseNextScheme(TS ts,PetscReal h,const PetscReal hmnorm[],PetscInt * next_scheme,PetscReal * next_h,PetscBool * finish)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;
795
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 }
820
TSGLLEGetMaxSizes(TS ts,PetscInt * max_r,PetscInt * max_s)821 static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s)
822 {
823 TS_GLLE *gl = (TS_GLLE *)ts->data;
824
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 }
830
TSSolve_GLLE(TS ts)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;
837
838 PetscFunctionBegin;
839 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
840
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));
845
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;
851
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));
858
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 }
867
868 PetscCheck(gl->current_scheme >= 0, PETSC_COMM_SELF, PETSC_ERR_ORDER, "A starting scheme has not been provided");
869
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];
877
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;
887
888 if (ts->ptime > ts->max_time) break;
889
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));
895
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;
904
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 */
915
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 }
928
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 }
942
943 gl->stage_time = ts->ptime + ts->time_step;
944
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);
968
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 */
973
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 }
983
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));
988
989 PetscCall(TSGLLEUpdateWRMS(ts));
990
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++;
995
996 PetscCall(TSPostEvaluate(ts));
997 PetscCall(TSPostStep(ts));
998 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
999
1000 gl->current_scheme = next_scheme;
1001 ts->time_step = next_h;
1002 }
1003 PetscFunctionReturn(PETSC_SUCCESS);
1004 }
1005
TSReset_GLLE(TS ts)1006 static PetscErrorCode TSReset_GLLE(TS ts)
1007 {
1008 TS_GLLE *gl = (TS_GLLE *)ts->data;
1009 PetscInt max_r, max_s;
1010
1011 PetscFunctionBegin;
1012 if (gl->setupcalled) {
1013 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
1014 PetscCall(VecDestroyVecs(max_r, &gl->Xold));
1015 PetscCall(VecDestroyVecs(max_r, &gl->X));
1016 PetscCall(VecDestroyVecs(max_s, &gl->Ydot));
1017 PetscCall(VecDestroyVecs(3, &gl->himom));
1018 PetscCall(VecDestroy(&gl->W));
1019 PetscCall(VecDestroy(&gl->Y));
1020 PetscCall(VecDestroy(&gl->Z));
1021 }
1022 gl->setupcalled = PETSC_FALSE;
1023 PetscFunctionReturn(PETSC_SUCCESS);
1024 }
1025
TSDestroy_GLLE(TS ts)1026 static PetscErrorCode TSDestroy_GLLE(TS ts)
1027 {
1028 TS_GLLE *gl = (TS_GLLE *)ts->data;
1029
1030 PetscFunctionBegin;
1031 PetscCall(TSReset_GLLE(ts));
1032 if (ts->dm) {
1033 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
1034 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
1035 }
1036 if (gl->adapt) PetscCall(TSGLLEAdaptDestroy(&gl->adapt));
1037 if (gl->Destroy) PetscCall((*gl->Destroy)(gl));
1038 PetscCall(PetscFree(ts->data));
1039 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL));
1040 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL));
1041 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL));
1042 PetscFunctionReturn(PETSC_SUCCESS);
1043 }
1044
1045 /*
1046 This defines the nonlinear equation that is to be solved with SNES
1047 g(x) = f(t,x,z+shift*x) = 0
1048 */
SNESTSFormFunction_GLLE(SNES snes,Vec x,Vec f,TS ts)1049 static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts)
1050 {
1051 TS_GLLE *gl = (TS_GLLE *)ts->data;
1052 Vec Z, Ydot;
1053 DM dm, dmsave;
1054
1055 PetscFunctionBegin;
1056 PetscCall(SNESGetDM(snes, &dm));
1057 PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
1058 PetscCall(VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z));
1059 dmsave = ts->dm;
1060 ts->dm = dm;
1061 PetscCall(TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE));
1062 ts->dm = dmsave;
1063 PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
1064 PetscFunctionReturn(PETSC_SUCCESS);
1065 }
1066
SNESTSFormJacobian_GLLE(SNES snes,Vec x,Mat A,Mat B,TS ts)1067 static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts)
1068 {
1069 TS_GLLE *gl = (TS_GLLE *)ts->data;
1070 Vec Z, Ydot;
1071 DM dm, dmsave;
1072
1073 PetscFunctionBegin;
1074 PetscCall(SNESGetDM(snes, &dm));
1075 PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
1076 dmsave = ts->dm;
1077 ts->dm = dm;
1078 /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
1079 PetscCall(TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE));
1080 ts->dm = dmsave;
1081 PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
1082 PetscFunctionReturn(PETSC_SUCCESS);
1083 }
1084
TSSetUp_GLLE(TS ts)1085 static PetscErrorCode TSSetUp_GLLE(TS ts)
1086 {
1087 TS_GLLE *gl = (TS_GLLE *)ts->data;
1088 PetscInt max_r, max_s;
1089 DM dm;
1090
1091 PetscFunctionBegin;
1092 gl->setupcalled = PETSC_TRUE;
1093 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
1094 PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->X));
1095 PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold));
1096 PetscCall(VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot));
1097 PetscCall(VecDuplicateVecs(ts->vec_sol, 3, &gl->himom));
1098 PetscCall(VecDuplicate(ts->vec_sol, &gl->W));
1099 PetscCall(VecDuplicate(ts->vec_sol, &gl->Y));
1100 PetscCall(VecDuplicate(ts->vec_sol, &gl->Z));
1101
1102 /* Default acceptance tests and adaptivity */
1103 if (!gl->Accept) PetscCall(TSGLLESetAcceptType(ts, TSGLLEACCEPT_ALWAYS));
1104 if (!gl->adapt) PetscCall(TSGLLEGetAdapt(ts, &gl->adapt));
1105
1106 if (gl->current_scheme < 0) {
1107 PetscInt i;
1108 for (i = 0;; i++) {
1109 if (gl->schemes[i]->p == gl->start_order) break;
1110 PetscCheck(i + 1 != gl->nschemes, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No schemes available with requested start order %" PetscInt_FMT, i);
1111 }
1112 gl->current_scheme = i;
1113 }
1114 PetscCall(TSGetDM(ts, &dm));
1115 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
1116 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
1117 PetscFunctionReturn(PETSC_SUCCESS);
1118 }
1119
TSSetFromOptions_GLLE(TS ts,PetscOptionItems PetscOptionsObject)1120 static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems PetscOptionsObject)
1121 {
1122 TS_GLLE *gl = (TS_GLLE *)ts->data;
1123 char tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify";
1124
1125 PetscFunctionBegin;
1126 PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options");
1127 {
1128 PetscBool flg;
1129 PetscCall(PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg));
1130 if (flg || !gl->type_name[0]) PetscCall(TSGLLESetType(ts, tname));
1131 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));
1132 PetscCall(PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL));
1133 PetscCall(PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL));
1134 PetscCall(PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL));
1135 PetscCall(PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL));
1136 PetscCall(PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL));
1137 PetscCall(PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL));
1138 PetscCall(PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL));
1139 PetscCall(PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg));
1140 if (flg) {
1141 PetscBool match1, match2;
1142 PetscCall(PetscStrcmp(completef, "rescale", &match1));
1143 PetscCall(PetscStrcmp(completef, "rescale-and-modify", &match2));
1144 if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale;
1145 else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
1146 else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef);
1147 }
1148 {
1149 char type[256] = TSGLLEACCEPT_ALWAYS;
1150 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));
1151 if (flg || !gl->accept_name[0]) PetscCall(TSGLLESetAcceptType(ts, type));
1152 }
1153 {
1154 TSGLLEAdapt adapt;
1155 PetscCall(TSGLLEGetAdapt(ts, &adapt));
1156 PetscCall(TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject));
1157 }
1158 }
1159 PetscOptionsHeadEnd();
1160 PetscFunctionReturn(PETSC_SUCCESS);
1161 }
1162
TSView_GLLE(TS ts,PetscViewer viewer)1163 static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer)
1164 {
1165 TS_GLLE *gl = (TS_GLLE *)ts->data;
1166 PetscInt i;
1167 PetscBool isascii, details;
1168
1169 PetscFunctionBegin;
1170 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1171 if (isascii) {
1172 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));
1173 PetscCall(PetscViewerASCIIPrintf(viewer, " Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction]));
1174 PetscCall(PetscViewerASCIIPrintf(viewer, " Extrapolation: %s\n", gl->extrapolate ? "yes" : "no"));
1175 PetscCall(PetscViewerASCIIPrintf(viewer, " Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)"));
1176 PetscCall(PetscViewerASCIIPushTab(viewer));
1177 PetscCall(TSGLLEAdaptView(gl->adapt, viewer));
1178 PetscCall(PetscViewerASCIIPopTab(viewer));
1179 PetscCall(PetscViewerASCIIPrintf(viewer, " type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)"));
1180 PetscCall(PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes));
1181 details = PETSC_FALSE;
1182 PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL));
1183 PetscCall(PetscViewerASCIIPushTab(viewer));
1184 for (i = 0; i < gl->nschemes; i++) PetscCall(TSGLLESchemeView(gl->schemes[i], details, viewer));
1185 if (gl->View) PetscCall((*gl->View)(gl, viewer));
1186 PetscCall(PetscViewerASCIIPopTab(viewer));
1187 }
1188 PetscFunctionReturn(PETSC_SUCCESS);
1189 }
1190
1191 /*@C
1192 TSGLLERegister - adds a `TSGLLE` implementation
1193
1194 Not Collective, No Fortran Support
1195
1196 Input Parameters:
1197 + sname - name of user-defined general linear scheme
1198 - function - routine to create method context
1199
1200 Level: advanced
1201
1202 Note:
1203 `TSGLLERegister()` may be called multiple times to add several user-defined families.
1204
1205 Example Usage:
1206 .vb
1207 TSGLLERegister("my_scheme", MySchemeCreate);
1208 .ve
1209
1210 Then, your scheme can be chosen with the procedural interface via
1211 .vb
1212 TSGLLESetType(ts, "my_scheme")
1213 .ve
1214 or at runtime via the option
1215 .vb
1216 -ts_gl_type my_scheme
1217 .ve
1218
1219 .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`
1220 @*/
TSGLLERegister(const char sname[],PetscErrorCode (* function)(TS))1221 PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS))
1222 {
1223 PetscFunctionBegin;
1224 PetscCall(TSGLLEInitializePackage());
1225 PetscCall(PetscFunctionListAdd(&TSGLLEList, sname, function));
1226 PetscFunctionReturn(PETSC_SUCCESS);
1227 }
1228
1229 /*@C
1230 TSGLLEAcceptRegister - adds a `TSGLLE` acceptance scheme
1231
1232 Not Collective
1233
1234 Input Parameters:
1235 + sname - name of user-defined acceptance scheme
1236 - function - routine to create method context, see `TSGLLEAcceptFn` for the calling sequence
1237
1238 Level: advanced
1239
1240 Note:
1241 `TSGLLEAcceptRegister()` may be called multiple times to add several user-defined families.
1242
1243 Example Usage:
1244 .vb
1245 TSGLLEAcceptRegister("my_scheme", MySchemeCreate);
1246 .ve
1247
1248 Then, your scheme can be chosen with the procedural interface via
1249 .vb
1250 TSGLLESetAcceptType(ts, "my_scheme")
1251 .ve
1252 or at runtime via the option `-ts_gl_accept_type my_scheme`
1253
1254 .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLEAcceptFn`
1255 @*/
TSGLLEAcceptRegister(const char sname[],TSGLLEAcceptFn * function)1256 PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFn *function)
1257 {
1258 PetscFunctionBegin;
1259 PetscCall(PetscFunctionListAdd(&TSGLLEAcceptList, sname, function));
1260 PetscFunctionReturn(PETSC_SUCCESS);
1261 }
1262
1263 /*@C
1264 TSGLLERegisterAll - Registers all of the general linear methods in `TSGLLE`
1265
1266 Not Collective
1267
1268 Level: advanced
1269
1270 .seealso: [](ch_ts), `TSGLLE`, `TSGLLERegisterDestroy()`
1271 @*/
TSGLLERegisterAll(void)1272 PetscErrorCode TSGLLERegisterAll(void)
1273 {
1274 PetscFunctionBegin;
1275 if (TSGLLERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
1276 TSGLLERegisterAllCalled = PETSC_TRUE;
1277
1278 PetscCall(TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS));
1279 PetscCall(TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always));
1280 PetscFunctionReturn(PETSC_SUCCESS);
1281 }
1282
1283 /*@C
1284 TSGLLEInitializePackage - This function initializes everything in the `TSGLLE` package. It is called
1285 from `TSInitializePackage()`.
1286
1287 Level: developer
1288
1289 .seealso: [](ch_ts), `PetscInitialize()`, `TSInitializePackage()`, `TSGLLEFinalizePackage()`
1290 @*/
TSGLLEInitializePackage(void)1291 PetscErrorCode TSGLLEInitializePackage(void)
1292 {
1293 PetscFunctionBegin;
1294 if (TSGLLEPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
1295 TSGLLEPackageInitialized = PETSC_TRUE;
1296 PetscCall(TSGLLERegisterAll());
1297 PetscCall(PetscRegisterFinalize(TSGLLEFinalizePackage));
1298 PetscFunctionReturn(PETSC_SUCCESS);
1299 }
1300
1301 /*@C
1302 TSGLLEFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
1303 called from `PetscFinalize()`.
1304
1305 Level: developer
1306
1307 .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEInitializePackage()`, `TSInitializePackage()`
1308 @*/
TSGLLEFinalizePackage(void)1309 PetscErrorCode TSGLLEFinalizePackage(void)
1310 {
1311 PetscFunctionBegin;
1312 PetscCall(PetscFunctionListDestroy(&TSGLLEList));
1313 PetscCall(PetscFunctionListDestroy(&TSGLLEAcceptList));
1314 TSGLLEPackageInitialized = PETSC_FALSE;
1315 TSGLLERegisterAllCalled = PETSC_FALSE;
1316 PetscFunctionReturn(PETSC_SUCCESS);
1317 }
1318
1319 /*MC
1320 TSGLLE - DAE solver using implicit General Linear methods {cite}`butcher_2007` {cite}`butcher2016numerical`
1321
1322 Options Database Keys:
1323 + -ts_gl_type <type> - the class of general linear method (irks)
1324 . -ts_gl_rtol <tol> - relative error
1325 . -ts_gl_atol <tol> - absolute error
1326 . -ts_gl_min_order <p> - minimum order method to consider (default=1)
1327 . -ts_gl_max_order <p> - maximum order method to consider (default=3)
1328 . -ts_gl_start_order <p> - order of starting method (default=1)
1329 . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
1330 - -ts_adapt_type <method> - adaptive controller to use (none step both)
1331
1332 Level: beginner
1333
1334 Notes:
1335 These methods contain Runge-Kutta and multistep schemes as special cases. These special cases
1336 have some fundamental limitations. For example, diagonally implicit Runge-Kutta cannot have
1337 stage order greater than 1 which limits their applicability to very stiff systems.
1338 Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF are not
1339 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high
1340 stage order and reliable error estimates for both 1 and 2 orders higher to facilitate
1341 adaptive step sizes and adaptive order schemes. All this is possible while preserving a
1342 singly diagonally implicit structure.
1343
1344 This integrator can be applied to DAE.
1345
1346 Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit
1347 Runge-Kutta (DIRK). They are represented by the tableau
1348
1349 .vb
1350 A | U
1351 -------
1352 B | V
1353 .ve
1354
1355 combined with a vector c of abscissa. "Diagonally implicit" means that $A$ is lower
1356 triangular. A step of the general method reads
1357
1358 $$
1359 \begin{align*}
1360 [ Y ] = [A U] [ Y' ] \\
1361 [X^k] = [B V] [X^{k-1}]
1362 \end{align*}
1363 $$
1364
1365 where Y is the multivector of stage values, $Y'$ is the multivector of stage derivatives, $X^k$
1366 is the Nordsieck vector of the solution at step $k$. The Nordsieck vector consists of the first
1367 $r$ moments of the solution, given by
1368
1369 $$
1370 X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
1371 $$
1372
1373 If $A$ is lower triangular, we can solve the stages $(Y, Y')$ sequentially
1374
1375 $$
1376 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}
1377 $$
1378
1379 and then construct the pieces to carry to the next step
1380
1381 $$
1382 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}
1383 $$
1384
1385 Note that when the equations are cast in implicit form, we are using the stage equation to
1386 define $y'_i$ in terms of $y_i$ and known stuff ($y_j$ for $j<i$ and $x_j$ for all $j$).
1387
1388 Error estimation
1389
1390 At present, the most attractive GL methods for stiff problems are singly diagonally implicit
1391 schemes which posses Inherent Runge-Kutta Stability (`TSIRKS`). These methods have $r=s$, the
1392 number of items passed between steps is equal to the number of stages. The order and
1393 stage-order are one less than the number of stages. We use the error estimates in the 2007
1394 paper which provide the following estimates
1395
1396 $$
1397 \begin{align*}
1398 h^{p+1} X^{(p+1)} = \phi_0^T Y' + [0 \psi_0^T] Xold \\
1399 h^{p+2} X^{(p+2)} = \phi_1^T Y' + [0 \psi_1^T] Xold \\
1400 h^{p+2} (dx'/dx) X^{(p+1)} = \phi_2^T Y' + [0 \psi_2^T] Xold
1401 \end{align*}
1402 $$
1403
1404 These estimates are accurate to $ O(h^{p+3})$.
1405
1406 Changing the step size
1407
1408 Uses the generalized "rescale and modify" scheme, see equation (4.5) of {cite}`butcher_2007`.
1409
1410 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType`
1411 M*/
TSCreate_GLLE(TS ts)1412 PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
1413 {
1414 TS_GLLE *gl;
1415
1416 PetscFunctionBegin;
1417 PetscCall(TSGLLEInitializePackage());
1418
1419 PetscCall(PetscNew(&gl));
1420 ts->data = (void *)gl;
1421
1422 ts->ops->reset = TSReset_GLLE;
1423 ts->ops->destroy = TSDestroy_GLLE;
1424 ts->ops->view = TSView_GLLE;
1425 ts->ops->setup = TSSetUp_GLLE;
1426 ts->ops->solve = TSSolve_GLLE;
1427 ts->ops->setfromoptions = TSSetFromOptions_GLLE;
1428 ts->ops->snesfunction = SNESTSFormFunction_GLLE;
1429 ts->ops->snesjacobian = SNESTSFormJacobian_GLLE;
1430
1431 ts->usessnes = PETSC_TRUE;
1432
1433 gl->max_step_rejections = 1;
1434 gl->min_order = 1;
1435 gl->max_order = 3;
1436 gl->start_order = 1;
1437 gl->current_scheme = -1;
1438 gl->extrapolate = PETSC_FALSE;
1439
1440 gl->wrms_atol = 1e-8;
1441 gl->wrms_rtol = 1e-5;
1442
1443 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE));
1444 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE));
1445 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE));
1446 PetscFunctionReturn(PETSC_SUCCESS);
1447 }
1448