Actual source code: glee.c

  1: /*
  2:   Code for time stepping with the General Linear with Error Estimation method

  4:   Notes:
  5:   The general system is written as

  7:   Udot = F(t,U)

  9: */
 10: #include <petsc/private/tsimpl.h>
 11: #include <petscdm.h>

 13: static PetscBool  cited = PETSC_FALSE;
 14: static const char citation[] =
 15:   "@ARTICLE{Constantinescu_TR2016b,\n"
 16:   " author = {Constantinescu, E.M.},\n"
 17:   " title = {Estimating Global Errors in Time Stepping},\n"
 18:   " journal = {ArXiv e-prints},\n"
 19:   " year = 2016,\n"
 20:   " adsurl = {http://adsabs.harvard.edu/abs/2015arXiv150305166C}\n}\n";

 22: static TSGLEEType   TSGLEEDefaultType = TSGLEE35;
 23: static PetscBool    TSGLEERegisterAllCalled;
 24: static PetscBool    TSGLEEPackageInitialized;
 25: static PetscInt     explicit_stage_time_id;

 27: typedef struct _GLEETableau *GLEETableau;
 28: struct _GLEETableau {
 29:   char      *name;
 30:   PetscInt   order;               /* Classical approximation order of the method i*/
 31:   PetscInt   s;                   /* Number of stages */
 32:   PetscInt   r;                   /* Number of steps */
 33:   PetscReal  gamma;               /* LTE ratio */
 34:   PetscReal *A,*B,*U,*V,*S,*F,*c; /* Tableau */
 35:   PetscReal *Fembed;              /* Embedded final method coefficients */
 36:   PetscReal *Ferror;              /* Coefficients for computing error   */
 37:   PetscReal *Serror;              /* Coefficients for initializing the error   */
 38:   PetscInt   pinterp;             /* Interpolation order */
 39:   PetscReal *binterp;             /* Interpolation coefficients */
 40:   PetscReal  ccfl;                /* Placeholder for CFL coefficient relative to forward Euler  */
 41: };
 42: typedef struct _GLEETableauLink *GLEETableauLink;
 43: struct _GLEETableauLink {
 44:   struct _GLEETableau tab;
 45:   GLEETableauLink     next;
 46: };
 47: static GLEETableauLink GLEETableauList;

 49: typedef struct {
 50:   GLEETableau  tableau;
 51:   Vec          *Y;         /* Solution vector (along with auxiliary solution y~ or eps) */
 52:   Vec          *X;         /* Temporary solution vector */
 53:   Vec          *YStage;    /* Stage values */
 54:   Vec          *YdotStage; /* Stage right hand side */
 55:   Vec          W;          /* Right-hand-side for implicit stage solve */
 56:   Vec          Ydot;       /* Work vector holding Ydot during residual evaluation */
 57:   Vec          yGErr;      /* Vector holding the global error after a step is completed */
 58:   PetscScalar  *swork;     /* Scalar work (size of the number of stages)*/
 59:   PetscScalar  *rwork;     /* Scalar work (size of the number of steps)*/
 60:   PetscReal    scoeff;     /* shift = scoeff/dt */
 61:   PetscReal    stage_time;
 62:   TSStepStatus status;
 63: } TS_GLEE;

 65: /*MC
 66:      TSGLEE23 - Second order three stage GLEE method

 68:      This method has three stages.
 69:      s = 3, r = 2

 71:      Level: advanced

 73: .seealso: TSGLEE
 74: M*/
 75: /*MC
 76:      TSGLEE24 - Second order four stage GLEE method

 78:      This method has four stages.
 79:      s = 4, r = 2

 81:      Level: advanced

 83: .seealso: TSGLEE
 84: M*/
 85: /*MC
 86:      TSGLEE25i - Second order five stage GLEE method

 88:      This method has five stages.
 89:      s = 5, r = 2

 91:      Level: advanced

 93: .seealso: TSGLEE
 94: M*/
 95: /*MC
 96:      TSGLEE35  - Third order five stage GLEE method

 98:      This method has five stages.
 99:      s = 5, r = 2

101:      Level: advanced

103: .seealso: TSGLEE
104: M*/
105: /*MC
106:      TSGLEEEXRK2A  - Second order six stage GLEE method

108:      This method has six stages.
109:      s = 6, r = 2

111:      Level: advanced

113: .seealso: TSGLEE
114: M*/
115: /*MC
116:      TSGLEERK32G1  - Third order eight stage GLEE method

118:      This method has eight stages.
119:      s = 8, r = 2

121:      Level: advanced

123: .seealso: TSGLEE
124: M*/
125: /*MC
126:      TSGLEERK285EX  - Second order nine stage GLEE method

128:      This method has nine stages.
129:      s = 9, r = 2

131:      Level: advanced

133: .seealso: TSGLEE
134: M*/

