Actual source code: glle.c


  2: #include <../src/ts/impls/implicit/glle/glle.h>
  3: #include <petscdm.h>
  4: #include <petscblaslapack.h>

  6: static const char       *TSGLLEErrorDirections[] = {"FORWARD", "BACKWARD", "TSGLLEErrorDirection", "TSGLLEERROR_", NULL};
  7: static PetscFunctionList TSGLLEList;
  8: static PetscFunctionList TSGLLEAcceptList;
  9: static PetscBool         TSGLLEPackageInitialized;
 10: static PetscBool         TSGLLERegisterAllCalled;

 12: /* This function is pure */
 13: static PetscScalar Factorial(PetscInt n)
 14: {
 15:   PetscInt i;
 16:   if (n < 12) { /* Can compute with 32-bit integers */
 17:     PetscInt f = 1;
 18:     for (i = 2; i <= n; i++) f *= i;
 19:     return (PetscScalar)f;
 20:   } else {
 21:     PetscScalar f = 1.;
 22:     for (i = 2; i <= n; i++) f *= (PetscScalar)i;
 23:     return f;
 24:   }
 25: }

 27: /* This function is pure */
 28: static PetscScalar CPowF(PetscScalar c, PetscInt p)
 29: {
 30:   return PetscPowRealInt(PetscRealPart(c), p) / Factorial(p);
 31: }

 33: static PetscErrorCode TSGLLEGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
 34: {
 35:   TS_GLLE *gl = (TS_GLLE *)ts->data;

 37:   if (Z) {
 38:     if (dm && dm != ts->dm) {
 39:       DMGetNamedGlobalVector(dm, "TSGLLE_Z", Z);
 40:     } else *Z = gl->Z;
 41:   }
 42:   if (Ydotstage) {
 43:     if (dm && dm != ts->dm) {
 44:       DMGetNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage);
 45:     } else *Ydotstage = gl->Ydot[gl->stage];
 46:   }
 47:   return 0;
 48: }

 50: static PetscErrorCode TSGLLERestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
 51: {
 52:   if (Z) {
 53:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSGLLE_Z", Z);
 54:   }
 55:   if (Ydotstage) {
 56:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage);
 57:   }
 58:   return 0;
 59: }

 61: static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, void *ctx)
 62: {
 63:   return 0;
 64: }

 66: static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
 67: {
 68:   TS  ts = (TS)ctx;
 69:   Vec Ydot, Ydot_c;

 71:   TSGLLEGetVecs(ts, fine, NULL, &Ydot);
 72:   TSGLLEGetVecs(ts, coarse, NULL, &Ydot_c);
 73:   MatRestrict(restrct, Ydot, Ydot_c);
 74:   VecPointwiseMult(Ydot_c, rscale, Ydot_c);
 75:   TSGLLERestoreVecs(ts, fine, NULL, &Ydot);
 76:   TSGLLERestoreVecs(ts, coarse, NULL, &Ydot_c);
 77:   return 0;
 78: }

 80: static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, void *ctx)
 81: {
 82:   return 0;
 83: }

 85: static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
 86: {
 87:   TS  ts = (TS)ctx;
 88:   Vec Ydot, Ydot_s;

 90:   TSGLLEGetVecs(ts, dm, NULL, &Ydot);
 91:   TSGLLEGetVecs(ts, subdm, NULL, &Ydot_s);

 93:   VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD);
 94:   VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD);

 96:   TSGLLERestoreVecs(ts, dm, NULL, &Ydot);
 97:   TSGLLERestoreVecs(ts, subdm, NULL, &Ydot_s);
 98:   return 0;
 99: }

101: static PetscErrorCode TSGLLESchemeCreate(PetscInt p, PetscInt q, PetscInt r, PetscInt s, const PetscScalar *c, const PetscScalar *a, const PetscScalar *b, const PetscScalar *u, const PetscScalar *v, TSGLLEScheme *inscheme)
102: {
103:   TSGLLEScheme scheme;
104:   PetscInt     j;

110:   *inscheme = NULL;
111:   PetscNew(&scheme);
112:   scheme->p = p;
113:   scheme->q = q;
114:   scheme->r = r;
115:   scheme->s = s;

117:   PetscMalloc5(s, &scheme->c, s * s, &scheme->a, r * s, &scheme->b, r * s, &scheme->u, r * r, &scheme->v);
118:   PetscArraycpy(scheme->c, c, s);
119:   for (j = 0; j < s * s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j];
120:   for (j = 0; j < r * s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j];
121:   for (j = 0; j < s * r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j];
122:   for (j = 0; j < r * r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j];

124:   PetscMalloc6(r, &scheme->alpha, r, &scheme->beta, r, &scheme->gamma, 3 * s, &scheme->phi, 3 * r, &scheme->psi, r, &scheme->stage_error);
125:   {
126:     PetscInt     i, j, k, ss = s + 2;
127:     PetscBLASInt m, n, one = 1, *ipiv, lwork = 4 * ((s + 3) * 3 + 3), info, ldb;
128:     PetscReal    rcond, *sing, *workreal;
129:     PetscScalar *ImV, *H, *bmat, *workscalar, *c = scheme->c, *a = scheme->a, *b = scheme->b, *u = scheme->u, *v = scheme->v;
130:     PetscBLASInt rank;
131:     PetscMalloc7(PetscSqr(r), &ImV, 3 * s, &H, 3 * ss, &bmat, lwork, &workscalar, 5 * (3 + r), &workreal, r + s, &sing, r + s, &ipiv);

133:     /* column-major input */
134:     for (i = 0; i < r - 1; i++) {
135:       for (j = 0; j < r - 1; j++) ImV[i + j * r] = 1.0 * (i == j) - v[(i + 1) * r + j + 1];
136:     }
137:     /* Build right hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */
138:     for (i = 1; i < r; i++) {
139:       scheme->alpha[i] = 1. / Factorial(p + 1 - i);
140:       for (j = 0; j < s; j++) scheme->alpha[i] -= b[i * s + j] * CPowF(c[j], p);
141:     }
142:     PetscBLASIntCast(r - 1, &m);
143:     PetscBLASIntCast(r, &n);
144:     PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&m, &one, ImV, &n, ipiv, scheme->alpha + 1, &n, &info));

148:     /* Build right hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */
149:     for (i = 1; i < r; i++) {
150:       scheme->beta[i] = 1. / Factorial(p + 2 - i) - scheme->alpha[i];
151:       for (j = 0; j < s; j++) scheme->beta[i] -= b[i * s + j] * CPowF(c[j], p + 1);
152:     }
153:     PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->beta + 1, &n, &info));

157:     /* Build stage_error vector
158:            xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha;
159:     */
160:     for (i = 0; i < s; i++) {
161:       scheme->stage_error[i] = CPowF(c[i], p + 1);
162:       for (j = 0; j < s; j++) scheme->stage_error[i] -= a[i * s + j] * CPowF(c[j], p);
163:       for (j = 1; j < r; j++) scheme->stage_error[i] += u[i * r + j] * scheme->alpha[j];
164:     }

166:     /* alpha[0] (epsilon in B,J,W 2007)
167:            epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha;
168:     */
169:     scheme->alpha[0] = 1. / Factorial(p + 1);
170:     for (j = 0; j < s; j++) scheme->alpha[0] -= b[0 * s + j] * CPowF(c[j], p);
171:     for (j = 1; j < r; j++) scheme->alpha[0] += v[0 * r + j] * scheme->alpha[j];

