Actual source code: glee.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1: /*
  2:   Code for time stepping with the General Linear with Error Estimation method


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

  8:   Udot = F(t,U)

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

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

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

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

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

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

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

 72:      Level: advanced

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

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

 82:      Level: advanced

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

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

 92:      Level: advanced

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

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

102:      Level: advanced

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

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

112:      Level: advanced

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

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

122:      Level: advanced

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

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

132:      Level: advanced

134: .seealso: TSGLEE
135: M*/

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

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

142:   Level: advanced

144: .keywords: TS, TSGLEE, register, all

146: .seealso:  TSGLEERegisterDestroy()
147: @*/
148: PetscErrorCode TSGLEERegisterAll(void)
149: {

153:   if (TSGLEERegisterAllCalled) return(0);
154:   TSGLEERegisterAllCalled = PETSC_TRUE;

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

357:   return(0);
358: }

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

363:    Not Collective

365:    Level: advanced

367: .keywords: TSGLEE, register, destroy
368: .seealso: TSGLEERegister(), TSGLEERegisterAll()
369: @*/
370: PetscErrorCode TSGLEERegisterDestroy(void)
371: {
373:   GLEETableauLink link;

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

392: /*@C
393:   TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called
394:   from TSInitializePackage().

396:   Level: developer

398: .keywords: TS, TSGLEE, initialize, package
399: .seealso: PetscInitialize()
400: @*/
401: PetscErrorCode TSGLEEInitializePackage(void)
402: {

406:   if (TSGLEEPackageInitialized) return(0);
407:   TSGLEEPackageInitialized = PETSC_TRUE;
408:   TSGLEERegisterAll();
409:   PetscObjectComposedDataRegister(&explicit_stage_time_id);
410:   PetscRegisterFinalize(TSGLEEFinalizePackage);
411:   return(0);
412: }

414: /*@C
415:   TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is
416:   called from PetscFinalize().

418:   Level: developer

420: .keywords: Petsc, destroy, package
421: .seealso: PetscFinalize()
422: @*/
423: PetscErrorCode TSGLEEFinalizePackage(void)
424: {

428:   TSGLEEPackageInitialized = PETSC_FALSE;
429:   TSGLEERegisterDestroy();
430:   return(0);
431: }

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

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

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

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

460:    Level: advanced

462: .keywords: TS, register

464: .seealso: TSGLEE
465: @*/
466: PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r,
467:                               PetscReal gamma,
468:                               const PetscReal A[],const PetscReal B[],
469:                               const PetscReal U[],const PetscReal V[],
470:                               const PetscReal S[],const PetscReal F[],
471:                               const PetscReal c[],
472:                               const PetscReal Fembed[],const PetscReal Ferror[],
473:                               const PetscReal Serror[],
474:                               PetscInt pinterp, const PetscReal binterp[])
475: {
476:   PetscErrorCode    ierr;
477:   GLEETableauLink   link;
478:   GLEETableau       t;
479:   PetscInt          i,j;

482:   TSGLEEInitializePackage();
483:   PetscMalloc(sizeof(*link),&link);
484:   PetscMemzero(link,sizeof(*link));
485:   t        = &link->tab;
486:   PetscStrallocpy(name,&t->name);
487:   t->order = order;
488:   t->s     = s;
489:   t->r     = r;
490:   t->gamma = gamma;
491:   PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U);
492:   PetscMalloc2(r,&t->S,r,&t->F);
493:   PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
494:   PetscMemcpy(t->B,B,r*s*sizeof(B[0]));
495:   PetscMemcpy(t->U,U,s*r*sizeof(U[0]));
496:   PetscMemcpy(t->V,V,r*r*sizeof(V[0]));
497:   PetscMemcpy(t->S,S,r  *sizeof(S[0]));
498:   PetscMemcpy(t->F,F,r  *sizeof(F[0]));
499:   if (c) {
500:     PetscMemcpy(t->c,c,s*sizeof(c[0]));
501:   } else {
502:     for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
503:   }
504:   PetscMalloc1(r,&t->Fembed);
505:   PetscMalloc1(r,&t->Ferror);
506:   PetscMalloc1(r,&t->Serror);
507:   PetscMemcpy(t->Fembed,Fembed,r*sizeof(Fembed[0]));
508:   PetscMemcpy(t->Ferror,Ferror,r*sizeof(Ferror[0]));
509:   PetscMemcpy(t->Serror,Serror,r*sizeof(Serror[0]));
510:   t->pinterp = pinterp;
511:   PetscMalloc1(s*pinterp,&t->binterp);
512:   PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));