136: /*@C
137:   TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in TSGLEE

139:   Not Collective, but should be called by all processes which will need the schemes to be registered

141:   Level: advanced

143: .seealso:  TSGLEERegisterDestroy()
144: @*/
145: PetscErrorCode TSGLEERegisterAll(void)
146: {

150:   if (TSGLEERegisterAllCalled) return(0);
151:   TSGLEERegisterAllCalled = PETSC_TRUE;

153:   {
154: #define GAMMA 0.5
155:     /* y-eps form */
156:     const PetscInt
157:       p = 1,
158:       s = 3,
159:       r = 2;
160:     const PetscReal
161:       A[3][3]   = {{1.0,0,0},{0,0.5,0},{0,0.5,0.5}},
162:       B[2][3]   = {{1.0,0,0},{-2.0,1.0,1.0}},
163:       U[3][2]   = {{1.0,0},{1.0,0.5},{1.0,0.5}},
164:       V[2][2]   = {{1,0},{0,1}},
165:       S[2]      = {1,0},
166:       F[2]      = {1,0},
167:       Fembed[2] = {1,1-GAMMA},
168:       Ferror[2] = {0,1},
169:       Serror[2] = {1,0};
170:     TSGLEERegister(TSGLEEi1,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
171:   }
172:   {
173: #undef GAMMA
174: #define GAMMA 0.0
175:     /* y-eps form */
176:     const PetscInt
177:       p = 2,
178:       s = 3,
179:       r = 2;
180:     const PetscReal
181:       A[3][3]   = {{0,0,0},{1,0,0},{0.25,0.25,0}},
182:       B[2][3]   = {{1.0/12.0,1.0/12.0,5.0/6.0},{1.0/12.0,1.0/12.0,-1.0/6.0}},
183:       U[3][2]   = {{1,0},{1,10},{1,-1}},
184:       V[2][2]   = {{1,0},{0,1}},
185:       S[2]      = {1,0},
186:       F[2]      = {1,0},
187:       Fembed[2] = {1,1-GAMMA},
188:       Ferror[2] = {0,1},
189:       Serror[2] = {1,0};
190:     TSGLEERegister(TSGLEE23,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
191:   }
192:   {
193: #undef GAMMA
194: #define GAMMA 0.0
195:     /* y-y~ form */
196:     const PetscInt
197:       p = 2,
198:       s = 4,
199:       r = 2;
200:     const PetscReal
201:       A[4][4]   = {{0,0,0,0},{0.75,0,0,0},{0.25,29.0/60.0,0,0},{-21.0/44.0,145.0/44.0,-20.0/11.0,0}},
202:       B[2][4]   = {{109.0/275.0,58.0/75.0,-37.0/110.0,1.0/6.0},{3.0/11.0,0,75.0/88.0,-1.0/8.0}},
203:       U[4][2]   = {{0,1},{75.0/58.0,-17.0/58.0},{0,1},{0,1}},
204:       V[2][2]   = {{1,0},{0,1}},
205:       S[2]      = {1,1},
206:       F[2]      = {1,0},
207:       Fembed[2] = {0,1},
208:       Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
209:       Serror[2] = {1.0-GAMMA,1.0};
210:       TSGLEERegister(TSGLEE24,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
211:   }
212:   {
213: #undef GAMMA
214: #define GAMMA 0.0
215:     /* y-y~ form */
216:     const PetscInt
217:       p = 2,
218:       s = 5,
219:       r = 2;
220:     const PetscReal
221:       A[5][5]   = {{0,0,0,0,0},
222:                    {-0.94079244066783383269,0,0,0,0},
223:                    {0.64228187778301907108,0.10915356933958500042,0,0,0},
224:                    {-0.51764297742287450812,0.74414270351096040738,-0.71404164927824538121,0,0},
225:                    {-0.44696561556825969206,-0.76768425657590196518,0.20111608138142987881,0.93828186737840469796,0}},
226:       B[2][5]   = {{-0.029309178948150356153,-0.49671981884013874923,0.34275801517650053274,0.32941112623949194988,0.85385985637229662276},
227:                    {0.78133219686062535272,0.074238691892675897635,0.57957363498384957966,-0.24638502829674959968,-0.18875949544040123033}},
228:       U[5][2]   = {{0.16911424754448327735,0.83088575245551672265},
229:                    {0.53638465733199574340,0.46361534266800425660},
230:                    {0.39901579167169582526,0.60098420832830417474},
231:                    {0.87689005530618575480,0.12310994469381424520},
232:                    {0.99056100455550913009,0.0094389954444908699092}},
233:       V[2][2]   = {{1,0},{0,1}},
234:       S[2]      = {1,1},
235:       F[2]      = {1,0},
236:       Fembed[2] = {0,1},
237:       Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
238:       Serror[2] = {1.0-GAMMA,1.0};
239:     TSGLEERegister(TSGLEE25I,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
240:   }
241:   {
242: #undef GAMMA
243: #define GAMMA 0.0
244:     /* y-y~ form */
245:     const PetscInt
246:       p = 3,
247:       s = 5,
248:       r = 2;
249:     const PetscReal
250:       A[5][5]   = {{0,0,0,0,0},
251:                    {- 2169604947363702313.0 /  24313474998937147335.0,0,0,0,0},
252:                    {46526746497697123895.0 /  94116917485856474137.0,-10297879244026594958.0 /  49199457603717988219.0,0,0,0},
253:                    {23364788935845982499.0 /  87425311444725389446.0,-79205144337496116638.0 / 148994349441340815519.0,40051189859317443782.0 /  36487615018004984309.0,0,0},
254:                    {42089522664062539205.0 / 124911313006412840286.0,-15074384760342762939.0 / 137927286865289746282.0,-62274678522253371016.0 / 125918573676298591413.0,13755475729852471739.0 /  79257927066651693390.0,0}},
255:       B[2][5]   = {{61546696837458703723.0 /  56982519523786160813.0,-55810892792806293355.0 / 206957624151308356511.0,24061048952676379087.0 / 158739347956038723465.0,3577972206874351339.0 /   7599733370677197135.0,-59449832954780563947.0 / 137360038685338563670.0},
256:                    {- 9738262186984159168.0 /  99299082461487742983.0,-32797097931948613195.0 /  61521565616362163366.0,42895514606418420631.0 /  71714201188501437336.0,22608567633166065068.0 /  55371917805607957003.0,94655809487476459565.0 / 151517167160302729021.0}},
257:       U[5][2]   = {{70820309139834661559.0 /  80863923579509469826.0,10043614439674808267.0 /  80863923579509469826.0},
258:                    {161694774978034105510.0 / 106187653640211060371.0,-55507121337823045139.0 / 106187653640211060371.0},
259:                    {78486094644566264568.0 /  88171030896733822981.0,9684936252167558413.0 /  88171030896733822981.0},
260:                    {65394922146334854435.0 /  84570853840405479554.0,19175931694070625119.0 /  84570853840405479554.0},
261:                    {8607282770183754108.0 / 108658046436496925911.0,100050763666313171803.0 / 108658046436496925911.0}},
262:       V[2][2]   = {{1,0},{0,1}},
263:       S[2]      = {1,1},
264:       F[2]      = {1,0},
265:       Fembed[2] = {0,1},
266:       Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
267:       Serror[2] = {1.0-GAMMA,1.0};
268:     TSGLEERegister(TSGLEE35,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
269:   }
270:   {
271: #undef GAMMA
272: #define GAMMA 0.25
273:     /* y-eps form */
274:     const PetscInt
275:       p = 2,
276:       s = 6,
277:       r = 2;
278:     const PetscReal
279:       A[6][6]   = {{0,0,0,0,0,0},
280:                    {1,0,0,0,0,0},
281:                    {0,0,0,0,0,0},
282:                    {0,0,0.5,0,0,0},
283:                    {0,0,0.25,0.25,0,0},
284:                    {0,0,0.25,0.25,0.5,0}},
285:       B[2][6]   = {{0.5,0.5,0,0,0,0},
286:                    {-2.0/3.0,-2.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0}},
287:       U[6][2]   = {{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75}},
288:       V[2][2]   = {{1,0},{0,1}},
289:       S[2]      = {1,0},
290:       F[2]      = {1,0},
291:       Fembed[2] = {1,1-GAMMA},
292:       Ferror[2] = {0,1},
293:       Serror[2] = {1,0};
294:     TSGLEERegister(TSGLEEEXRK2A,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
295:   }
296:   {
297: #undef GAMMA
298: #define GAMMA 0.0
299:     /* y-eps form */
300:     const PetscInt
301:       p = 3,
302:       s = 8,
303:       r = 2;
304:     const PetscReal
305:       A[8][8]   = {{0,0,0,0,0,0,0,0},
306:                    {0.5,0,0,0,0,0,0,0},
307:                    {-1,2,0,0,0,0,0,0},
308:                    {1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
309:                    {0,0,0,0,0,0,0,0},
310:                    {-7.0/24.0,1.0/3.0,1.0/12.0,-1.0/8.0,0.5,0,0,0},
311:                    {7.0/6.0,-4.0/3.0,-1.0/3.0,0.5,-1.0,2.0,0,0},
312:                    {0,0,0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}},
313:       B[2][8]   = {{1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
314:                    {-1.0/6.0,-2.0/3.0,-1.0/6.0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}},
315:       U[8][2]   = {{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,1}},
316:       V[2][2]   = {{1,0},{0,1}},
317:       S[2]      = {1,0},
318:       F[2]      = {1,0},
319:       Fembed[2] = {1,1-GAMMA},
320:       Ferror[2] = {0,1},
321:       Serror[2] = {1,0};
322:     TSGLEERegister(TSGLEERK32G1,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
323:   }
324:   {
325: #undef GAMMA
326: #define GAMMA 0.25
327:     /* y-eps form */
328:     const PetscInt
329:       p = 2,
330:       s = 9,
331:       r = 2;
332:     const PetscReal
333:       A[9][9]   = {{0,0,0,0,0,0,0,0,0},
334:                    {0.585786437626904966,0,0,0,0,0,0,0,0},
335:                    {0.149999999999999994,0.849999999999999978,0,0,0,0,0,0,0},
336:                    {0,0,0,0,0,0,0,0,0},
337:                    {0,0,0,0.292893218813452483,0,0,0,0,0},
338:                    {0,0,0,0.0749999999999999972,0.424999999999999989,0,0,0,0},
339:                    {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0,0,0},
340:                    {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.292893218813452483,0,0},
341:                    {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.0749999999999999972,0.424999999999999989,0}},
342:       B[2][9]   = {{0.353553390593273786,0.353553390593273786,0.292893218813452483,0,0,0,0,0,0},
343:                    {-0.471404520791031678,-0.471404520791031678,-0.390524291751269959,0.235702260395515839,0.235702260395515839,0.195262145875634979,0.235702260395515839,0.235702260395515839,0.195262145875634979}},
344:       U[9][2]   = {{1,0},{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75},{1,0.75},{1,0.75}},
345:       V[2][2]   = {{1,0},{0,1}},
346:       S[2]      = {1,0},
347:       F[2]      = {1,0},
348:       Fembed[2] = {1,1-GAMMA},
349:       Ferror[2] = {0,1},
350:       Serror[2] = {1,0};
351:     TSGLEERegister(TSGLEERK285EX,p,s,r,GAMMA,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);
352:   }

354:   return(0);
355: }

357: /*@C
358:    TSGLEERegisterDestroy - Frees the list of schemes that were registered by TSGLEERegister().

360:    Not Collective

362:    Level: advanced

364: .seealso: TSGLEERegister(), TSGLEERegisterAll()
365: @*/
366: PetscErrorCode TSGLEERegisterDestroy(void)
367: {
369:   GLEETableauLink link;

372:   while ((link = GLEETableauList)) {
373:     GLEETableau t = &link->tab;
374:     GLEETableauList = link->next;
375:     PetscFree5(t->A,t->B,t->U,t->V,t->c);
376:     PetscFree2(t->S,t->F);               
377:     PetscFree (t->Fembed);               
378:     PetscFree (t->Ferror);               
379:     PetscFree (t->Serror);               
380:     PetscFree (t->binterp);              
381:     PetscFree (t->name);                 
382:     PetscFree (link);                    
383:   }
384:   TSGLEERegisterAllCalled = PETSC_FALSE;
385:   return(0);
386: }

388: /*@C
389:   TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called
390:   from TSInitializePackage().

392:   Level: developer

394: .seealso: PetscInitialize()
395: @*/
396: PetscErrorCode TSGLEEInitializePackage(void)
397: {

401:   if (TSGLEEPackageInitialized) return(0);
402:   TSGLEEPackageInitialized = PETSC_TRUE;
403:   TSGLEERegisterAll();
404:   PetscObjectComposedDataRegister(&explicit_stage_time_id);
405:   PetscRegisterFinalize(TSGLEEFinalizePackage);
406:   return(0);
407: }

409: /*@C
410:   TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is
411:   called from PetscFinalize().

413:   Level: developer

415: .seealso: PetscFinalize()
416: @*/
417: PetscErrorCode TSGLEEFinalizePackage(void)
418: {

422:   TSGLEEPackageInitialized = PETSC_FALSE;
423:   TSGLEERegisterDestroy();
424:   return(0);
425: }

427: /*@C
428:    TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau

430:    Not Collective, but the same schemes should be registered on all processes on which they will be used

432:    Input Parameters:
433: +  name   - identifier for method
434: .  order  - order of method
435: .  s - number of stages
436: .  r - number of steps
437: .  gamma - LTE ratio
438: .  A - stage coefficients (dimension s*s, row-major)
439: .  B - step completion coefficients (dimension r*s, row-major)
440: .  U - method coefficients (dimension s*r, row-major)
441: .  V - method coefficients (dimension r*r, row-major)
442: .  S - starting coefficients
443: .  F - finishing coefficients
444: .  c - abscissa (dimension s; NULL to use row sums of A)
445: .  Fembed - step completion coefficients for embedded method
446: .  Ferror - error computation coefficients
447: .  Serror - error initialization coefficients
448: .  pinterp - order of interpolation (0 if unavailable)
449: -  binterp - array of interpolation coefficients (NULL if unavailable)

451:    Notes:
452:    Several GLEE methods are provided, this function is only needed to create new methods.

454:    Level: advanced

456: .seealso: TSGLEE
457: @*/
458: PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r,
459:                               PetscReal gamma,
460:                               const PetscReal A[],const PetscReal B[],
461:                               const PetscReal U[],const PetscReal V[],
462:                               const PetscReal S[],const PetscReal F[],
463:                               const PetscReal c[],
464:                               const PetscReal Fembed[],const PetscReal Ferror[],
465:                               const PetscReal Serror[],
466:                               PetscInt pinterp, const PetscReal binterp[])
467: {
468:   PetscErrorCode    ierr;
469:   GLEETableauLink   link;
470:   GLEETableau       t;
471:   PetscInt          i,j;

474:   TSGLEEInitializePackage();
475:   PetscNew(&link);
476:   t        = &link->tab;
477:   PetscStrallocpy(name,&t->name);
478:   t->order = order;
479:   t->s     = s;
480:   t->r     = r;
481:   t->gamma = gamma;
482:   PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U);
483:   PetscMalloc2(r,&t->S,r,&t->F);
484:   PetscArraycpy(t->A,A,s*s);
485:   PetscArraycpy(t->B,B,r*s);
486:   PetscArraycpy(t->U,U,s*r);
487:   PetscArraycpy(t->V,V,r*r);
488:   PetscArraycpy(t->S,S,r);
489:   PetscArraycpy(t->F,F,r);
490:   if (c) {
491:     PetscArraycpy(t->c,c,s);
492:   } else {
493:     for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
494:   }
495:   PetscMalloc1(r,&t->Fembed);
496:   PetscMalloc1(r,&t->Ferror);
497:   PetscMalloc1(r,&t->Serror);
498:   PetscArraycpy(t->Fembed,Fembed,r);
499:   PetscArraycpy(t->Ferror,Ferror,r);
500:   PetscArraycpy(t->Serror,Serror,r);
501:   t->pinterp = pinterp;
502:   PetscMalloc1(s*pinterp,&t->binterp);
503:   PetscArraycpy(t->binterp,binterp,s*pinterp);

505:   link->next      = GLEETableauList;
506:   GLEETableauList = link;
507:   return(0);
508: }

510: static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done)
511: {
512:   TS_GLEE         *glee = (TS_GLEE*) ts->data;
513:   GLEETableau     tab = glee->tableau;
514:   PetscReal       h, *B = tab->B, *V = tab->V,
515:                   *F = tab->F,
516:                   *Fembed = tab->Fembed;
517:   PetscInt        s = tab->s, r = tab->r, i, j;
518:   Vec             *Y = glee->Y, *YdotStage = glee->YdotStage;
519:   PetscScalar     *ws = glee->swork, *wr = glee->rwork;
520:   PetscErrorCode  ierr;


524:   switch (glee->status) {
525:     case TS_STEP_INCOMPLETE:
526:     case TS_STEP_PENDING:
527:       h = ts->time_step; break;
528:     case TS_STEP_COMPLETE:
529:       h = ts->ptime - ts->ptime_prev; break;
530:     default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
531:   }

533:   if (order == tab->order) {

535:     /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
536:              or TS_STEP_COMPLETE, glee->X has the solution at the
537:              beginning of the time step. So no need to roll-back.
538:     */
539:     if (glee->status == TS_STEP_INCOMPLETE) {
540:       for (i=0; i<r; i++) {
541:         VecZeroEntries(Y[i]);
542:         for (j=0; j<r; j++) wr[j] = V[i*r+j];
543:         VecMAXPY(Y[i],r,wr,glee->X);
544:         for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
545:         VecMAXPY(Y[i],s,ws,YdotStage);
546:       }
547:       VecZeroEntries(X);
548:       for (j=0; j<r; j++) wr[j] = F[j];
549:       VecMAXPY(X,r,wr,Y);
550:     } else {VecCopy(ts->vec_sol,X);}
551:     return(0);

553:   } else if (order == tab->order-1) {

555:     /* Complete with the embedded method (Fembed) */
556:     for (i=0; i<r; i++) {
557:       VecZeroEntries(Y[i]);
558:       for (j=0; j<r; j++) wr[j] = V[i*r+j];
559:       VecMAXPY(Y[i],r,wr,glee->X);
560:       for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
561:       VecMAXPY(Y[i],s,ws,YdotStage);
562:     }
563:     VecZeroEntries(X);
564:     for (j=0; j<r; j++) wr[j] = Fembed[j];
565:     VecMAXPY(X,r,wr,Y);

567:     if (done) *done = PETSC_TRUE;
568:     return(0);
569:   }
570:   if (done) *done = PETSC_FALSE;
571:   else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
572:   return(0);
573: }

575: static PetscErrorCode TSStep_GLEE(TS ts)
576: {
577:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
578:   GLEETableau     tab = glee->tableau;
579:   const PetscInt  s = tab->s, r = tab->r;
580:   PetscReal       *A = tab->A, *U = tab->U,
581:                   *F = tab->F,
582:                   *c = tab->c;
583:   Vec             *Y = glee->Y, *X = glee->X,
584:                   *YStage = glee->YStage,
585:                   *YdotStage = glee->YdotStage,
586:                   W = glee->W;
587:   SNES            snes;
588:   PetscScalar     *ws = glee->swork, *wr = glee->rwork;
589:   TSAdapt         adapt;
590:   PetscInt        i,j,reject,next_scheme,its,lits;
591:   PetscReal       next_time_step;
592:   PetscReal       t;
593:   PetscBool       accept;
594:   PetscErrorCode  ierr;

597:   PetscCitationsRegister(citation,&cited);

599:   for (i=0; i<r; i++) { VecCopy(Y[i],X[i]); }

601:   TSGetSNES(ts,&snes);
602:   next_time_step = ts->time_step;
603:   t              = ts->ptime;
604:   accept         = PETSC_TRUE;
605:   glee->status   = TS_STEP_INCOMPLETE;

607:   for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {

609:     PetscReal h = ts->time_step;
610:     TSPreStep(ts);

612:     for (i=0; i<s; i++) {

614:       glee->stage_time = t + h*c[i];
615:       TSPreStage(ts,glee->stage_time);

617:       if (A[i*s+i] == 0) {  /* Explicit stage */
618:         VecZeroEntries(YStage[i]);
619:         for (j=0; j<r; j++) wr[j] = U[i*r+j];
620:         VecMAXPY(YStage[i],r,wr,X);
621:         for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
622:         VecMAXPY(YStage[i],i,ws,YdotStage);
623:       } else {              /* Implicit stage */
624:         glee->scoeff = 1.0/A[i*s+i];
625:         /* compute right-hand-side */
626:         VecZeroEntries(W);
627:         for (j=0; j<r; j++) wr[j] = U[i*r+j];
628:         VecMAXPY(W,r,wr,X);
629:         for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
630:         VecMAXPY(W,i,ws,YdotStage);
631:         VecScale(W,glee->scoeff/h);
632:         /* set initial guess */
633:         VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i]);
634:         /* solve for this stage */
635:         SNESSolve(snes,W,YStage[i]);
636:         SNESGetIterationNumber(snes,&its);
637:         SNESGetLinearSolveIterations(snes,&lits);
638:         ts->snes_its += its; ts->ksp_its += lits;
639:       }
640:       TSGetAdapt(ts,&adapt);
641:       TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept);
642:       if (!accept) goto reject_step;
643:       TSPostStage(ts,glee->stage_time,i,YStage);
644:       TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i]);
645:     }
646:     TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
647:     glee->status = TS_STEP_PENDING;

649:     /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
650:     TSGetAdapt(ts,&adapt);
651:     TSAdaptCandidatesClear(adapt);
652:     TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);
653:     TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
654:     if (accept) {
655:       /* ignore next_scheme for now */
656:       ts->ptime     += ts->time_step;
657:       ts->time_step  = next_time_step;
658:       glee->status = TS_STEP_COMPLETE;
659:       /* compute and store the global error */
660:       /* Note: this is not needed if TSAdaptGLEE is not used */
661:       TSGetTimeError(ts,0,&(glee->yGErr));
662:       PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
663:       break;
664:     } else {                    /* Roll back the current step */
665:       for (j=0; j<r; j++) wr[j] = F[j];
666:       VecMAXPY(ts->vec_sol,r,wr,X);
667:       ts->time_step = next_time_step;
668:       glee->status  = TS_STEP_INCOMPLETE;
669:     }
670: reject_step: continue;
671:   }
672:   if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
673:   return(0);
674: }