173:     /* right hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */
174:     for (i = 1; i < r; i++) {
175:       scheme->gamma[i] = (i == 1 ? -1. : 0) * scheme->alpha[0];
176:       for (j = 0; j < s; j++) scheme->gamma[i] += b[i * s + j] * scheme->stage_error[j];
177:     }
178:     PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->gamma + 1, &n, &info));

182:     /* beta[0] (rho in B,J,W 2007)
183:         e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ...
184:             + glm.V(1,2:end)*e.beta;% - e.epsilon;
185:     % Note: The paper (B,J,W 2007) includes the last term in their definition
186:     * */
187:     scheme->beta[0] = 1. / Factorial(p + 2);
188:     for (j = 0; j < s; j++) scheme->beta[0] -= b[0 * s + j] * CPowF(c[j], p + 1);
189:     for (j = 1; j < r; j++) scheme->beta[0] += v[0 * r + j] * scheme->beta[j];

191:     /* gamma[0] (sigma in B,J,W 2007)
192:     *   e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma;
193:     * */
194:     scheme->gamma[0] = 0.0;
195:     for (j = 0; j < s; j++) scheme->gamma[0] += b[0 * s + j] * scheme->stage_error[j];
196:     for (j = 1; j < r; j++) scheme->gamma[0] += v[0 * s + j] * scheme->gamma[j];

198:     /* Assemble H
199:     *    % " PetscInt_FMT "etermine the error estimators phi
200:        H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ...
201:                [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]';
202:     % Paper has formula above without the 0, but that term must be left
203:     % out to satisfy the conditions they propose and to make the
204:     % example schemes work
205:     e.H = H;
206:     e.phi = (H \ [1 0 0;1 1 0;0 0 -1])';
207:     e.psi = -e.phi*C;
208:     * */
209:     for (j = 0; j < s; j++) {
210:       H[0 + j * 3] = CPowF(c[j], p);
211:       H[1 + j * 3] = CPowF(c[j], p + 1);
212:       H[2 + j * 3] = scheme->stage_error[j];
213:       for (k = 1; k < r; k++) {
214:         H[0 + j * 3] += CPowF(c[j], k - 1) * scheme->alpha[k];
215:         H[1 + j * 3] += CPowF(c[j], k - 1) * scheme->beta[k];
216:         H[2 + j * 3] -= CPowF(c[j], k - 1) * scheme->gamma[k];
217:       }
218:     }
219:     bmat[0 + 0 * ss] = 1.;
220:     bmat[0 + 1 * ss] = 0.;
221:     bmat[0 + 2 * ss] = 0.;
222:     bmat[1 + 0 * ss] = 1.;
223:     bmat[1 + 1 * ss] = 1.;
224:     bmat[1 + 2 * ss] = 0.;
225:     bmat[2 + 0 * ss] = 0.;
226:     bmat[2 + 1 * ss] = 0.;
227:     bmat[2 + 2 * ss] = -1.;
228:     m                = 3;
229:     PetscBLASIntCast(s, &n);
230:     PetscBLASIntCast(ss, &ldb);
231:     rcond = 1e-12;
232: #if defined(PETSC_USE_COMPLEX)
233:     /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
234:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, workreal, &info));
235: #else
236:     /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
237:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, &info));
238: #endif

242:     for (j = 0; j < 3; j++) {
243:       for (k = 0; k < s; k++) scheme->phi[k + j * s] = bmat[k + j * ss];
244:     }