514:   link->next      = GLEETableauList;
515:   GLEETableauList = link;
516:   return(0);
517: }

519: static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done)
520: {
521:   TS_GLEE         *glee = (TS_GLEE*) ts->data;
522:   GLEETableau     tab = glee->tableau;
523:   PetscReal       h, *B = tab->B, *V = tab->V,
524:                   *F = tab->F,
525:                   *Fembed = tab->Fembed;
526:   PetscInt        s = tab->s, r = tab->r, i, j;
527:   Vec             *Y = glee->Y, *YdotStage = glee->YdotStage;
528:   PetscScalar     *ws = glee->swork, *wr = glee->rwork;
529:   PetscErrorCode  ierr;


533:   switch (glee->status) {
534:     case TS_STEP_INCOMPLETE:
535:     case TS_STEP_PENDING:
536:       h = ts->time_step; break;
537:     case TS_STEP_COMPLETE:
538:       h = ts->ptime - ts->ptime_prev; break;
539:     default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
540:   }

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

544:     /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
545:              or TS_STEP_COMPLETE, glee->X has the solution at the
546:              beginning of the time step. So no need to roll-back.
547:     */
548:     if (glee->status == TS_STEP_INCOMPLETE) {
549:       for (i=0; i<r; i++) {
550:         VecZeroEntries(Y[i]);
551:         for (j=0; j<r; j++) wr[j] = V[i*r+j];
552:         VecMAXPY(Y[i],r,wr,glee->X);
553:         for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
554:         VecMAXPY(Y[i],s,ws,YdotStage);
555:       }
556:       VecZeroEntries(X);
557:       for (j=0; j<r; j++) wr[j] = F[j];
558:       VecMAXPY(X,r,wr,Y);
559:     } else {VecCopy(ts->vec_sol,X);}
560:     return(0);

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

564:     /* Complete with the embedded method (Fembed) */
565:     for (i=0; i<r; i++) {
566:       VecZeroEntries(Y[i]);
567:       for (j=0; j<r; j++) wr[j] = V[i*r+j];
568:       VecMAXPY(Y[i],r,wr,glee->X);
569:       for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
570:       VecMAXPY(Y[i],s,ws,YdotStage);
571:     }
572:     VecZeroEntries(X);
573:     for (j=0; j<r; j++) wr[j] = Fembed[j];
574:     VecMAXPY(X,r,wr,Y);

576:     if (done) *done = PETSC_TRUE;
577:     return(0);
578:   }
579:   if (done) *done = PETSC_FALSE;
580:   else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
581:   return(0);
582: }

584: static PetscErrorCode TSStep_GLEE(TS ts)
585: {
586:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
587:   GLEETableau     tab = glee->tableau;
588:   const PetscInt  s = tab->s, r = tab->r;
589:   PetscReal       *A = tab->A, *U = tab->U,
590:                   *F = tab->F,
591:                   *c = tab->c;
592:   Vec             *Y = glee->Y, *X = glee->X,
593:                   *YStage = glee->YStage,
594:                   *YdotStage = glee->YdotStage,
595:                   W = glee->W;
596:   SNES            snes;
597:   PetscScalar     *ws = glee->swork, *wr = glee->rwork;
598:   TSAdapt         adapt;
599:   PetscInt        i,j,reject,next_scheme,its,lits;
600:   PetscReal       next_time_step;
601:   PetscReal       t;
602:   PetscBool       accept;
603:   PetscErrorCode  ierr;

606:   PetscCitationsRegister(citation,&cited);

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

610:   TSGetSNES(ts,&snes);
611:   next_time_step = ts->time_step;
612:   t              = ts->ptime;
613:   accept         = PETSC_TRUE;
614:   glee->status   = TS_STEP_INCOMPLETE;

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

618:     PetscReal h = ts->time_step;
619:     TSPreStep(ts);

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

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

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

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

685: static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X)
686: {
687:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
688:   PetscInt        s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j;
689:   PetscReal       h,tt,t;
690:   PetscScalar     *b;
691:   const PetscReal *B = glee->tableau->binterp;
692:   PetscErrorCode  ierr;

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

721: /*------------------------------------------------------------*/
722: static PetscErrorCode TSReset_GLEE(TS ts)
723: {
724:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
725:   PetscInt       s, r;

729:   if (!glee->tableau) return(0);
730:   s    = glee->tableau->s;
731:   r    = glee->tableau->r;
732:   VecDestroyVecs(r,&glee->Y);
733:   VecDestroyVecs(r,&glee->X);
734:   VecDestroyVecs(s,&glee->YStage);
735:   VecDestroyVecs(s,&glee->YdotStage);
736:   VecDestroy(&glee->Ydot);
737:   VecDestroy(&glee->yGErr);
738:   VecDestroy(&glee->W);
739:   PetscFree2(glee->swork,glee->rwork);
740:   return(0);
741: }

743: static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot)
744: {
745:   TS_GLEE     *glee = (TS_GLEE*)ts->data;

749:   if (Ydot) {
750:     if (dm && dm != ts->dm) {
751:       DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
752:     } else *Ydot = glee->Ydot;
753:   }
754:   return(0);
755: }


