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) {
 54:       DMRestoreNamedGlobalVector(dm,"TSGLLE_Z",Z);
 55:     }
 56:   }
 57:   if (Ydotstage) {

 59:     if (dm && dm != ts->dm) {
 60:       DMRestoreNamedGlobalVector(dm,"TSGLLE_Ydot",Ydotstage);
 61:     }
 62:   }
 63:   return 0;
 64: }

 66: static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine,DM coarse,void *ctx)
 67: {
 68:   return 0;
 69: }

 71: static PetscErrorCode DMRestrictHook_TSGLLE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
 72: {
 73:   TS             ts = (TS)ctx;
 74:   Vec            Ydot,Ydot_c;

 76:   TSGLLEGetVecs(ts,fine,NULL,&Ydot);
 77:   TSGLLEGetVecs(ts,coarse,NULL,&Ydot_c);
 78:   MatRestrict(restrct,Ydot,Ydot_c);
 79:   VecPointwiseMult(Ydot_c,rscale,Ydot_c);
 80:   TSGLLERestoreVecs(ts,fine,NULL,&Ydot);
 81:   TSGLLERestoreVecs(ts,coarse,NULL,&Ydot_c);
 82:   return 0;
 83: }

 85: static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm,DM subdm,void *ctx)
 86: {
 87:   return 0;
 88: }

 90: static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm,VecScatter gscat, VecScatter lscat,DM subdm,void *ctx)
 91: {
 92:   TS             ts = (TS)ctx;
 93:   Vec            Ydot,Ydot_s;

 95:   TSGLLEGetVecs(ts,dm,NULL,&Ydot);
 96:   TSGLLEGetVecs(ts,subdm,NULL,&Ydot_s);

 98:   VecScatterBegin(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);
 99:   VecScatterEnd(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);

101:   TSGLLERestoreVecs(ts,dm,NULL,&Ydot);
102:   TSGLLERestoreVecs(ts,subdm,NULL,&Ydot_s);
103:   return 0;
104: }

106: static PetscErrorCode TSGLLESchemeCreate(PetscInt p,PetscInt q,PetscInt r,PetscInt s,const PetscScalar *c,
107:                                        const PetscScalar *a,const PetscScalar *b,const PetscScalar *u,const PetscScalar *v,TSGLLEScheme *inscheme)
108: {
109:   TSGLLEScheme     scheme;
110:   PetscInt       j;

116:   *inscheme = NULL;
117:   PetscNew(&scheme);
118:   scheme->p = p;
119:   scheme->q = q;
120:   scheme->r = r;
121:   scheme->s = s;

123:   PetscMalloc5(s,&scheme->c,s*s,&scheme->a,r*s,&scheme->b,r*s,&scheme->u,r*r,&scheme->v);
124:   PetscArraycpy(scheme->c,c,s);
125:   for (j=0; j<s*s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j];
126:   for (j=0; j<r*s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j];
127:   for (j=0; j<s*r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j];
128:   for (j=0; j<r*r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j];

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

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

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:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->beta+1,&n,&info));

163:     /* Build stage_error vector
164:            xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha;
165:     */
166:     for (i=0; i<s; i++) {
167:       scheme->stage_error[i] = CPowF(c[i],p+1);
168:       for (j=0; j<s; j++) scheme->stage_error[i] -= a[i*s+j]*CPowF(c[j],p);
169:       for (j=1; j<r; j++) scheme->stage_error[i] += u[i*r+j]*scheme->alpha[j];
170:     }

172:     /* alpha[0] (epsilon in B,J,W 2007)
173:            epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha;
174:     */
175:     scheme->alpha[0] = 1./Factorial(p+1);
176:     for (j=0; j<s; j++) scheme->alpha[0] -= b[0*s+j]*CPowF(c[j],p);
177:     for (j=1; j<r; j++) scheme->alpha[0] += v[0*r+j]*scheme->alpha[j];

179:     /* right hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */
180:     for (i=1; i<r; i++) {
181:       scheme->gamma[i] = (i==1 ? -1. : 0)*scheme->alpha[0];
182:       for (j=0; j<s; j++) scheme->gamma[i] += b[i*s+j]*scheme->stage_error[j];
183:     }
184:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->gamma+1,&n,&info));

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

197:     /* gamma[0] (sigma in B,J,W 2007)
198:     *   e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma;
199:     * */
200:     scheme->gamma[0] = 0.0;
201:     for (j=0; j<s; j++) scheme->gamma[0] += b[0*s+j]*scheme->stage_error[j];
202:     for (j=1; j<r; j++) scheme->gamma[0] += v[0*s+j]*scheme->gamma[j];