246:     /* the other part of the error estimator, psi in B,J,W 2007 */
247:     scheme->psi[0 * r + 0] = 0.;
248:     scheme->psi[1 * r + 0] = 0.;
249:     scheme->psi[2 * r + 0] = 0.;
250:     for (j = 1; j < r; j++) {
251:       scheme->psi[0 * r + j] = 0.;
252:       scheme->psi[1 * r + j] = 0.;
253:       scheme->psi[2 * r + j] = 0.;
254:       for (k = 0; k < s; k++) {
255:         scheme->psi[0 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[0 * s + k];
256:         scheme->psi[1 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[1 * s + k];
257:         scheme->psi[2 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[2 * s + k];
258:       }
259:     }
260:     PetscFree7(ImV, H, bmat, workscalar, workreal, sing, ipiv);
261:   }
262:   /* Check which properties are satisfied */
263:   scheme->stiffly_accurate = PETSC_TRUE;
264:   if (scheme->c[s - 1] != 1.) scheme->stiffly_accurate = PETSC_FALSE;
265:   for (j = 0; j < s; j++)
266:     if (a[(s - 1) * s + j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE;
267:   for (j = 0; j < r; j++)
268:     if (u[(s - 1) * r + j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE;
269:   scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */
270:   for (j = 0; j < s - 1; j++)
271:     if (r > 1 && b[1 * s + j] != 0.) scheme->fsal = PETSC_FALSE;
272:   if (b[1 * s + r - 1] != 1.) scheme->fsal = PETSC_FALSE;
273:   for (j = 0; j < r; j++)
274:     if (r > 1 && v[1 * r + j] != 0.) scheme->fsal = PETSC_FALSE;

276:   *inscheme = scheme;
277:   return 0;
278: }

280: static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc)
281: {
282:   PetscFree5(sc->c, sc->a, sc->b, sc->u, sc->v);
283:   PetscFree6(sc->alpha, sc->beta, sc->gamma, sc->phi, sc->psi, sc->stage_error);
284:   PetscFree(sc);
285:   return 0;
286: }

288: static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
289: {
290:   PetscInt i;

292:   for (i = 0; i < gl->nschemes; i++) {
293:     if (gl->schemes[i]) TSGLLESchemeDestroy(gl->schemes[i]);
294:   }
295:   PetscFree(gl->schemes);
296:   gl->nschemes = 0;
297:   PetscMemzero(gl->type_name, sizeof(gl->type_name));
298:   return 0;
299: }

301: static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer, PetscInt m, PetscInt n, const PetscScalar a[], const char name[])
302: {
303:   PetscBool iascii;
304:   PetscInt  i, j;

306:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
307:   if (iascii) {
308:     PetscViewerASCIIPrintf(viewer, "%30s = [", name);
309:     for (i = 0; i < m; i++) {
310:       if (i) PetscViewerASCIIPrintf(viewer, "%30s   [", "");
311:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
312:       for (j = 0; j < n; j++) PetscViewerASCIIPrintf(viewer, " %12.8g", (double)PetscRealPart(a[i * n + j]));
313:       PetscViewerASCIIPrintf(viewer, "]\n");
314:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
315:     }
316:   }
317:   return 0;
318: }

320: static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer)
321: {
322:   PetscBool iascii;

324:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
325:   if (iascii) {
326:     PetscViewerASCIIPrintf(viewer, "GL scheme p,q,r,s = %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "\n", sc->p, sc->q, sc->r, sc->s);
327:     PetscViewerASCIIPushTab(viewer);
328:     PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s,  FSAL: %s\n", sc->stiffly_accurate ? "yes" : "no", sc->fsal ? "yes" : "no");
329:     PetscViewerASCIIPrintf(viewer, "Leading error constants: %10.3e  %10.3e  %10.3e\n", (double)PetscRealPart(sc->alpha[0]), (double)PetscRealPart(sc->beta[0]), (double)PetscRealPart(sc->gamma[0]));
330:     TSGLLEViewTable_Private(viewer, 1, sc->s, sc->c, "Abscissas c");
331:     if (view_details) {
332:       TSGLLEViewTable_Private(viewer, sc->s, sc->s, sc->a, "A");
333:       TSGLLEViewTable_Private(viewer, sc->r, sc->s, sc->b, "B");
334:       TSGLLEViewTable_Private(viewer, sc->s, sc->r, sc->u, "U");
335:       TSGLLEViewTable_Private(viewer, sc->r, sc->r, sc->v, "V");

337:       TSGLLEViewTable_Private(viewer, 3, sc->s, sc->phi, "Error estimate phi");
338:       TSGLLEViewTable_Private(viewer, 3, sc->r, sc->psi, "Error estimate psi");
339:       TSGLLEViewTable_Private(viewer, 1, sc->r, sc->alpha, "Modify alpha");
340:       TSGLLEViewTable_Private(viewer, 1, sc->r, sc->beta, "Modify beta");
341:       TSGLLEViewTable_Private(viewer, 1, sc->r, sc->gamma, "Modify gamma");
342:       TSGLLEViewTable_Private(viewer, 1, sc->s, sc->stage_error, "Stage error xi");
343:     }
344:     PetscViewerASCIIPopTab(viewer);
345:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported", ((PetscObject)viewer)->type_name);
346:   return 0;
347: }

349: static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[])
350: {
351:   PetscInt i;

354:   /* build error vectors*/
355:   for (i = 0; i < 3; i++) {
356:     PetscScalar phih[64];
357:     PetscInt    j;
358:     for (j = 0; j < sc->s; j++) phih[j] = sc->phi[i * sc->s + j] * h;
359:     VecZeroEntries(hm[i]);
360:     VecMAXPY(hm[i], sc->s, phih, Ydot);
361:     VecMAXPY(hm[i], sc->r, &sc->psi[i * sc->r], Xold);
362:   }
363:   return 0;
364: }

366: static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
367: {
368:   PetscScalar brow[32], vrow[32];
369:   PetscInt    i, j, r, s;

371:   /* Build the new solution from (X,Ydot) */
372:   r = sc->r;
373:   s = sc->s;
374:   for (i = 0; i < r; i++) {
375:     VecZeroEntries(X[i]);
376:     for (j = 0; j < s; j++) brow[j] = h * sc->b[i * s + j];
377:     VecMAXPY(X[i], s, brow, Ydot);
378:     for (j = 0; j < r; j++) vrow[j] = sc->v[i * r + j];
379:     VecMAXPY(X[i], r, vrow, Xold);
380:   }
381:   return 0;
382: }

384: static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
385: {
386:   PetscScalar brow[32], vrow[32];
387:   PetscReal   ratio;
388:   PetscInt    i, j, p, r, s;

390:   /* Build the new solution from (X,Ydot) */
391:   p     = sc->p;
392:   r     = sc->r;
393:   s     = sc->s;
394:   ratio = next_h / h;
395:   for (i = 0; i < r; i++) {
396:     VecZeroEntries(X[i]);
397:     for (j = 0; j < s; j++) {
398:       brow[j] = h * (PetscPowRealInt(ratio, i) * sc->b[i * s + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->phi[0 * s + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->phi[1 * s + j] + sc->gamma[i] * sc->phi[2 * s + j]));
399:     }
400:     VecMAXPY(X[i], s, brow, Ydot);
401:     for (j = 0; j < r; j++) {
402:       vrow[j] = (PetscPowRealInt(ratio, i) * sc->v[i * r + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->psi[0 * r + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->psi[1 * r + j] + sc->gamma[i] * sc->psi[2 * r + j]));
403:     }
404:     VecMAXPY(X[i], r, vrow, Xold);
405:   }
406:   if (r < next_sc->r) {
408:     VecZeroEntries(X[r]);
409:     for (j = 0; j < s; j++) brow[j] = h * PetscPowRealInt(ratio, p + 1) * sc->phi[0 * s + j];
410:     VecMAXPY(X[r], s, brow, Ydot);
411:     for (j = 0; j < r; j++) vrow[j] = PetscPowRealInt(ratio, p + 1) * sc->psi[0 * r + j];
412:     VecMAXPY(X[r], r, vrow, Xold);
413:   }
414:   return 0;
415: }

417: static PetscErrorCode TSGLLECreate_IRKS(TS ts)
418: {
419:   TS_GLLE *gl = (TS_GLLE *)ts->data;

421:   gl->Destroy               = TSGLLEDestroy_Default;
422:   gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default;
423:   gl->CompleteStep          = TSGLLECompleteStep_RescaleAndModify;
424:   PetscMalloc1(10, &gl->schemes);
425:   gl->nschemes = 0;

427:   {
428:     /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3
429:     * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE.
430:     * irks(0.3,0,[.3,1],[1],1)
431:     * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2)
432:     * but doing so would sacrifice the error estimator.
433:     */
434:     const PetscScalar c[2]    = {3. / 10., 1.};
435:     const PetscScalar a[2][2] = {
436:       {3. / 10., 0       },
437:       {7. / 10., 3. / 10.}
438:     };
439:     const PetscScalar b[2][2] = {
440:       {7. / 10., 3. / 10.},
441:       {0,        1       }
442:     };
443:     const PetscScalar u[2][2] = {
444:       {1, 0},
445:       {1, 0}
446:     };
447:     const PetscScalar v[2][2] = {
448:       {1, 0},
449:       {0, 0}
450:     };
451:     TSGLLESchemeCreate(1, 1, 2, 2, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
452:   }

454:   {
455:     /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */
456:     /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */
457:     const PetscScalar c[3]    = {1. / 3., 2. / 3., 1};
458:     const PetscScalar a[3][3] = {
459:       {4. / 9.,              0,                     0      },
460:       {1.03750643704090e+00, 4. / 9.,               0      },
461:       {7.67024779410304e-01, -3.81140216918943e-01, 4. / 9.}
462:     };
463:     const PetscScalar b[3][3] = {
464:       {0.767024779410304,  -0.381140216918943, 4. / 9.          },
465:       {0.000000000000000,  0.000000000000000,  1.000000000000000},
466:       {-2.075048385225385, 0.621728385225383,  1.277197204924873}
467:     };
468:     const PetscScalar u[3][3] = {
469:       {1.0000000000000000, -0.1111111111111109, -0.0925925925925922},
470:       {1.0000000000000000, -0.8152842148186744, -0.4199095530877056},
471:       {1.0000000000000000, 0.1696709930641948,  0.0539741070314165 }
472:     };
473:     const PetscScalar v[3][3] = {
474:       {1.0000000000000000, 0.1696709930641948, 0.0539741070314165},
475:       {0.000000000000000,  0.000000000000000,  0.000000000000000 },
476:       {0.000000000000000,  0.176122795075129,  0.000000000000000 }
477:     };
478:     TSGLLESchemeCreate(2, 2, 3, 3, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
479:   }
480:   {
481:     /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */
482:     const PetscScalar c[4]    = {0.25, 0.5, 0.75, 1.0};
483:     const PetscScalar a[4][4] = {
484:       {9. / 40.,             0,                     0,                 0       },
485:       {2.11286958887701e-01, 9. / 40.,              0,                 0       },
486:       {9.46338294287584e-01, -3.42942861246094e-01, 9. / 40.,          0       },
487:       {0.521490453970721,    -0.662474225622980,    0.490476425459734, 9. / 40.}
488:     };
489:     const PetscScalar b[4][4] = {
490:       {0.521490453970721,  -0.662474225622980, 0.490476425459734,  9. / 40.         },
491:       {0.000000000000000,  0.000000000000000,  0.000000000000000,  1.000000000000000},
492:       {-0.084677029310348, 1.390757514776085,  -1.568157386206001, 2.023192696767826},
493:       {0.465383797936408,  1.478273530625148,  -1.930836081010182, 1.644872111193354}
494:     };
495:     const PetscScalar u[4][4] = {
496:       {1.00000000000000000, 0.02500000000001035,  -0.02499999999999053, -0.00442708333332865},
497:       {1.00000000000000000, 0.06371304111232945,  -0.04032173972189845, -0.01389438413189452},
498:       {1.00000000000000000, -0.07839543304147778, 0.04738685705116663,  0.02032603595928376 },
499:       {1.00000000000000000, 0.42550734619251651,  0.10800718022400080,  -0.01726712647760034}
500:     };
501:     const PetscScalar v[4][4] = {
502:       {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034},
503:       {0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000   },
504:       {0.000000000000000,   -1.761115796027561,  -0.521284157173780,  0.258249384305463   },
505:       {0.000000000000000,   -1.657693358744728,  -1.052227765232394,  0.521284157173780   }
506:     };
507:     TSGLLESchemeCreate(3, 3, 4, 4, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
508:   }
509:   {
510:     /* p=q=4, r=s=5:
511:           irks(3/11,0,[1:5]/5, [0.1715   -0.1238    0.6617],...
512:           [ -0.0812    0.4079    1.0000
513:              1.0000         0         0
514:              0.8270    1.0000         0])
515:     */
516:     const PetscScalar c[5]    = {0.2, 0.4, 0.6, 0.8, 1.0};
517:     const PetscScalar a[5][5] = {
518:       {2.72727272727352e-01,  0.00000000000000e+00,  0.00000000000000e+00,  0.00000000000000e+00, 0.00000000000000e+00},
519:       {-1.03980153733431e-01, 2.72727272727405e-01,  0.00000000000000e+00,  0.00000000000000e+00, 0.00000000000000e+00},
520:       {-1.58615400341492e+00, 7.44168951881122e-01,  2.72727272727309e-01,  0.00000000000000e+00, 0.00000000000000e+00},
521:       {-8.73658042865628e-01, 5.37884671894595e-01,  -1.63298538799523e-01, 2.72727272726996e-01, 0.00000000000000e+00},
522:       {2.95489397443992e-01,  -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01}
523:     };
524:     const PetscScalar b[5][5] = {
525:       {2.95489397443992e-01,  -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00,  2.72727272727288e-01},
526:       {0.00000000000000e+00,  1.11022302462516e-16,  -2.22044604925031e-16, 0.00000000000000e+00,  1.00000000000000e+00},
527:       {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00,  6.32331093108427e-01},
528:       {8.35690179937017e+00,  -2.26640927349732e+00, 6.86647884973826e+00,  -5.22595158025740e+00, 4.50893068837431e+00},
529:       {1.27656267027479e+01,  2.80882153840821e+00,  8.91173096522890e+00,  -1.07936444078906e+01, 4.82534148988854e+00}
530:     };
531:     const PetscScalar u[5][5] = {
532:       {1.00000000000000e+00, -7.27272727273551e-02, -3.45454545454419e-02, -4.12121212119565e-03, -2.96969696964014e-04},
533:       {1.00000000000000e+00, 2.31252881006154e-01,  -8.29487834416481e-03, -9.07191207681020e-03, -1.70378403743473e-03},
534:       {1.00000000000000e+00, 1.16925777880663e+00,  3.59268562942635e-02,  -4.09013451730615e-02, -1.02411119670164e-02},
535:       {1.00000000000000e+00, 1.02634463704356e+00,  1.59375044913405e-01,  1.89673015035370e-03,  -4.89987231897569e-03},
536:       {1.00000000000000e+00, 1.27746320298021e+00,  2.37186008132728e-01,  -8.28694373940065e-02, -5.34396510196430e-02}
537:     };
538:     const PetscScalar v[5][5] = {
539:       {1.00000000000000e+00, 1.27746320298021e+00,  2.37186008132728e-01,  -8.28694373940065e-02, -5.34396510196430e-02},
540:       {0.00000000000000e+00, -1.77635683940025e-15, -1.99840144432528e-15, -9.99200722162641e-16, -3.33066907387547e-16},
541:       {0.00000000000000e+00, 4.37280081906924e+00,  5.49221645016377e-02,  -8.88913177394943e-02, 1.12879077989154e-01 },
542:       {0.00000000000000e+00, -1.22399504837280e+01, -5.21287338448645e+00, -8.03952325565291e-01, 4.60298678047147e-01 },
543:       {0.00000000000000e+00, -1.85178762883829e+01, -5.21411849862624e+00, -1.04283436528809e+00, 7.49030161063651e-01 }
544:     };
545:     TSGLLESchemeCreate(4, 4, 5, 5, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
546:   }
547:   {
548:     /* p=q=5, r=s=6;
549:        irks(1/3,0,[1:6]/6,...
550:           [-0.0489    0.4228   -0.8814    0.9021],...
551:           [-0.3474   -0.6617    0.6294    0.2129
552:             0.0044   -0.4256   -0.1427   -0.8936
553:            -0.8267    0.4821    0.1371   -0.2557
554:            -0.4426   -0.3855   -0.7514    0.3014])
555:     */
556:     const PetscScalar c[6]    = {1. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1.};
557:     const PetscScalar a[6][6] = {
558:       {3.33333333333940e-01,  0,                     0,                     0,                     0,                    0                   },
559:       {-8.64423857333350e-02, 3.33333333332888e-01,  0,                     0,                     0,                    0                   },
560:       {-2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01,  0,                     0,                    0                   },
561:       {-4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01,  0,                    0                   },
562:       {-6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01,  -4.48352364517632e-01, 3.33333333328483e-01, 0                   },
563:       {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00,  -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}
564:     };
565:     const PetscScalar b[6][6] = {
566:       {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00,  -2.23969848002481e+00, 6.62807710124007e-01,  3.33333333335440e-01 },
567:       {-8.88178419700125e-16, 4.44089209850063e-16,  -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00,  1.00000000000001e+00 },
568:       {-2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01,  2.56943874812797e+01,  -3.06702268304488e+01, 6.68067773510103e+00 },
569:       {5.47971245256474e+01,  6.80366875868284e+01,  -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01,  -1.17819043489036e+01},
570:       {-2.33332114788869e+02, 6.12942539462634e+01,  -4.91850135865944e+01, 1.82716844135480e+02,  -1.29788173979395e+02, 3.09968095651099e+01 },
571:       {-1.72049132343751e+02, 8.60194713593999e+00,  7.98154219170200e-01,  1.50371386053218e+02,  -1.18515423962066e+02, 2.50898277784663e+01 }
572:     };
573:     const PetscScalar u[6][6] = {
574:       {1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06},
575:       {1.00000000000000e+00, 8.64423857327162e-02,  -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04},
576:       {1.00000000000000e+00, 4.57135912953434e+00,  1.06514719719137e+00,  1.33517564218007e-01,  1.11365952968659e-02,  6.12382756769504e-04 },
577:       {1.00000000000000e+00, 9.23391519753404e+00,  2.22431212392095e+00,  2.91823807741891e-01,  2.52058456411084e-02,  1.22800542949647e-03 },
578:       {1.00000000000000e+00, 1.48175480533865e+01,  3.73439117461835e+00,  5.14648336541804e-01,  4.76430038853402e-02,  2.56798515502156e-03 },
579:       {1.00000000000000e+00, 1.50512347758335e+01,  4.10099701165164e+00,  5.66039141003603e-01,  3.91213893800891e-02,  -2.99136269067853e-03}
580:     };
581:     const PetscScalar v[6][6] = {
582:       {1.00000000000000e+00, 1.50512347758335e+01,  4.10099701165164e+00,  5.66039141003603e-01,  3.91213893800891e-02,  -2.99136269067853e-03},
583:       {0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16},
584:       {0.00000000000000e+00, 1.22250171233141e+01,  -1.77150760606169e+00, 3.54516769879390e-01,  6.22298845883398e-01,  2.31647447450276e-01 },
585:       {0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01,  6.55727990241799e-02,  1.63175368287079e-01 },
586:       {0.00000000000000e+00, 1.37297394708005e+02,  -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01,  9.16629423682464e-01 },
587:       {0.00000000000000e+00, 1.05703241119022e+02,  -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00 }
588:     };
589:     TSGLLESchemeCreate(5, 5, 6, 6, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]);
590:   }
591:   return 0;
592: }

594: /*@C
595:    TSGLLESetType - sets the class of general linear method to use for time-stepping

597:    Collective on TS

599:    Input Parameters:
600: +  ts - the TS context
601: -  type - a method

603:    Options Database Key:
604: .  -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks)

606:    Notes:
607:    See "petsc/include/petscts.h" for available methods (for instance)
608: .    TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems)

610:    Normally, it is best to use the TSSetFromOptions() command and
611:    then set the TSGLLE type from the options database rather than by using
612:    this routine.  Using the options database provides the user with
613:    maximum flexibility in evaluating the many different solvers.
614:    The TSGLLESetType() routine is provided for those situations where it
615:    is necessary to set the timestepping solver independently of the
616:    command line or options database.  This might be the case, for example,
617:    when the choice of solver changes during the execution of the
618:    program, and the user's application is taking responsibility for
619:    choosing the appropriate method.

621:    Level: intermediate

623: @*/
624: PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type)
625: {
628:   PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type));
629:   return 0;
630: }

632: /*@C
633:    TSGLLESetAcceptType - sets the acceptance test

635:    Time integrators that need to control error must have the option to reject a time step based on local error
636:    estimates.  This function allows different schemes to be set.

638:    Logically Collective on TS

640:    Input Parameters:
641: +  ts - the TS context
642: -  type - the type

644:    Options Database Key:
645: .  -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step

647:    Level: intermediate

649: .seealso: `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt`, `set` `type`
650: @*/
651: PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type)
652: {
655:   PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type));
656:   return 0;
657: }

659: /*@C
660:    TSGLLEGetAdapt - gets the TSGLLEAdapt object from the TS

662:    Not Collective

664:    Input Parameter:
665: .  ts - the TS context

667:    Output Parameter:
668: .  adapt - the TSGLLEAdapt context

670:    Notes:
671:    This allows the user set options on the TSGLLEAdapt object.  Usually it is better to do this using the options
672:    database, so this function is rarely needed.

674:    Level: advanced

676: .seealso: `TSGLLEAdapt`, `TSGLLEAdaptRegister()`
677: @*/
678: PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt)
679: {
682:   PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt));
683:   return 0;
684: }