676: static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X)
677: {
678:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
679:   PetscInt        s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j;
680:   PetscReal       h,tt,t;
681:   PetscScalar     *b;
682:   const PetscReal *B = glee->tableau->binterp;
683:   PetscErrorCode  ierr;

686:   if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name);
687:   switch (glee->status) {
688:     case TS_STEP_INCOMPLETE:
689:     case TS_STEP_PENDING:
690:       h = ts->time_step;
691:       t = (itime - ts->ptime)/h;
692:       break;
693:     case TS_STEP_COMPLETE:
694:       h = ts->ptime - ts->ptime_prev;
695:       t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
696:       break;
697:     default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
698:   }
699:   PetscMalloc1(s,&b);
700:   for (i=0; i<s; i++) b[i] = 0;
701:   for (j=0,tt=t; j<pinterp; j++,tt*=t) {
702:     for (i=0; i<s; i++) {
703:       b[i]  += h * B[i*pinterp+j] * tt;
704:     }
705:   }
706:   VecCopy(glee->YStage[0],X);
707:   VecMAXPY(X,s,b,glee->YdotStage);
708:   PetscFree(b);
709:   return(0);
710: }

712: /*------------------------------------------------------------*/
713: static PetscErrorCode TSReset_GLEE(TS ts)
714: {
715:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
716:   PetscInt       s, r;

720:   if (!glee->tableau) return(0);
721:   s    = glee->tableau->s;
722:   r    = glee->tableau->r;
723:   VecDestroyVecs(r,&glee->Y);
724:   VecDestroyVecs(r,&glee->X);
725:   VecDestroyVecs(s,&glee->YStage);
726:   VecDestroyVecs(s,&glee->YdotStage);
727:   VecDestroy(&glee->Ydot);
728:   VecDestroy(&glee->yGErr);
729:   VecDestroy(&glee->W);
730:   PetscFree2(glee->swork,glee->rwork);
731:   return(0);
732: }