758: static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot)
759: {

763:   if (Ydot) {
764:     if (dm && dm != ts->dm) {
765:       DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
766:     }
767:   }
768:   return(0);
769: }

771: /*
772:   This defines the nonlinear equation that is to be solved with SNES
773: */
774: static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts)
775: {
776:   TS_GLEE       *glee = (TS_GLEE*)ts->data;
777:   DM             dm,dmsave;
778:   Vec            Ydot;
779:   PetscReal      shift = glee->scoeff / ts->time_step;

783:   SNESGetDM(snes,&dm);
784:   TSGLEEGetVecs(ts,dm,&Ydot);
785:   /* Set Ydot = shift*X */
786:   VecCopy(X,Ydot);
787:   VecScale(Ydot,shift);
788:   dmsave = ts->dm;
789:   ts->dm = dm;

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

793:   ts->dm = dmsave;
794:   TSGLEERestoreVecs(ts,dm,&Ydot);
795:   return(0);
796: }

798: static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts)
799: {
800:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
801:   DM             dm,dmsave;
802:   Vec            Ydot;
803:   PetscReal      shift = glee->scoeff / ts->time_step;

807:   SNESGetDM(snes,&dm);
808:   TSGLEEGetVecs(ts,dm,&Ydot);
809:   /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
810:   dmsave = ts->dm;
811:   ts->dm = dm;

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

815:   ts->dm = dmsave;
816:   TSGLEERestoreVecs(ts,dm,&Ydot);
817:   return(0);
818: }

820: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx)
821: {
823:   return(0);
824: }

826: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
827: {
829:   return(0);
830: }


833: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx)
834: {
836:   return(0);
837: }

839: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
840: {
842:   return(0);
843: }

845: static PetscErrorCode TSSetUp_GLEE(TS ts)
846: {
847:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
848:   GLEETableau    tab;
849:   PetscInt       s,r;
851:   DM             dm;

854:   if (!glee->tableau) {
855:     TSGLEESetType(ts,TSGLEEDefaultType);
856:   }
857:   tab  = glee->tableau;
858:   s    = tab->s;
859:   r    = tab->r;
860:   VecDuplicateVecs(ts->vec_sol,r,&glee->Y);
861:   VecDuplicateVecs(ts->vec_sol,r,&glee->X);
862:   VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);
863:   VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);
864:   VecDuplicate(ts->vec_sol,&glee->Ydot);
865:   VecDuplicate(ts->vec_sol,&glee->yGErr);
866:   VecZeroEntries(glee->yGErr);
867:   VecDuplicate(ts->vec_sol,&glee->W);
868:   PetscMalloc2(s,&glee->swork,r,&glee->rwork);
869:   TSGetDM(ts,&dm);
870:   DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
871:   DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
872:   return(0);
873: }

875: PetscErrorCode TSStartingMethod_GLEE(TS ts)
876: {
877:   TS_GLEE        *glee = (TS_GLEE*)ts->data;
878:   GLEETableau    tab  = glee->tableau;
879:   PetscInt       r=tab->r,i;
880:   PetscReal      *S=tab->S;

884:   for (i=0; i<r; i++) {
885:     VecZeroEntries(glee->Y[i]);
886:     VecAXPY(glee->Y[i],S[i],ts->vec_sol);
887:   }

889:   return(0);
890: }


893: /*------------------------------------------------------------*/