686: static PetscErrorCode TSGLLEAccept_Always(TS ts, PetscReal tleft, PetscReal h, const PetscReal enorms[], PetscBool *accept)
687: {
688:   *accept = PETSC_TRUE;
689:   return 0;
690: }

692: static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
693: {
694:   TS_GLLE     *gl = (TS_GLLE *)ts->data;
695:   PetscScalar *x, *w;
696:   PetscInt     n, i;

698:   VecGetArray(gl->X[0], &x);
699:   VecGetArray(gl->W, &w);
700:   VecGetLocalSize(gl->W, &n);
701:   for (i = 0; i < n; i++) w[i] = 1. / (gl->wrms_atol + gl->wrms_rtol * PetscAbsScalar(x[i]));
702:   VecRestoreArray(gl->X[0], &x);
703:   VecRestoreArray(gl->W, &w);
704:   return 0;
705: }

707: static PetscErrorCode TSGLLEVecNormWRMS(TS ts, Vec X, PetscReal *nrm)
708: {
709:   TS_GLLE     *gl = (TS_GLLE *)ts->data;
710:   PetscScalar *x, *w;
711:   PetscReal    sum = 0.0, gsum;
712:   PetscInt     n, N, i;

714:   VecGetArray(X, &x);
715:   VecGetArray(gl->W, &w);
716:   VecGetLocalSize(gl->W, &n);
717:   for (i = 0; i < n; i++) sum += PetscAbsScalar(PetscSqr(x[i] * w[i]));
718:   VecRestoreArray(X, &x);
719:   VecRestoreArray(gl->W, &w);
720:   MPIU_Allreduce(&sum, &gsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));
721:   VecGetSize(gl->W, &N);
722:   *nrm = PetscSqrtReal(gsum / (1. * N));
723:   return 0;
724: }