204:     /* Assemble H
205:     *    % Determine 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.;  bmat[0+1*ss] = 0.;  bmat[0+2*ss] = 0.;
226:     bmat[1+0*ss] = 1.;  bmat[1+1*ss] = 1.;  bmat[1+2*ss] = 0.;
227:     bmat[2+0*ss] = 0.;  bmat[2+1*ss] = 0.;  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:     PetscStackCallBLAS("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:     PetscStackCallBLAS("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++) if (a[(s-1)*s+j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE;
266:   for (j=0; j<r; j++) if (u[(s-1)*r+j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE;
267:   scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */
268:   for (j=0; j<s-1; j++) if (r>1 && b[1*s+j] != 0.) scheme->fsal = PETSC_FALSE;
269:   if (b[1*s+r-1] != 1.) scheme->fsal = PETSC_FALSE;
270:   for (j=0; j<r; j++) if (r>1 && v[1*r+j] != 0.) scheme->fsal = PETSC_FALSE;

272:   *inscheme = scheme;
273:   return 0;
274: }

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

284: static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
285: {
286:   PetscInt       i;

288:   for (i=0; i<gl->nschemes; i++) {
289:     if (gl->schemes[i]) TSGLLESchemeDestroy(gl->schemes[i]);
290:   }
291:   PetscFree(gl->schemes);
292:   gl->nschemes = 0;
293:   PetscMemzero(gl->type_name,sizeof(gl->type_name));
294:   return 0;
295: }

297: static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer,PetscInt m,PetscInt n,const PetscScalar a[],const char name[])
298: {
299:   PetscBool      iascii;
300:   PetscInt       i,j;

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

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

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

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

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

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

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

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

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

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

422: static PetscErrorCode TSGLLECreate_IRKS(TS ts)
423: {
424:   TS_GLLE        *gl = (TS_GLLE*)ts->data;

426:   gl->Destroy               = TSGLLEDestroy_Default;
427:   gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default;
428:   gl->CompleteStep          = TSGLLECompleteStep_RescaleAndModify;
429:   PetscMalloc1(10,&gl->schemes);
430:   gl->nschemes = 0;

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

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

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

558:    Collective on TS

560:    Input Parameters:
561: +  ts - the TS context
562: -  type - a method

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

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

571:    Normally, it is best to use the TSSetFromOptions() command and
572:    then set the TSGLLE type from the options database rather than by using
573:    this routine.  Using the options database provides the user with
574:    maximum flexibility in evaluating the many different solvers.
575:    The TSGLLESetType() routine is provided for those situations where it
576:    is necessary to set the timestepping solver independently of the
577:    command line or options database.  This might be the case, for example,
578:    when the choice of solver changes during the execution of the
579:    program, and the user's application is taking responsibility for
580:    choosing the appropriate method.

582:    Level: intermediate

584: @*/
585: PetscErrorCode  TSGLLESetType(TS ts,TSGLLEType type)
586: {
589:   PetscTryMethod(ts,"TSGLLESetType_C",(TS,TSGLLEType),(ts,type));
590:   return 0;
591: }

593: /*@C
594:    TSGLLESetAcceptType - sets the acceptance test

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

599:    Logically Collective on TS

601:    Input Parameters:
602: +  ts - the TS context
603: -  type - the type

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

608:    Level: intermediate

610: .seealso: TS, TSGLLE, TSGLLEAcceptRegister(), TSGLLEAdapt, set type
611: @*/
612: PetscErrorCode  TSGLLESetAcceptType(TS ts,TSGLLEAcceptType type)
613: {
616:   PetscTryMethod(ts,"TSGLLESetAcceptType_C",(TS,TSGLLEAcceptType),(ts,type));
617:   return 0;
618: }

620: /*@C
621:    TSGLLEGetAdapt - gets the TSGLLEAdapt object from the TS

623:    Not Collective

625:    Input Parameter:
626: .  ts - the TS context

628:    Output Parameter:
629: .  adapt - the TSGLLEAdapt context

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

635:    Level: advanced

637: .seealso: TSGLLEAdapt, TSGLLEAdaptRegister()
638: @*/
639: PetscErrorCode  TSGLLEGetAdapt(TS ts,TSGLLEAdapt *adapt)
640: {
643:   PetscUseMethod(ts,"TSGLLEGetAdapt_C",(TS,TSGLLEAdapt*),(ts,adapt));
644:   return 0;
645: }

647: static PetscErrorCode TSGLLEAccept_Always(TS ts,PetscReal tleft,PetscReal h,const PetscReal enorms[],PetscBool  *accept)
648: {
649:   *accept = PETSC_TRUE;
650:   return 0;
651: }

653: static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
654: {
655:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
656:   PetscScalar    *x,*w;
657:   PetscInt       n,i;

659:   VecGetArray(gl->X[0],&x);
660:   VecGetArray(gl->W,&w);
661:   VecGetLocalSize(gl->W,&n);
662:   for (i=0; i<n; i++) w[i] = 1./(gl->wrms_atol + gl->wrms_rtol*PetscAbsScalar(x[i]));
663:   VecRestoreArray(gl->X[0],&x);
664:   VecRestoreArray(gl->W,&w);
665:   return 0;
666: }

668: static PetscErrorCode TSGLLEVecNormWRMS(TS ts,Vec X,PetscReal *nrm)
669: {
670:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
671:   PetscScalar    *x,*w;
672:   PetscReal      sum = 0.0,gsum;
673:   PetscInt       n,N,i;

675:   VecGetArray(X,&x);
676:   VecGetArray(gl->W,&w);
677:   VecGetLocalSize(gl->W,&n);
678:   for (i=0; i<n; i++) sum += PetscAbsScalar(PetscSqr(x[i]*w[i]));
679:   VecRestoreArray(X,&x);
680:   VecRestoreArray(gl->W,&w);
681:   MPIU_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
682:   VecGetSize(gl->W,&N);
683:   *nrm = PetscSqrtReal(gsum/(1.*N));
684:   return 0;
685: }

687: static PetscErrorCode TSGLLESetType_GLLE(TS ts,TSGLLEType type)
688: {
689:   PetscBool      same;
690:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
691:   PetscErrorCode (*r)(TS);

693:   if (gl->type_name[0]) {
694:     PetscStrcmp(gl->type_name,type,&same);
695:     if (same) return 0;
696:     (*gl->Destroy)(gl);
697:   }

699:   PetscFunctionListFind(TSGLLEList,type,&r);
701:   (*r)(ts);
702:   PetscStrcpy(gl->type_name,type);
703:   return 0;
704: }

706: static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts,TSGLLEAcceptType type)
707: {
708:   TSGLLEAcceptFunction r;
709:   TS_GLLE              *gl = (TS_GLLE*)ts->data;

711:   PetscFunctionListFind(TSGLLEAcceptList,type,&r);
713:   gl->Accept = r;
714:   PetscStrncpy(gl->accept_name,type,sizeof(gl->accept_name));
715:   return 0;
716: }

718: static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts,TSGLLEAdapt *adapt)
719: {
720:   TS_GLLE        *gl = (TS_GLLE*)ts->data;

722:   if (!gl->adapt) {
723:     TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts),&gl->adapt);
724:     PetscObjectIncrementTabLevel((PetscObject)gl->adapt,(PetscObject)ts,1);
725:     PetscLogObjectParent((PetscObject)ts,(PetscObject)gl->adapt);
726:   }
727:   *adapt = gl->adapt;
728:   return 0;
729: }

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