734: static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot)
735: {
736:   TS_GLEE     *glee = (TS_GLEE*)ts->data;

740:   if (Ydot) {
741:     if (dm && dm != ts->dm) {
742:       DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
743:     } else *Ydot = glee->Ydot;
744:   }
745:   return(0);
746: }

748: static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot)
749: {

753:   if (Ydot) {
754:     if (dm && dm != ts->dm) {
755:       DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
756:     }
757:   }
758:   return(0);
759: }

761: /*
762:   This defines the nonlinear equation that is to be solved with SNES
763: */
764: static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts)
765: {
766:   TS_GLEE       *glee = (TS_GLEE*)ts->data;
767:   DM             dm,dmsave;
768:   Vec            Ydot;
769:   PetscReal      shift = glee->scoeff / ts->time_step;

773:   SNESGetDM(snes,&dm);
774:   TSGLEEGetVecs(ts,dm,&Ydot);
775:   /* Set Ydot = shift*X */
776:   VecCopy(X,Ydot);
777:   VecScale(Ydot,shift);
778:   dmsave = ts->dm;
779:   ts->dm = dm;

781:   TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);

783:   ts->dm = dmsave;
784:   TSGLEERestoreVecs(ts,dm,&Ydot);
785:   return(0);
786: }

788: static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts)
789: {
790:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
791:   DM             dm,dmsave;
792:   Vec            Ydot;
793:   PetscReal      shift = glee->scoeff / ts->time_step;

797:   SNESGetDM(snes,&dm);
798:   TSGLEEGetVecs(ts,dm,&Ydot);
799:   /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
800:   dmsave = ts->dm;
801:   ts->dm = dm;

803:   TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);