726: static PetscErrorCode TSGLLESetType_GLLE(TS ts, TSGLLEType type)
727: {
728:   PetscBool same;
729:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
730:   PetscErrorCode (*r)(TS);

732:   if (gl->type_name[0]) {
733:     PetscStrcmp(gl->type_name, type, &same);
734:     if (same) return 0;
735:     (*gl->Destroy)(gl);
736:   }

738:   PetscFunctionListFind(TSGLLEList, type, &r);
740:   (*r)(ts);
741:   PetscStrcpy(gl->type_name, type);
742:   return 0;
743: }

745: static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type)
746: {
747:   TSGLLEAcceptFunction r;
748:   TS_GLLE             *gl = (TS_GLLE *)ts->data;

750:   PetscFunctionListFind(TSGLLEAcceptList, type, &r);
752:   gl->Accept = r;
753:   PetscStrncpy(gl->accept_name, type, sizeof(gl->accept_name));
754:   return 0;
755: }

757: static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt)
758: {
759:   TS_GLLE *gl = (TS_GLLE *)ts->data;

761:   if (!gl->adapt) {
762:     TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts), &gl->adapt);
763:     PetscObjectIncrementTabLevel((PetscObject)gl->adapt, (PetscObject)ts, 1);
764:   }
765:   *adapt = gl->adapt;
766:   return 0;
767: }

769: static PetscErrorCode TSGLLEChooseNextScheme(TS ts, PetscReal h, const PetscReal hmnorm[], PetscInt *next_scheme, PetscReal *next_h, PetscBool *finish)
770: {
771:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
772:   PetscInt  i, n, cur_p, cur, next_sc, candidates[64], orders[64];
773:   PetscReal errors[64], costs[64], tleft;

775:   cur   = -1;
776:   cur_p = gl->schemes[gl->current_scheme]->p;
777:   tleft = ts->max_time - (ts->ptime + ts->time_step);
778:   for (i = 0, n = 0; i < gl->nschemes; i++) {
779:     TSGLLEScheme sc = gl->schemes[i];
780:     if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
781:     if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[0];
782:     else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[1];
783:     else if (sc->p == cur_p + 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * (hmnorm[2] + hmnorm[3]);
784:     else continue;
785:     candidates[n] = i;
786:     orders[n]     = PetscMin(sc->p, sc->q); /* order of global truncation error */
787:     costs[n]      = sc->s;                  /* estimate the cost as the number of stages */
788:     if (i == gl->current_scheme) cur = n;
789:     n++;
790:   }
792:   TSGLLEAdaptChoose(gl->adapt, n, orders, errors, costs, cur, h, tleft, &next_sc, next_h, finish);
793:   *next_scheme = candidates[next_sc];
794:   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,
795:                       gl->schemes[*next_scheme]->r, gl->schemes[*next_scheme]->s, (double)*next_h, PetscBools[*finish]));
796:   return 0;
797: }