895: static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts)
896: {
898:   char           gleetype[256];

901:   PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");
902:   {
903:     GLEETableauLink link;
904:     PetscInt        count,choice;
905:     PetscBool       flg;
906:     const char      **namelist;

908:     PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));
909:     for (link=GLEETableauList,count=0; link; link=link->next,count++) ;
910:     PetscMalloc1(count,(char***)&namelist);
911:     for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
912:     PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);
913:     TSGLEESetType(ts,flg ? namelist[choice] : gleetype);
914:     PetscFree(namelist);
915:   }
916:   PetscOptionsTail();
917:   return(0);
918: }

920: static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer)
921: {
922:   TS_GLEE        *glee   = (TS_GLEE*)ts->data;
923:   GLEETableau    tab  = glee->tableau;
924:   PetscBool      iascii;

928:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
929:   if (iascii) {
930:     TSGLEEType gleetype;
931:     char   buf[512];
932:     TSGLEEGetType(ts,&gleetype);
933:     PetscViewerASCIIPrintf(viewer,"  GLEE type %s\n",gleetype);
934:     PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
935:     PetscViewerASCIIPrintf(viewer,"  Abscissa     c = %s\n",buf);
936:     /* Note: print out r as well */
937:   }
938:   return(0);
939: }

941: static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer)
942: {
944:   SNES           snes;
945:   TSAdapt        tsadapt;

948:   TSGetAdapt(ts,&tsadapt);
949:   TSAdaptLoad(tsadapt,viewer);
950:   TSGetSNES(ts,&snes);
951:   SNESLoad(snes,viewer);
952:   /* function and Jacobian context for SNES when used with TS is always ts object */
953:   SNESSetFunction(snes,NULL,NULL,ts);
954:   SNESSetJacobian(snes,NULL,NULL,NULL,ts);
955:   return(0);
956: }

958: /*@C
959:   TSGLEESetType - Set the type of GLEE scheme

961:   Logically collective

963:   Input Parameter:
964: +  ts - timestepping context
965: -  gleetype - type of GLEE-scheme

967:   Level: intermediate

969: .seealso: TSGLEEGetType(), TSGLEE
970: @*/
971: PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype)
972: {

978:   PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));
979:   return(0);
980: }

982: /*@C
983:   TSGLEEGetType - Get the type of GLEE scheme

985:   Logically collective

987:   Input Parameter:
988: .  ts - timestepping context

990:   Output Parameter:
991: .  gleetype - type of GLEE-scheme

993:   Level: intermediate

995: .seealso: TSGLEESetType()
996: @*/
997: PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype)
998: {

1003:   PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));
1004:   return(0);
1005: }

1007: PetscErrorCode  TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype)
1008: {
1009:   TS_GLEE     *glee = (TS_GLEE*)ts->data;

1013:   if (!glee->tableau) {
1014:     TSGLEESetType(ts,TSGLEEDefaultType);
1015:   }
1016:   *gleetype = glee->tableau->name;
1017:   return(0);
1018: }
1019: PetscErrorCode  TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype)
1020: {
1021:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1022:   PetscErrorCode  ierr;
1023:   PetscBool       match;
1024:   GLEETableauLink link;

1027:   if (glee->tableau) {
1028:     PetscStrcmp(glee->tableau->name,gleetype,&match);
1029:     if (match) return(0);
1030:   }
1031:   for (link = GLEETableauList; link; link=link->next) {
1032:     PetscStrcmp(link->tab.name,gleetype,&match);
1033:     if (match) {
1034:       TSReset_GLEE(ts);
1035:       glee->tableau = &link->tab;
1036:       return(0);
1037:     }
1038:   }
1039:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype);
1040:   return(0);
1041: }

1043: static PetscErrorCode  TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y)
1044: {
1045:   TS_GLEE *glee = (TS_GLEE*)ts->data;

1048:   if (ns) *ns = glee->tableau->s;
1049:   if (Y)  *Y  = glee->YStage;
1050:   return(0);
1051: }

1053: PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y)
1054: {
1055:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1056:   GLEETableau     tab   = glee->tableau;
1057:   PetscErrorCode  ierr;

1060:   if (!Y) *n = tab->r;
1061:   else {
1062:     if ((*n >= 0) && (*n < tab->r)) {
1063:       VecCopy(glee->Y[*n],*Y);
1064:     } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1);
1065:   }
1066:   return(0);
1067: }