805:   ts->dm = dmsave;
806:   TSGLEERestoreVecs(ts,dm,&Ydot);
807:   return(0);
808: }

810: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx)
811: {
813:   return(0);
814: }

816: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
817: {
819:   return(0);
820: }

822: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx)
823: {
825:   return(0);
826: }

828: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
829: {
831:   return(0);
832: }

834: static PetscErrorCode TSSetUp_GLEE(TS ts)
835: {
836:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
837:   GLEETableau    tab;
838:   PetscInt       s,r;
840:   DM             dm;

843:   if (!glee->tableau) {
844:     TSGLEESetType(ts,TSGLEEDefaultType);
845:   }
846:   tab  = glee->tableau;
847:   s    = tab->s;
848:   r    = tab->r;
849:   VecDuplicateVecs(ts->vec_sol,r,&glee->Y);
850:   VecDuplicateVecs(ts->vec_sol,r,&glee->X);
851:   VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);
852:   VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);
853:   VecDuplicate(ts->vec_sol,&glee->Ydot);
854:   VecDuplicate(ts->vec_sol,&glee->yGErr);
855:   VecZeroEntries(glee->yGErr);
856:   VecDuplicate(ts->vec_sol,&glee->W);
857:   PetscMalloc2(s,&glee->swork,r,&glee->rwork);
858:   TSGetDM(ts,&dm);
859:   DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
860:   DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
861:   return(0);
862: }