799: static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s)
800: {
801:   TS_GLLE *gl = (TS_GLLE *)ts->data;

803:   *max_r = gl->schemes[gl->nschemes - 1]->r;
804:   *max_s = gl->schemes[gl->nschemes - 1]->s;
805:   return 0;
806: }

808: static PetscErrorCode TSSolve_GLLE(TS ts)
809: {
810:   TS_GLLE            *gl = (TS_GLLE *)ts->data;
811:   PetscInt            i, k, its, lits, max_r, max_s;
812:   PetscBool           final_step, finish;
813:   SNESConvergedReason snesreason;

815:   TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);

817:   TSGLLEGetMaxSizes(ts, &max_r, &max_s);
818:   VecCopy(ts->vec_sol, gl->X[0]);
819:   for (i = 1; i < max_r; i++) VecZeroEntries(gl->X[i]);
820:   TSGLLEUpdateWRMS(ts);

822:   if (0) {
823:     /* Find consistent initial data for DAE */
824:     gl->stage_time = ts->ptime + ts->time_step;
825:     gl->scoeff     = 1.;
826:     gl->stage      = 0;

828:     VecCopy(ts->vec_sol, gl->Z);
829:     VecCopy(ts->vec_sol, gl->Y);
830:     SNESSolve(ts->snes, NULL, gl->Y);
831:     SNESGetIterationNumber(ts->snes, &its);
832:     SNESGetLinearSolveIterations(ts->snes, &lits);
833:     SNESGetConvergedReason(ts->snes, &snesreason);

835:     ts->snes_its += its;
836:     ts->ksp_its += lits;
837:     if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
838:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
839:       PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures);
840:       return 0;
841:     }
842:   }


846:   for (k = 0, final_step = PETSC_FALSE, finish = PETSC_FALSE; k < ts->max_steps && !finish; k++) {
847:     PetscInt           j, r, s, next_scheme = 0, rejections;
848:     PetscReal          h, hmnorm[4], enorm[3], next_h;
849:     PetscBool          accept;
850:     const PetscScalar *c, *a, *u;
851:     Vec               *X, *Ydot, Y;
852:     TSGLLEScheme       scheme = gl->schemes[gl->current_scheme];

854:     r    = scheme->r;
855:     s    = scheme->s;
856:     c    = scheme->c;
857:     a    = scheme->a;
858:     u    = scheme->u;
859:     h    = ts->time_step;
860:     X    = gl->X;
861:     Ydot = gl->Ydot;
862:     Y    = gl->Y;

864:     if (ts->ptime > ts->max_time) break;

866:     /*
867:       We only call PreStep at the start of each STEP, not each STAGE.  This is because it is
868:       possible to fail (have to restart a step) after multiple stages.
869:     */
870:     TSPreStep(ts);

872:     rejections = 0;
873:     while (1) {
874:       for (i = 0; i < s; i++) {
875:         PetscScalar shift;
876:         gl->scoeff     = 1. / PetscRealPart(a[i * s + i]);
877:         shift          = gl->scoeff / ts->time_step;
878:         gl->stage      = i;
879:         gl->stage_time = ts->ptime + PetscRealPart(c[i]) * h;

881:         /*
882:         * Stage equation: Y = h A Y' + U X
883:         * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
884:         * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
885:         * Then y'_i = z + 1/(h a_ii) y_i
886:         */
887:         VecZeroEntries(gl->Z);
888:         for (j = 0; j < r; j++) VecAXPY(gl->Z, -shift * u[i * r + j], X[j]);
889:         for (j = 0; j < i; j++) VecAXPY(gl->Z, -shift * h * a[i * s + j], Ydot[j]);
890:         /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */

892:         /* Compute an estimate of Y to start Newton iteration */
893:         if (gl->extrapolate) {
894:           if (i == 0) {
895:             /* Linear extrapolation on the first stage */
896:             VecWAXPY(Y, c[i] * h, X[1], X[0]);
897:           } else {
898:             /* Linear extrapolation from the last stage */
899:             VecAXPY(Y, (c[i] - c[i - 1]) * h, Ydot[i - 1]);
900:           }
901:         } else if (i == 0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
902:           VecCopy(X[0], Y);
903:         }

905:         /* Solve this stage (Ydot[i] is computed during function evaluation) */
906:         SNESSolve(ts->snes, NULL, Y);
907:         SNESGetIterationNumber(ts->snes, &its);
908:         SNESGetLinearSolveIterations(ts->snes, &lits);
909:         SNESGetConvergedReason(ts->snes, &snesreason);
910:         ts->snes_its += its;
911:         ts->ksp_its += lits;
912:         if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
913:           ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
914:           PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures);
915:           return 0;
916:         }
917:       }

919:       gl->stage_time = ts->ptime + ts->time_step;

921:       (*gl->EstimateHigherMoments)(scheme, h, Ydot, gl->X, gl->himom);
922:       /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */
923:       for (i = 0; i < 3; i++) TSGLLEVecNormWRMS(ts, gl->himom[i], &hmnorm[i + 1]);
924:       enorm[0] = PetscRealPart(scheme->alpha[0]) * hmnorm[1];
925:       enorm[1] = PetscRealPart(scheme->beta[0]) * hmnorm[2];
926:       enorm[2] = PetscRealPart(scheme->gamma[0]) * hmnorm[3];
927:       (*gl->Accept)(ts, ts->max_time - gl->stage_time, h, enorm, &accept);
928:       if (accept) goto accepted;
929:       rejections++;
930:       PetscInfo(ts, "Step %" PetscInt_FMT " (t=%g) not accepted, rejections=%" PetscInt_FMT "\n", k, (double)gl->stage_time, rejections);
931:       if (rejections > gl->max_step_rejections) break;
932:       /*
933:         There are lots of reasons why a step might be rejected, including solvers not converging and other factors that
934:         TSGLLEChooseNextScheme does not support.  Additionally, the error estimates may be very screwed up, so I'm not
935:         convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor
936:         (the adaptor interface probably has to change).  Here we make an arbitrary and naive choice.  This assumes that
937:         steps were written in Nordsieck form.  The "correct" method would be to re-complete the previous time step with
938:         the correct "next" step size.  It is unclear to me whether the present ad-hoc method of rescaling X is stable.
939:       */
940:       h *= 0.5;
941:       for (i = 1; i < scheme->r; i++) VecScale(X[i], PetscPowRealInt(0.5, i));
942:     }
943:     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);

945:   accepted:
946:     /* This term is not error, but it *would* be the leading term for a lower order method */
947:     TSGLLEVecNormWRMS(ts, gl->X[scheme->r - 1], &hmnorm[0]);
948:     /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */

950:     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]);
951:     if (!final_step) {
952:       TSGLLEChooseNextScheme(ts, h, hmnorm, &next_scheme, &next_h, &final_step);
953:     } else {
954:       /* Dummy values to complete the current step in a consistent manner */
955:       next_scheme = gl->current_scheme;
956:       next_h      = h;
957:       finish      = PETSC_TRUE;
958:     }

960:     X        = gl->Xold;
961:     gl->Xold = gl->X;
962:     gl->X    = X;
963:     (*gl->CompleteStep)(scheme, h, gl->schemes[next_scheme], next_h, Ydot, gl->Xold, gl->X);

965:     TSGLLEUpdateWRMS(ts);

967:     /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */
968:     VecCopy(gl->X[0], ts->vec_sol);
969:     ts->ptime += h;
970:     ts->steps++;

972:     TSPostEvaluate(ts);
973:     TSPostStep(ts);
974:     TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);

976:     gl->current_scheme = next_scheme;
977:     ts->time_step      = next_h;
978:   }
979:   return 0;
980: }