737:   cur   = -1;
738:   cur_p = gl->schemes[gl->current_scheme]->p;
739:   tleft = ts->max_time - (ts->ptime + ts->time_step);
740:   for (i=0,n=0; i<gl->nschemes; i++) {
741:     TSGLLEScheme sc = gl->schemes[i];
742:     if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
743:     if (sc->p == cur_p - 1)    errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[0];
744:     else if (sc->p == cur_p)   errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[1];
745:     else if (sc->p == cur_p+1) errors[n] = PetscAbsScalar(sc->alpha[0])*(hmnorm[2]+hmnorm[3]);
746:     else continue;
747:     candidates[n] = i;
748:     orders[n]     = PetscMin(sc->p,sc->q); /* order of global truncation error */
749:     costs[n]      = sc->s;                 /* estimate the cost as the number of stages */
750:     if (i == gl->current_scheme) cur = n;
751:     n++;
752:   }
754:   TSGLLEAdaptChoose(gl->adapt,n,orders,errors,costs,cur,h,tleft,&next_sc,next_h,finish);
755:   *next_scheme = candidates[next_sc];
756:   PetscInfo(ts,"Adapt chose scheme %d (%d,%d,%d,%d) with step size %6.2e, finish=%d\n",*next_scheme,gl->schemes[*next_scheme]->p,gl->schemes[*next_scheme]->q,gl->schemes[*next_scheme]->r,gl->schemes[*next_scheme]->s,*next_h,*finish);
757:   return 0;
758: }

760: static PetscErrorCode TSGLLEGetMaxSizes(TS ts,PetscInt *max_r,PetscInt *max_s)
761: {
762:   TS_GLLE *gl = (TS_GLLE*)ts->data;

764:   *max_r = gl->schemes[gl->nschemes-1]->r;
765:   *max_s = gl->schemes[gl->nschemes-1]->s;
766:   return 0;
767: }