864: PetscErrorCode TSStartingMethod_GLEE(TS ts)
865: {
866:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
867:   GLEETableau    tab  = glee->tableau;
868:   PetscInt       r=tab->r,i;
869:   PetscReal      *S=tab->S;

873:   for (i=0; i<r; i++) {
874:     VecZeroEntries(glee->Y[i]);
875:     VecAXPY(glee->Y[i],S[i],ts->vec_sol);
876:   }

878:   return(0);
879: }

881: /*------------------------------------------------------------*/

883: static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts)
884: {
886:   char           gleetype[256];

889:   PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");
890:   {
891:     GLEETableauLink link;
892:     PetscInt        count,choice;
893:     PetscBool       flg;
894:     const char      **namelist;

896:     PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));
897:     for (link=GLEETableauList,count=0; link; link=link->next,count++) ;
898:     PetscMalloc1(count,(char***)&namelist);
899:     for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
900:     PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);
901:     TSGLEESetType(ts,flg ? namelist[choice] : gleetype);
902:     PetscFree(namelist);
903:   }
904:   PetscOptionsTail();
905:   return(0);
906: }

908: static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer)
909: {
910:   TS_GLEE        *glee   = (TS_GLEE*)ts->data;
911:   GLEETableau    tab  = glee->tableau;
912:   PetscBool      iascii;

916:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
917:   if (iascii) {
918:     TSGLEEType gleetype;
919:     char       buf[512];
920:     TSGLEEGetType(ts,&gleetype);
921:     PetscViewerASCIIPrintf(viewer,"  GLEE type %s\n",gleetype);
922:     PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
923:     PetscViewerASCIIPrintf(viewer,"  Abscissa     c = %s\n",buf);
924:     /* Note: print out r as well */
925:   }
926:   return(0);
927: }