982: /*------------------------------------------------------------*/

984: static PetscErrorCode TSReset_GLLE(TS ts)
985: {
986:   TS_GLLE *gl = (TS_GLLE *)ts->data;
987:   PetscInt max_r, max_s;

989:   if (gl->setupcalled) {
990:     TSGLLEGetMaxSizes(ts, &max_r, &max_s);
991:     VecDestroyVecs(max_r, &gl->Xold);
992:     VecDestroyVecs(max_r, &gl->X);
993:     VecDestroyVecs(max_s, &gl->Ydot);
994:     VecDestroyVecs(3, &gl->himom);
995:     VecDestroy(&gl->W);
996:     VecDestroy(&gl->Y);
997:     VecDestroy(&gl->Z);
998:   }
999:   gl->setupcalled = PETSC_FALSE;
1000:   return 0;
1001: }

1003: static PetscErrorCode TSDestroy_GLLE(TS ts)
1004: {
1005:   TS_GLLE *gl = (TS_GLLE *)ts->data;

1007:   TSReset_GLLE(ts);
1008:   if (ts->dm) {
1009:     DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts);
1010:     DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts);
1011:   }
1012:   if (gl->adapt) TSGLLEAdaptDestroy(&gl->adapt);
1013:   if (gl->Destroy) (*gl->Destroy)(gl);
1014:   PetscFree(ts->data);
1015:   PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL);
1016:   PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL);
1017:   PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL);
1018:   return 0;
1019: }

1021: /*
1022:     This defines the nonlinear equation that is to be solved with SNES
1023:     g(x) = f(t,x,z+shift*x) = 0
1024: */
1025: static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts)
1026: {
1027:   TS_GLLE *gl = (TS_GLLE *)ts->data;
1028:   Vec      Z, Ydot;
1029:   DM       dm, dmsave;

1031:   SNESGetDM(snes, &dm);
1032:   TSGLLEGetVecs(ts, dm, &Z, &Ydot);
1033:   VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z);
1034:   dmsave = ts->dm;
1035:   ts->dm = dm;
1036:   TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE);
1037:   ts->dm = dmsave;
1038:   TSGLLERestoreVecs(ts, dm, &Z, &Ydot);
1039:   return 0;
1040: }

1042: static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts)
1043: {
1044:   TS_GLLE *gl = (TS_GLLE *)ts->data;
1045:   Vec      Z, Ydot;
1046:   DM       dm, dmsave;

1048:   SNESGetDM(snes, &dm);
1049:   TSGLLEGetVecs(ts, dm, &Z, &Ydot);
1050:   dmsave = ts->dm;
1051:   ts->dm = dm;
1052:   /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
1053:   TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE);
1054:   ts->dm = dmsave;
1055:   TSGLLERestoreVecs(ts, dm, &Z, &Ydot);
1056:   return 0;
1057: }

1059: static PetscErrorCode TSSetUp_GLLE(TS ts)
1060: {
1061:   TS_GLLE *gl = (TS_GLLE *)ts->data;
1062:   PetscInt max_r, max_s;
1063:   DM       dm;

1065:   gl->setupcalled = PETSC_TRUE;
1066:   TSGLLEGetMaxSizes(ts, &max_r, &max_s);
1067:   VecDuplicateVecs(ts->vec_sol, max_r, &gl->X);
1068:   VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold);
1069:   VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot);
1070:   VecDuplicateVecs(ts->vec_sol, 3, &gl->himom);
1071:   VecDuplicate(ts->vec_sol, &gl->W);
1072:   VecDuplicate(ts->vec_sol, &gl->Y);
1073:   VecDuplicate(ts->vec_sol, &gl->Z);

1075:   /* Default acceptance tests and adaptivity */
1076:   if (!gl->Accept) TSGLLESetAcceptType(ts, TSGLLEACCEPT_ALWAYS);
1077:   if (!gl->adapt) TSGLLEGetAdapt(ts, &gl->adapt);

1079:   if (gl->current_scheme < 0) {
1080:     PetscInt i;
1081:     for (i = 0;; i++) {
1082:       if (gl->schemes[i]->p == gl->start_order) break;
1084:     }
1085:     gl->current_scheme = i;
1086:   }
1087:   TSGetDM(ts, &dm);
1088:   DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts);
1089:   DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts);
1090:   return 0;
1091: }
1092: /*------------------------------------------------------------*/

1094: static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems *PetscOptionsObject)
1095: {
1096:   TS_GLLE *gl         = (TS_GLLE *)ts->data;
1097:   char     tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify";

1099:   PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options");
1100:   {
1101:     PetscBool flg;
1102:     PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg);
1103:     if (flg || !gl->type_name[0]) TSGLLESetType(ts, tname);
1104:     PetscOptionsInt("-ts_gl_max_step_rejections", "Maximum number of times to attempt a step", "None", gl->max_step_rejections, &gl->max_step_rejections, NULL);
1105:     PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL);
1106:     PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL);
1107:     PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL);
1108:     PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL);
1109:     PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL);
1110:     PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL);
1111:     PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL);
1112:     PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg);
1113:     if (flg) {
1114:       PetscBool match1, match2;
1115:       PetscStrcmp(completef, "rescale", &match1);
1116:       PetscStrcmp(completef, "rescale-and-modify", &match2);
1117:       if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale;
1118:       else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
1119:       else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef);
1120:     }
1121:     {
1122:       char type[256] = TSGLLEACCEPT_ALWAYS;
1123:       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);
1124:       if (flg || !gl->accept_name[0]) TSGLLESetAcceptType(ts, type);
1125:     }
1126:     {
1127:       TSGLLEAdapt adapt;
1128:       TSGLLEGetAdapt(ts, &adapt);
1129:       TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject);
1130:     }
1131:   }
1132:   PetscOptionsHeadEnd();
1133:   return 0;
1134: }

1136: static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer)
1137: {
1138:   TS_GLLE  *gl = (TS_GLLE *)ts->data;
1139:   PetscInt  i;
1140:   PetscBool iascii, details;

1142:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1143:   if (iascii) {
1144:     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);
1145:     PetscViewerASCIIPrintf(viewer, "  Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction]);
1146:     PetscViewerASCIIPrintf(viewer, "  Extrapolation: %s\n", gl->extrapolate ? "yes" : "no");
1147:     PetscViewerASCIIPrintf(viewer, "  Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)");
1148:     PetscViewerASCIIPushTab(viewer);
1149:     TSGLLEAdaptView(gl->adapt, viewer);
1150:     PetscViewerASCIIPopTab(viewer);
1151:     PetscViewerASCIIPrintf(viewer, "  type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)");
1152:     PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes);
1153:     details = PETSC_FALSE;
1154:     PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL);
1155:     PetscViewerASCIIPushTab(viewer);
1156:     for (i = 0; i < gl->nschemes; i++) TSGLLESchemeView(gl->schemes[i], details, viewer);
1157:     if (gl->View) (*gl->View)(gl, viewer);
1158:     PetscViewerASCIIPopTab(viewer);
1159:   }
1160:   return 0;
1161: }

1163: /*@C
1164:    TSGLLERegister -  adds a TSGLLE implementation

1166:    Not Collective

1168:    Input Parameters:
1169: +  name_scheme - name of user-defined general linear scheme
1170: -  routine_create - routine to create method context

1172:    Notes:
1173:    TSGLLERegister() may be called multiple times to add several user-defined families.

1175:    Sample usage:
1176: .vb
1177:    TSGLLERegister("my_scheme",MySchemeCreate);
1178: .ve

1180:    Then, your scheme can be chosen with the procedural interface via
1181: $     TSGLLESetType(ts,"my_scheme")
1182:    or at runtime via the option
1183: $     -ts_gl_type my_scheme

1185:    Level: advanced

1187: .seealso: `TSGLLERegisterAll()`
1188: @*/
1189: PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS))
1190: {
1191:   TSGLLEInitializePackage();
1192:   PetscFunctionListAdd(&TSGLLEList, sname, function);
1193:   return 0;
1194: }