1069: PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X)
1070: {
1071:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1072:   GLEETableau     tab   = glee->tableau;
1073:   PetscReal       *F    = tab->Fembed;
1074:   PetscInt        r     = tab->r;
1075:   Vec             *Y    = glee->Y;
1076:   PetscScalar     *wr   = glee->rwork;
1077:   PetscInt        i;
1078:   PetscErrorCode  ierr;

1081:   VecZeroEntries(*X);
1082:   for (i=0; i<r; i++) wr[i] = F[i];
1083:   VecMAXPY((*X),r,wr,Y);
1084:   return(0);
1085: }

1087: PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X)
1088: {
1089:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1090:   GLEETableau     tab   = glee->tableau;
1091:   PetscReal       *F    = tab->Ferror;
1092:   PetscInt        r     = tab->r;
1093:   Vec             *Y    = glee->Y;
1094:   PetscScalar     *wr   = glee->rwork;
1095:   PetscInt        i;
1096:   PetscErrorCode  ierr;

1099:   VecZeroEntries(*X);
1100:   if(n==0){
1101:     for (i=0; i<r; i++) wr[i] = F[i];
1102:     VecMAXPY((*X),r,wr,Y);
1103:   } else if(n==-1) {
1104:     *X=glee->yGErr;
1105:   }
1106:   return(0);
1107: }

1109: PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X)
1110: {
1111:   TS_GLEE         *glee = (TS_GLEE*)ts->data;
1112:   GLEETableau     tab   = glee->tableau;
1113:   PetscReal       *S    = tab->Serror;
1114:   PetscInt        r     = tab->r,i;
1115:   Vec             *Y    = glee->Y;
1116:   PetscErrorCode  ierr;

1119:   if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r);
1120:   for (i=1; i<r; i++) {
1121:     VecCopy(ts->vec_sol,Y[i]);
1122:     VecAXPBY(Y[i],S[0],S[1],X);
1123:     VecCopy(X,glee->yGErr);
1124:   }
1125:   return(0);
1126: }

1128: static PetscErrorCode TSDestroy_GLEE(TS ts)
1129: {

1133:   TSReset_GLEE(ts);
1134:   if (ts->dm) {
1135:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
1136:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
1137:   }
1138:   PetscFree(ts->data);
1139:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);
1140:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);
1141:   return(0);
1142: }

1144: /* ------------------------------------------------------------ */
1145: /*MC
1146:       TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes

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

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

1154:   Level: beginner

1156: .seealso:  TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(),
1157:            TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 
1158:            TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister()

1160: M*/
1161: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
1162: {
1163:   TS_GLEE         *th;
1164:   PetscErrorCode  ierr;

1167: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
1168:   TSGLEEInitializePackage();
1169: #endif

1171:   ts->ops->reset                  = TSReset_GLEE;
1172:   ts->ops->destroy                = TSDestroy_GLEE;
1173:   ts->ops->view                   = TSView_GLEE;
1174:   ts->ops->load                   = TSLoad_GLEE;
1175:   ts->ops->setup                  = TSSetUp_GLEE;
1176:   ts->ops->step                   = TSStep_GLEE;
1177:   ts->ops->interpolate            = TSInterpolate_GLEE;
1178:   ts->ops->evaluatestep           = TSEvaluateStep_GLEE;
1179:   ts->ops->setfromoptions         = TSSetFromOptions_GLEE;
1180:   ts->ops->getstages              = TSGetStages_GLEE;
1181:   ts->ops->snesfunction           = SNESTSFormFunction_GLEE;
1182:   ts->ops->snesjacobian           = SNESTSFormJacobian_GLEE;
1183:   ts->ops->getsolutioncomponents  = TSGetSolutionComponents_GLEE;
1184:   ts->ops->getauxsolution         = TSGetAuxSolution_GLEE;
1185:   ts->ops->gettimeerror           = TSGetTimeError_GLEE;
1186:   ts->ops->settimeerror           = TSSetTimeError_GLEE;
1187:   ts->ops->startingmethod         = TSStartingMethod_GLEE;
1188:   ts->default_adapt_type          = TSADAPTGLEE;

1190:   ts->usessnes = PETSC_TRUE;

1192:   PetscNewLog(ts,&th);
1193:   ts->data = (void*)th;

1195:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);
1196:   PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);
1197:   return(0);
1198: }