769: static PetscErrorCode TSSolve_GLLE(TS ts)
770: {
771:   TS_GLLE             *gl = (TS_GLLE*)ts->data;
772:   PetscInt            i,k,its,lits,max_r,max_s;
773:   PetscBool           final_step,finish;
774:   SNESConvergedReason snesreason;

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

778:   TSGLLEGetMaxSizes(ts,&max_r,&max_s);
779:   VecCopy(ts->vec_sol,gl->X[0]);
780:   for (i=1; i<max_r; i++) {
781:     VecZeroEntries(gl->X[i]);
782:   }
783:   TSGLLEUpdateWRMS(ts);

785:   if (0) {
786:     /* Find consistent initial data for DAE */
787:     gl->stage_time = ts->ptime + ts->time_step;
788:     gl->scoeff = 1.;
789:     gl->stage  = 0;

791:     VecCopy(ts->vec_sol,gl->Z);
792:     VecCopy(ts->vec_sol,gl->Y);
793:     SNESSolve(ts->snes,NULL,gl->Y);
794:     SNESGetIterationNumber(ts->snes,&its);
795:     SNESGetLinearSolveIterations(ts->snes,&lits);
796:     SNESGetConvergedReason(ts->snes,&snesreason);

798:     ts->snes_its += its; ts->ksp_its += lits;
799:     if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
800:       ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
801:       PetscInfo(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
802:       return 0;
803:     }
804:   }


808:   for (k=0,final_step=PETSC_FALSE,finish=PETSC_FALSE; k<ts->max_steps && !finish; k++) {
809:     PetscInt          j,r,s,next_scheme = 0,rejections;
810:     PetscReal         h,hmnorm[4],enorm[3],next_h;
811:     PetscBool         accept;
812:     const PetscScalar *c,*a,*u;
813:     Vec               *X,*Ydot,Y;
814:     TSGLLEScheme        scheme = gl->schemes[gl->current_scheme];

816:     r = scheme->r; s = scheme->s;
817:     c = scheme->c;
818:     a = scheme->a; u = scheme->u;
819:     h = ts->time_step;
820:     X = gl->X; Ydot = gl->Ydot; Y = gl->Y;

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

824:     /*
825:       We only call PreStep at the start of each STEP, not each STAGE.  This is because it is
826:       possible to fail (have to restart a step) after multiple stages.
827:     */
828:     TSPreStep(ts);

830:     rejections = 0;
831:     while (1) {
832:       for (i=0; i<s; i++) {
833:         PetscScalar shift;
834:         gl->scoeff     = 1./PetscRealPart(a[i*s+i]);
835:         shift          = gl->scoeff/ts->time_step;
836:         gl->stage      = i;
837:         gl->stage_time = ts->ptime + PetscRealPart(c[i])*h;

839:         /*
840:         * Stage equation: Y = h A Y' + U X
841:         * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
842:         * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
843:         * Then y'_i = z + 1/(h a_ii) y_i
844:         */
845:         VecZeroEntries(gl->Z);
846:         for (j=0; j<r; j++) {
847:           VecAXPY(gl->Z,-shift*u[i*r+j],X[j]);
848:         }
849:         for (j=0; j<i; j++) {
850:           VecAXPY(gl->Z,-shift*h*a[i*s+j],Ydot[j]);
851:         }
852:         /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */

854:         /* Compute an estimate of Y to start Newton iteration */
855:         if (gl->extrapolate) {
856:           if (i==0) {
857:             /* Linear extrapolation on the first stage */
858:             VecWAXPY(Y,c[i]*h,X[1],X[0]);
859:           } else {
860:             /* Linear extrapolation from the last stage */
861:             VecAXPY(Y,(c[i]-c[i-1])*h,Ydot[i-1]);
862:           }
863:         } else if (i==0) {        /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
864:           VecCopy(X[0],Y);
865:         }

867:         /* Solve this stage (Ydot[i] is computed during function evaluation) */
868:         SNESSolve(ts->snes,NULL,Y);
869:         SNESGetIterationNumber(ts->snes,&its);
870:         SNESGetLinearSolveIterations(ts->snes,&lits);
871:         SNESGetConvergedReason(ts->snes,&snesreason);
872:         ts->snes_its += its; ts->ksp_its += lits;
873:         if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
874:           ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
875:           PetscInfo(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
876:           return 0;
877:         }
878:       }

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

882:       (*gl->EstimateHigherMoments)(scheme,h,Ydot,gl->X,gl->himom);
883:       /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */
884:       for (i=0; i<3; i++) {
885:         TSGLLEVecNormWRMS(ts,gl->himom[i],&hmnorm[i+1]);
886:       }
887:       enorm[0] = PetscRealPart(scheme->alpha[0])*hmnorm[1];
888:       enorm[1] = PetscRealPart(scheme->beta[0]) *hmnorm[2];
889:       enorm[2] = PetscRealPart(scheme->gamma[0])*hmnorm[3];
890:       (*gl->Accept)(ts,ts->max_time-gl->stage_time,h,enorm,&accept);
891:       if (accept) goto accepted;
892:       rejections++;
893:       PetscInfo(ts,"Step %D (t=%g) not accepted, rejections=%D\n",k,gl->stage_time,rejections);
894:       if (rejections > gl->max_step_rejections) break;
895:       /*
896:         There are lots of reasons why a step might be rejected, including solvers not converging and other factors that
897:         TSGLLEChooseNextScheme does not support.  Additionally, the error estimates may be very screwed up, so I'm not
898:         convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor
899:         (the adaptor interface probably has to change).  Here we make an arbitrary and naive choice.  This assumes that
900:         steps were written in Nordsieck form.  The "correct" method would be to re-complete the previous time step with
901:         the correct "next" step size.  It is unclear to me whether the present ad-hoc method of rescaling X is stable.
902:       */
903:       h *= 0.5;
904:       for (i=1; i<scheme->r; i++) {
905:         VecScale(X[i],PetscPowRealInt(0.5,i));
906:       }
907:     }
908:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Time step %D (t=%g) not accepted after %D failures",k,gl->stage_time,rejections);

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

915:     PetscInfo(ts,"Last moment norm %10.2e, estimated error norms %10.2e %10.2e %10.2e\n",hmnorm[0],enorm[0],enorm[1],enorm[2]);
916:     if (!final_step) {
917:       TSGLLEChooseNextScheme(ts,h,hmnorm,&next_scheme,&next_h,&final_step);
918:     } else {
919:       /* Dummy values to complete the current step in a consistent manner */
920:       next_scheme = gl->current_scheme;
921:       next_h      = h;
922:       finish      = PETSC_TRUE;
923:     }

925:     X        = gl->Xold;
926:     gl->Xold = gl->X;
927:     gl->X    = X;
928:     (*gl->CompleteStep)(scheme,h,gl->schemes[next_scheme],next_h,Ydot,gl->Xold,gl->X);

930:     TSGLLEUpdateWRMS(ts);

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

937:     TSPostEvaluate(ts);
938:     TSPostStep(ts);
939:     TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);

941:     gl->current_scheme = next_scheme;
942:     ts->time_step      = next_h;
943:   }
944:   return 0;
945: }