929: static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer)
930: {
932:   SNES           snes;
933:   TSAdapt        tsadapt;

936:   TSGetAdapt(ts,&tsadapt);
937:   TSAdaptLoad(tsadapt,viewer);
938:   TSGetSNES(ts,&snes);
939:   SNESLoad(snes,viewer);
940:   /* function and Jacobian context for SNES when used with TS is always ts object */
941:   SNESSetFunction(snes,NULL,NULL,ts);
942:   SNESSetJacobian(snes,NULL,NULL,NULL,ts);
943:   return(0);
944: }

946: /*@C
947:   TSGLEESetType - Set the type of GLEE scheme

949:   Logically collective

951:   Input Parameter:
952: +  ts - timestepping context
953: -  gleetype - type of GLEE-scheme

955:   Level: intermediate

957: .seealso: TSGLEEGetType(), TSGLEE
958: @*/
959: PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype)
960: {

966:   PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));
967:   return(0);
968: }

970: /*@C
971:   TSGLEEGetType - Get the type of GLEE scheme

973:   Logically collective

975:   Input Parameter:
976: .  ts - timestepping context

978:   Output Parameter:
979: .  gleetype - type of GLEE-scheme

981:   Level: intermediate

983: .seealso: TSGLEESetType()
984: @*/
985: PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype)
986: {

991:   PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));
992:   return(0);
993: }

995: PetscErrorCode  TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype)
996: {
997:   TS_GLEE     *glee = (TS_GLEE*)ts->data;

1001:   if (!glee->tableau) {
1002:     TSGLEESetType(ts,TSGLEEDefaultType);
1003:   }
1004:   *gleetype = glee->tableau->name;
1005:   return(0);
1006: }
1007: PetscErrorCode  TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype)
1008: {
1009:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1010:   PetscErrorCode  ierr;
1011:   PetscBool       match;
1012:   GLEETableauLink link;

1015:   if (glee->tableau) {
1016:     PetscStrcmp(glee->tableau->name,gleetype,&match);
1017:     if (match) return(0);
1018:   }
1019:   for (link = GLEETableauList; link; link=link->next) {
1020:     PetscStrcmp(link->tab.name,gleetype,&match);
1021:     if (match) {
1022:       TSReset_GLEE(ts);
1023:       glee->tableau = &link->tab;
1024:       return(0);
1025:     }
1026:   }
1027:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype);
1028: }