1196: /*@C
1197:    TSGLLEAcceptRegister -  adds a TSGLLE acceptance scheme

1199:    Not Collective

1201:    Input Parameters:
1202: +  name_scheme - name of user-defined acceptance scheme
1203: -  routine_create - routine to create method context

1205:    Notes:
1206:    TSGLLEAcceptRegister() may be called multiple times to add several user-defined families.

1208:    Sample usage:
1209: .vb
1210:    TSGLLEAcceptRegister("my_scheme",MySchemeCreate);
1211: .ve

1213:    Then, your scheme can be chosen with the procedural interface via
1214: $     TSGLLESetAcceptType(ts,"my_scheme")
1215:    or at runtime via the option
1216: $     -ts_gl_accept_type my_scheme

1218:    Level: advanced

1220: .seealso: `TSGLLERegisterAll()`
1221: @*/
1222: PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFunction function)
1223: {
1224:   PetscFunctionListAdd(&TSGLLEAcceptList, sname, function);
1225:   return 0;
1226: }

1228: /*@C
1229:   TSGLLERegisterAll - Registers all of the general linear methods in TSGLLE

1231:   Not Collective

1233:   Level: advanced

1235: .seealso: `TSGLLERegisterDestroy()`
1236: @*/
1237: PetscErrorCode TSGLLERegisterAll(void)
1238: {
1239:   if (TSGLLERegisterAllCalled) return 0;
1240:   TSGLLERegisterAllCalled = PETSC_TRUE;

1242:   TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS);
1243:   TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always);
1244:   return 0;
1245: }

1247: /*@C
1248:   TSGLLEInitializePackage - This function initializes everything in the TSGLLE package. It is called
1249:   from TSInitializePackage().

1251:   Level: developer

1253: .seealso: `PetscInitialize()`
1254: @*/
1255: PetscErrorCode TSGLLEInitializePackage(void)
1256: {
1257:   if (TSGLLEPackageInitialized) return 0;
1258:   TSGLLEPackageInitialized = PETSC_TRUE;
1259:   TSGLLERegisterAll();
1260:   PetscRegisterFinalize(TSGLLEFinalizePackage);
1261:   return 0;
1262: }

1264: /*@C
1265:   TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
1266:   called from PetscFinalize().

1268:   Level: developer

1270: .seealso: `PetscFinalize()`
1271: @*/
1272: PetscErrorCode TSGLLEFinalizePackage(void)
1273: {
1274:   PetscFunctionListDestroy(&TSGLLEList);
1275:   PetscFunctionListDestroy(&TSGLLEAcceptList);
1276:   TSGLLEPackageInitialized = PETSC_FALSE;
1277:   TSGLLERegisterAllCalled  = PETSC_FALSE;
1278:   return 0;
1279: }

1281: /* ------------------------------------------------------------ */
1282: /*MC
1283:       TSGLLE - DAE solver using implicit General Linear methods

1285:   These methods contain Runge-Kutta and multistep schemes as special cases.  These special cases have some fundamental
1286:   limitations.  For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their
1287:   applicability to very stiff systems.  Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF
1288:   are not 0-stable for order greater than 6.  GL methods can be A- and L-stable with arbitrarily high stage order and
1289:   reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes.
1290:   All this is possible while preserving a singly diagonally implicit structure.

1292:   Options database keys:
1293: +  -ts_gl_type <type> - the class of general linear method (irks)
1294: .  -ts_gl_rtol <tol>  - relative error
1295: .  -ts_gl_atol <tol>  - absolute error
1296: .  -ts_gl_min_order <p> - minimum order method to consider (default=1)
1297: .  -ts_gl_max_order <p> - maximum order method to consider (default=3)
1298: .  -ts_gl_start_order <p> - order of starting method (default=1)
1299: .  -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
1300: -  -ts_adapt_type <method> - adaptive controller to use (none step both)

1302:   Notes:
1303:   This integrator can be applied to DAE.

1305:   Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK).
1306:   They are represented by the tableau

1308: .vb
1309:   A  |  U
1310:   -------
1311:   B  |  V
1312: .ve

1314:   combined with a vector c of abscissa.  "Diagonally implicit" means that A is lower triangular.
1315:   A step of the general method reads

1317: .vb
1318:   [ Y ] = [A  U] [  Y'   ]
1319:   [X^k] = [B  V] [X^{k-1}]
1320: .ve

1322:   where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of
1323:   the solution at step k.  The Nordsieck vector consists of the first r moments of the solution, given by

1325: .vb
1326:   X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
1327: .ve

1329:   If A is lower triangular, we can solve the stages (Y,Y') sequentially

1331: .vb
1332:   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}
1333: .ve

1335:   and then construct the pieces to carry to the next step

1337: .vb
1338:   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}
1339: .ve

1341:   Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i
1342:   in terms of y_i and known stuff (y_j for j<i and x_j for all j).

1344:   Error estimation

1346:   At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses
1347:   Inherent Runge-Kutta Stability (IRKS).  These methods have r=s, the number of items passed between steps is equal to
1348:   the number of stages.  The order and stage-order are one less than the number of stages.  We use the error estimates
1349:   in the 2007 paper which provide the following estimates

1351: .vb
1352:   h^{p+1} X^{(p+1)}          = phi_0^T Y' + [0 psi_0^T] Xold
1353:   h^{p+2} X^{(p+2)}          = phi_1^T Y' + [0 psi_1^T] Xold
1354:   h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold
1355: .ve

1357:   These estimates are accurate to O(h^{p+3}).

1359:   Changing the step size

1361:   We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper.

1363:   Level: beginner

1365:   References:
1366: + * - John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for
1367:   ordinary differential equations, Journal of Complexity, Vol 23, 2007.
1368: - * - John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009.

1370: .seealso: `TSCreate()`, `TS`, `TSSetType()`

1372: M*/
1373: PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
1374: {
1375:   TS_GLLE *gl;

1377:   TSGLLEInitializePackage();

1379:   PetscNew(&gl);
1380:   ts->data = (void *)gl;

1382:   ts->ops->reset          = TSReset_GLLE;
1383:   ts->ops->destroy        = TSDestroy_GLLE;
1384:   ts->ops->view           = TSView_GLLE;
1385:   ts->ops->setup          = TSSetUp_GLLE;
1386:   ts->ops->solve          = TSSolve_GLLE;
1387:   ts->ops->setfromoptions = TSSetFromOptions_GLLE;
1388:   ts->ops->snesfunction   = SNESTSFormFunction_GLLE;
1389:   ts->ops->snesjacobian   = SNESTSFormJacobian_GLLE;

1391:   ts->usessnes = PETSC_TRUE;

1393:   gl->max_step_rejections = 1;
1394:   gl->min_order           = 1;
1395:   gl->max_order           = 3;
1396:   gl->start_order         = 1;
1397:   gl->current_scheme      = -1;
1398:   gl->extrapolate         = PETSC_FALSE;

1400:   gl->wrms_atol = 1e-8;
1401:   gl->wrms_rtol = 1e-5;

1403:   PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE);
1404:   PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE);
1405:   PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE);
1406:   return 0;
1407: }