947: /*------------------------------------------------------------*/

949: static PetscErrorCode TSReset_GLLE(TS ts)
950: {
951:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
952:   PetscInt       max_r,max_s;

954:   if (gl->setupcalled) {
955:     TSGLLEGetMaxSizes(ts,&max_r,&max_s);
956:     VecDestroyVecs(max_r,&gl->Xold);
957:     VecDestroyVecs(max_r,&gl->X);
958:     VecDestroyVecs(max_s,&gl->Ydot);
959:     VecDestroyVecs(3,&gl->himom);
960:     VecDestroy(&gl->W);
961:     VecDestroy(&gl->Y);
962:     VecDestroy(&gl->Z);
963:   }
964:   gl->setupcalled = PETSC_FALSE;
965:   return 0;
966: }

968: static PetscErrorCode TSDestroy_GLLE(TS ts)
969: {
970:   TS_GLLE        *gl = (TS_GLLE*)ts->data;

972:   TSReset_GLLE(ts);
973:   if (ts->dm) {
974:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);
975:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);
976:   }
977:   if (gl->adapt) TSGLLEAdaptDestroy(&gl->adapt);
978:   if (gl->Destroy) (*gl->Destroy)(gl);
979:   PetscFree(ts->data);
980:   PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C",NULL);
981:   PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",NULL);
982:   PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C",NULL);
983:   return 0;
984: }

986: /*
987:     This defines the nonlinear equation that is to be solved with SNES
988:     g(x) = f(t,x,z+shift*x) = 0
989: */
990: static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes,Vec x,Vec f,TS ts)
991: {
992:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
993:   Vec            Z,Ydot;
994:   DM             dm,dmsave;

996:   SNESGetDM(snes,&dm);
997:   TSGLLEGetVecs(ts,dm,&Z,&Ydot);
998:   VecWAXPY(Ydot,gl->scoeff/ts->time_step,x,Z);
999:   dmsave = ts->dm;
1000:   ts->dm = dm;
1001:   TSComputeIFunction(ts,gl->stage_time,x,Ydot,f,PETSC_FALSE);
1002:   ts->dm = dmsave;
1003:   TSGLLERestoreVecs(ts,dm,&Z,&Ydot);
1004:   return 0;
1005: }