1030: static PetscErrorCode  TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y)
1031: {
1032:   TS_GLEE *glee = (TS_GLEE*)ts->data;

1035:   if (ns) *ns = glee->tableau->s;
1036:   if (Y)  *Y  = glee->YStage;
1037:   return(0);
1038: }

1040: PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y)
1041: {
1042:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1043:   GLEETableau     tab   = glee->tableau;
1044:   PetscErrorCode  ierr;

1047:   if (!Y) *n = tab->r;
1048:   else {
1049:     if ((*n >= 0) && (*n < tab->r)) {
1050:       VecCopy(glee->Y[*n],*Y);
1051:     } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1);
1052:   }
1053:   return(0);
1054: }

1056: PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X)
1057: {
1058:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1059:   GLEETableau     tab   = glee->tableau;
1060:   PetscReal       *F    = tab->Fembed;
1061:   PetscInt        r     = tab->r;
1062:   Vec             *Y    = glee->Y;
1063:   PetscScalar     *wr   = glee->rwork;
1064:   PetscInt        i;
1065:   PetscErrorCode  ierr;

1068:   VecZeroEntries(*X);
1069:   for (i=0; i<r; i++) wr[i] = F[i];
1070:   VecMAXPY((*X),r,wr,Y);
1071:   return(0);
1072: }

1074: PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X)
1075: {
1076:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1077:   GLEETableau     tab   = glee->tableau;
1078:   PetscReal       *F    = tab->Ferror;
1079:   PetscInt        r     = tab->r;
1080:   Vec             *Y    = glee->Y;
1081:   PetscScalar     *wr   = glee->rwork;
1082:   PetscInt        i;
1083:   PetscErrorCode  ierr;

1086:   VecZeroEntries(*X);
1087:   if (n==0) {
1088:     for (i=0; i<r; i++) wr[i] = F[i];
1089:     VecMAXPY((*X),r,wr,Y);
1090:   } else if (n==-1) {
1091:     *X=glee->yGErr;
1092:   }
1093:   return(0);
1094: }

1096: PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X)
1097: {
1098:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1099:   GLEETableau     tab   = glee->tableau;
1100:   PetscReal       *S    = tab->Serror;
1101:   PetscInt        r     = tab->r,i;
1102:   Vec             *Y    = glee->Y;
1103:   PetscErrorCode  ierr;

1106:   if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r);
1107:   for (i=1; i<r; i++) {
1108:     VecCopy(ts->vec_sol,Y[i]);
1109:     VecAXPBY(Y[i],S[0],S[1],X);
1110:     VecCopy(X,glee->yGErr);
1111:   }
1112:   return(0);
1113: }

1115: static PetscErrorCode TSDestroy_GLEE(TS ts)
1116: {

1120:   TSReset_GLEE(ts);
1121:   if (ts->dm) {
1122:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
1123:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
1124:   }
1125:   PetscFree(ts->data);
1126:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);
1127:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);
1128:   return(0);
1129: }

1131: /* ------------------------------------------------------------ */
1132: /*MC
1133:       TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes

1135:   The user should provide the right hand side of the equation
1136:   using TSSetRHSFunction().

1138:   Notes:
1139:   The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type

1141:   Level: beginner

1143: .seealso:  TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(),
1144:            TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A,
1145:            TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister()

1147: M*/
1148: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
1149: {
1150:   TS_GLEE         *th;
1151:   PetscErrorCode  ierr;

1154:   TSGLEEInitializePackage();

1156:   ts->ops->reset                  = TSReset_GLEE;
1157:   ts->ops->destroy                = TSDestroy_GLEE;
1158:   ts->ops->view                   = TSView_GLEE;
1159:   ts->ops->load                   = TSLoad_GLEE;
1160:   ts->ops->setup                  = TSSetUp_GLEE;
1161:   ts->ops->step                   = TSStep_GLEE;
1162:   ts->ops->interpolate            = TSInterpolate_GLEE;
1163:   ts->ops->evaluatestep           = TSEvaluateStep_GLEE;
1164:   ts->ops->setfromoptions         = TSSetFromOptions_GLEE;
1165:   ts->ops->getstages              = TSGetStages_GLEE;
1166:   ts->ops->snesfunction           = SNESTSFormFunction_GLEE;
1167:   ts->ops->snesjacobian           = SNESTSFormJacobian_GLEE;
1168:   ts->ops->getsolutioncomponents  = TSGetSolutionComponents_GLEE;
1169:   ts->ops->getauxsolution         = TSGetAuxSolution_GLEE;
1170:   ts->ops->gettimeerror           = TSGetTimeError_GLEE;
1171:   ts->ops->settimeerror           = TSSetTimeError_GLEE;
1172:   ts->ops->startingmethod         = TSStartingMethod_GLEE;
1173:   ts->default_adapt_type          = TSADAPTGLEE;

1175:   ts->usessnes = PETSC_TRUE;

1177:   PetscNewLog(ts,&th);
1178:   ts->data = (void*)th;

1180:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);
1181:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);
1182:   return(0);
1183: }