1007: static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes,Vec x,Mat A,Mat B,TS ts)
1008: {
1009:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
1010:   Vec            Z,Ydot;
1011:   DM             dm,dmsave;

1013:   SNESGetDM(snes,&dm);
1014:   TSGLLEGetVecs(ts,dm,&Z,&Ydot);
1015:   dmsave = ts->dm;
1016:   ts->dm = dm;
1017:   /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
1018:   TSComputeIJacobian(ts,gl->stage_time,x,gl->Ydot[gl->stage],gl->scoeff/ts->time_step,A,B,PETSC_FALSE);
1019:   ts->dm = dmsave;
1020:   TSGLLERestoreVecs(ts,dm,&Z,&Ydot);
1021:   return 0;
1022: }

1024: static PetscErrorCode TSSetUp_GLLE(TS ts)
1025: {
1026:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
1027:   PetscInt       max_r,max_s;
1028:   DM             dm;

1030:   gl->setupcalled = PETSC_TRUE;
1031:   TSGLLEGetMaxSizes(ts,&max_r,&max_s);
1032:   VecDuplicateVecs(ts->vec_sol,max_r,&gl->X);
1033:   VecDuplicateVecs(ts->vec_sol,max_r,&gl->Xold);
1034:   VecDuplicateVecs(ts->vec_sol,max_s,&gl->Ydot);
1035:   VecDuplicateVecs(ts->vec_sol,3,&gl->himom);
1036:   VecDuplicate(ts->vec_sol,&gl->W);
1037:   VecDuplicate(ts->vec_sol,&gl->Y);
1038:   VecDuplicate(ts->vec_sol,&gl->Z);

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

1044:   if (gl->current_scheme < 0) {
1045:     PetscInt i;
1046:     for (i=0;; i++) {
1047:       if (gl->schemes[i]->p == gl->start_order) break;
1049:     }
1050:     gl->current_scheme = i;
1051:   }
1052:   TSGetDM(ts,&dm);
1053:   DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);
1054:   DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);
1055:   return 0;
1056: }
1057: /*------------------------------------------------------------*/

1059: static PetscErrorCode TSSetFromOptions_GLLE(PetscOptionItems *PetscOptionsObject,TS ts)
1060: {
1061:   TS_GLLE        *gl        = (TS_GLLE*)ts->data;
1062:   char           tname[256] = TSGLLE_IRKS,completef[256] = "rescale-and-modify";

1064:   PetscOptionsHead(PetscOptionsObject,"General Linear ODE solver options");
1065:   {
1066:     PetscBool flg;
1067:     PetscOptionsFList("-ts_gl_type","Type of GL method","TSGLLESetType",TSGLLEList,gl->type_name[0] ? gl->type_name : tname,tname,sizeof(tname),&flg);
1068:     if (flg || !gl->type_name[0]) {
1069:       TSGLLESetType(ts,tname);
1070:     }
1071:     PetscOptionsInt("-ts_gl_max_step_rejections","Maximum number of times to attempt a step","None",gl->max_step_rejections,&gl->max_step_rejections,NULL);
1072:     PetscOptionsInt("-ts_gl_max_order","Maximum order to try","TSGLLESetMaxOrder",gl->max_order,&gl->max_order,NULL);
1073:     PetscOptionsInt("-ts_gl_min_order","Minimum order to try","TSGLLESetMinOrder",gl->min_order,&gl->min_order,NULL);
1074:     PetscOptionsInt("-ts_gl_start_order","Initial order to try","TSGLLESetMinOrder",gl->start_order,&gl->start_order,NULL);
1075:     PetscOptionsEnum("-ts_gl_error_direction","Which direction to look when estimating error","TSGLLESetErrorDirection",TSGLLEErrorDirections,(PetscEnum)gl->error_direction,(PetscEnum*)&gl->error_direction,NULL);
1076:     PetscOptionsBool("-ts_gl_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSGLLESetExtrapolate",gl->extrapolate,&gl->extrapolate,NULL);
1077:     PetscOptionsReal("-ts_gl_atol","Absolute tolerance","TSGLLESetTolerances",gl->wrms_atol,&gl->wrms_atol,NULL);
1078:     PetscOptionsReal("-ts_gl_rtol","Relative tolerance","TSGLLESetTolerances",gl->wrms_rtol,&gl->wrms_rtol,NULL);
1079:     PetscOptionsString("-ts_gl_complete","Method to use for completing the step","none",completef,completef,sizeof(completef),&flg);
1080:     if (flg) {
1081:       PetscBool match1,match2;
1082:       PetscStrcmp(completef,"rescale",&match1);
1083:       PetscStrcmp(completef,"rescale-and-modify",&match2);
1084:       if (match1)      gl->CompleteStep = TSGLLECompleteStep_Rescale;
1085:       else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
1086:       else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"%s",completef);
1087:     }
1088:     {
1089:       char type[256] = TSGLLEACCEPT_ALWAYS;
1090:       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);
1091:       if (flg || !gl->accept_name[0]) {
1092:         TSGLLESetAcceptType(ts,type);
1093:       }
1094:     }
1095:     {
1096:       TSGLLEAdapt adapt;
1097:       TSGLLEGetAdapt(ts,&adapt);
1098:       TSGLLEAdaptSetFromOptions(PetscOptionsObject,adapt);
1099:     }
1100:   }
1101:   PetscOptionsTail();
1102:   return 0;
1103: }

1105: static PetscErrorCode TSView_GLLE(TS ts,PetscViewer viewer)
1106: {
1107:   TS_GLLE        *gl = (TS_GLLE*)ts->data;
1108:   PetscInt       i;
1109:   PetscBool      iascii,details;

1111:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1112:   if (iascii) {
1113:     PetscViewerASCIIPrintf(viewer,"  min order %D, max order %D, current order %D\n",gl->min_order,gl->max_order,gl->schemes[gl->current_scheme]->p);
1114:     PetscViewerASCIIPrintf(viewer,"  Error estimation: %s\n",TSGLLEErrorDirections[gl->error_direction]);
1115:     PetscViewerASCIIPrintf(viewer,"  Extrapolation: %s\n",gl->extrapolate ? "yes" : "no");
1116:     PetscViewerASCIIPrintf(viewer,"  Acceptance test: %s\n",gl->accept_name[0] ? gl->accept_name : "(not yet set)");
1117:     PetscViewerASCIIPushTab(viewer);
1118:     TSGLLEAdaptView(gl->adapt,viewer);
1119:     PetscViewerASCIIPopTab(viewer);
1120:     PetscViewerASCIIPrintf(viewer,"  type: %s\n",gl->type_name[0] ? gl->type_name : "(not yet set)");
1121:     PetscViewerASCIIPrintf(viewer,"Schemes within family (%d):\n",gl->nschemes);
1122:     details = PETSC_FALSE;
1123:     PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_gl_view_detailed",&details,NULL);
1124:     PetscViewerASCIIPushTab(viewer);
1125:     for (i=0; i<gl->nschemes; i++) {
1126:       TSGLLESchemeView(gl->schemes[i],details,viewer);
1127:     }
1128:     if (gl->View) {
1129:       (*gl->View)(gl,viewer);
1130:     }
1131:     PetscViewerASCIIPopTab(viewer);
1132:   }
1133:   return 0;
1134: }

1136: /*@C
1137:    TSGLLERegister -  adds a TSGLLE implementation

1139:    Not Collective

1141:    Input Parameters:
1142: +  name_scheme - name of user-defined general linear scheme
1143: -  routine_create - routine to create method context

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

1148:    Sample usage:
1149: .vb
1150:    TSGLLERegister("my_scheme",MySchemeCreate);
1151: .ve

1153:    Then, your scheme can be chosen with the procedural interface via
1154: $     TSGLLESetType(ts,"my_scheme")
1155:    or at runtime via the option
1156: $     -ts_gl_type my_scheme

1158:    Level: advanced

1160: .seealso: TSGLLERegisterAll()
1161: @*/
1162: PetscErrorCode  TSGLLERegister(const char sname[],PetscErrorCode (*function)(TS))
1163: {
1164:   TSGLLEInitializePackage();
1165:   PetscFunctionListAdd(&TSGLLEList,sname,function);
1166:   return 0;
1167: }

1169: /*@C
1170:    TSGLLEAcceptRegister -  adds a TSGLLE acceptance scheme

1172:    Not Collective

1174:    Input Parameters:
1175: +  name_scheme - name of user-defined acceptance scheme
1176: -  routine_create - routine to create method context

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

1181:    Sample usage:
1182: .vb
1183:    TSGLLEAcceptRegister("my_scheme",MySchemeCreate);
1184: .ve

1186:    Then, your scheme can be chosen with the procedural interface via
1187: $     TSGLLESetAcceptType(ts,"my_scheme")
1188:    or at runtime via the option
1189: $     -ts_gl_accept_type my_scheme

1191:    Level: advanced

1193: .seealso: TSGLLERegisterAll()
1194: @*/
1195: PetscErrorCode  TSGLLEAcceptRegister(const char sname[],TSGLLEAcceptFunction function)
1196: {
1197:   PetscFunctionListAdd(&TSGLLEAcceptList,sname,function);
1198:   return 0;
1199: }

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

1204:   Not Collective

1206:   Level: advanced

1208: .seealso:  TSGLLERegisterDestroy()
1209: @*/
1210: PetscErrorCode  TSGLLERegisterAll(void)
1211: {
1212:   if (TSGLLERegisterAllCalled) return 0;
1213:   TSGLLERegisterAllCalled = PETSC_TRUE;

1215:   TSGLLERegister(TSGLLE_IRKS,              TSGLLECreate_IRKS);
1216:   TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS,TSGLLEAccept_Always);
1217:   return 0;
1218: }

1220: /*@C
1221:   TSGLLEInitializePackage - This function initializes everything in the TSGLLE package. It is called
1222:   from TSInitializePackage().

1224:   Level: developer

1226: .seealso: PetscInitialize()
1227: @*/
1228: PetscErrorCode  TSGLLEInitializePackage(void)
1229: {
1230:   if (TSGLLEPackageInitialized) return 0;
1231:   TSGLLEPackageInitialized = PETSC_TRUE;
1232:   TSGLLERegisterAll();
1233:   PetscRegisterFinalize(TSGLLEFinalizePackage);
1234:   return 0;
1235: }

1237: /*@C
1238:   TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
1239:   called from PetscFinalize().

1241:   Level: developer

1243: .seealso: PetscFinalize()
1244: @*/
1245: PetscErrorCode  TSGLLEFinalizePackage(void)
1246: {
1247:   PetscFunctionListDestroy(&TSGLLEList);
1248:   PetscFunctionListDestroy(&TSGLLEAcceptList);
1249:   TSGLLEPackageInitialized = PETSC_FALSE;
1250:   TSGLLERegisterAllCalled  = PETSC_FALSE;
1251:   return 0;
1252: }

1254: /* ------------------------------------------------------------ */
1255: /*MC
1256:       TSGLLE - DAE solver using implicit General Linear methods

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

1265:   Options database keys:
1266: +  -ts_gl_type <type> - the class of general linear method (irks)
1267: .  -ts_gl_rtol <tol>  - relative error
1268: .  -ts_gl_atol <tol>  - absolute error
1269: .  -ts_gl_min_order <p> - minimum order method to consider (default=1)
1270: .  -ts_gl_max_order <p> - maximum order method to consider (default=3)
1271: .  -ts_gl_start_order <p> - order of starting method (default=1)
1272: .  -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
1273: -  -ts_adapt_type <method> - adaptive controller to use (none step both)

1275:   Notes:
1276:   This integrator can be applied to DAE.

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

1281: .vb
1282:   A  |  U
1283:   -------
1284:   B  |  V
1285: .ve

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

1290: .vb
1291:   [ Y ] = [A  U] [  Y'   ]
1292:   [X^k] = [B  V] [X^{k-1}]
1293: .ve

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

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

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

1304: .vb
1305:   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}
1306: .ve

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

1310: .vb
1311:   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}
1312: .ve

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

1317:   Error estimation

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

1324: .vb
1325:   h^{p+1} X^{(p+1)}          = phi_0^T Y' + [0 psi_0^T] Xold
1326:   h^{p+2} X^{(p+2)}          = phi_1^T Y' + [0 psi_1^T] Xold
1327:   h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold
1328: .ve

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

1332:   Changing the step size

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

1336:   Level: beginner

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

1343: .seealso:  TSCreate(), TS, TSSetType()

1345: M*/
1346: PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
1347: {
1348:   TS_GLLE        *gl;

1350:   TSGLLEInitializePackage();

1352:   PetscNewLog(ts,&gl);
1353:   ts->data = (void*)gl;

1355:   ts->ops->reset          = TSReset_GLLE;
1356:   ts->ops->destroy        = TSDestroy_GLLE;
1357:   ts->ops->view           = TSView_GLLE;
1358:   ts->ops->setup          = TSSetUp_GLLE;
1359:   ts->ops->solve          = TSSolve_GLLE;
1360:   ts->ops->setfromoptions = TSSetFromOptions_GLLE;
1361:   ts->ops->snesfunction   = SNESTSFormFunction_GLLE;
1362:   ts->ops->snesjacobian   = SNESTSFormJacobian_GLLE;

1364:   ts->usessnes = PETSC_TRUE;

1366:   gl->max_step_rejections = 1;
1367:   gl->min_order           = 1;
1368:   gl->max_order           = 3;
1369:   gl->start_order         = 1;
1370:   gl->current_scheme      = -1;
1371:   gl->extrapolate         = PETSC_FALSE;

1373:   gl->wrms_atol = 1e-8;
1374:   gl->wrms_rtol = 1e-5;

1376:   PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C",      &TSGLLESetType_GLLE);
1377:   PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",&TSGLLESetAcceptType_GLLE);
1378:   PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C",     &TSGLLEGetAdapt_GLLE);
1379:   return 0;
1380: }