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[] = "@ARTICLE{Constantinescu_TR2016b,\n"
 15:                                " author = {Constantinescu, E.M.},\n"
 16:                                " title = {Estimating Global Errors in Time Stepping},\n"
 17:                                " journal = {ArXiv e-prints},\n"
 18:                                " year = 2016,\n"
 19:                                " adsurl = {http://adsabs.harvard.edu/abs/2015arXiv150305166C}\n}\n";

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

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

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

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

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

 70:      Level: advanced

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

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

 80:      Level: advanced

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

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

 90:      Level: advanced

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

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

100:      Level: advanced

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

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

110:      Level: advanced

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

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

120:      Level: advanced

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

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

130:      Level: advanced

132: .seealso: `TSGLEE`
133: M*/

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

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

140:   Level: advanced

142: .seealso: `TSGLEERegisterDestroy()`
143: @*/
144: PetscErrorCode TSGLEERegisterAll(void)
145: {
146:   if (TSGLEERegisterAllCalled) return 0;
147:   TSGLEERegisterAllCalled = PETSC_TRUE;

149:   {
150: #define GAMMA 0.5
151:     /* y-eps form */
152:     const PetscInt  p = 1, s = 3, r = 2;
153:     const PetscReal A[3][3] =
154:       {
155:         {1.0, 0,   0  },
156:         {0,   0.5, 0  },
157:         {0,   0.5, 0.5}
158:     },
159:                     B[2][3] = {{1.0, 0, 0}, {-2.0, 1.0, 1.0}}, U[3][2] = {{1.0, 0}, {1.0, 0.5}, {1.0, 0.5}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
160:     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);
161:   }
162:   {
163: #undef GAMMA
164: #define GAMMA 0.0
165:     /* y-eps form */
166:     const PetscInt  p = 2, s = 3, r = 2;
167:     const PetscReal A[3][3] =
168:       {
169:         {0,    0,    0},
170:         {1,    0,    0},
171:         {0.25, 0.25, 0}
172:     },
173:                     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}}, U[3][2] = {{1, 0}, {1, 10}, {1, -1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
174:     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);
175:   }
176:   {
177: #undef GAMMA
178: #define GAMMA 0.0
179:     /* y-y~ form */
180:     const PetscInt  p = 2, s = 4, r = 2;
181:     const PetscReal A[4][4] =
182:       {
183:         {0,            0,            0,            0},
184:         {0.75,         0,            0,            0},
185:         {0.25,         29.0 / 60.0,  0,            0},
186:         {-21.0 / 44.0, 145.0 / 44.0, -20.0 / 11.0, 0}
187:     },
188:                     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}}, U[4][2] = {{0, 1}, {75.0 / 58.0, -17.0 / 58.0}, {0, 1}, {0, 1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0};
189:     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);
190:   }
191:   {
192: #undef GAMMA
193: #define GAMMA 0.0
194:     /* y-y~ form */
195:     const PetscInt  p = 2, s = 5, r = 2;
196:     const PetscReal A[5][5] =
197:       {
198:         {0,                       0,                       0,                       0,                      0},
199:         {-0.94079244066783383269, 0,                       0,                       0,                      0},
200:         {0.64228187778301907108,  0.10915356933958500042,  0,                       0,                      0},
201:         {-0.51764297742287450812, 0.74414270351096040738,  -0.71404164927824538121, 0,                      0},
202:         {-0.44696561556825969206, -0.76768425657590196518, 0.20111608138142987881,  0.93828186737840469796, 0}
203:     },
204:                     B[2][5] = {{-0.029309178948150356153, -0.49671981884013874923, 0.34275801517650053274, 0.32941112623949194988, 0.85385985637229662276}, {0.78133219686062535272, 0.074238691892675897635, 0.57957363498384957966, -0.24638502829674959968, -0.18875949544040123033}}, U[5][2] = {{0.16911424754448327735, 0.83088575245551672265}, {0.53638465733199574340, 0.46361534266800425660}, {0.39901579167169582526, 0.60098420832830417474}, {0.87689005530618575480, 0.12310994469381424520}, {0.99056100455550913009, 0.0094389954444908699092}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0};
205:     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);
206:   }
207:   {
208: #undef GAMMA
209: #define GAMMA 0.0
210:     /* y-y~ form */
211:     const PetscInt  p = 3, s = 5, r = 2;
212:     const PetscReal A[5][5] =
213:       {
214:         {0,                                                0,                                                 0,                                                 0,                                               0},
215:         {-2169604947363702313.0 / 24313474998937147335.0,  0,                                                 0,                                                 0,                                               0},
216:         {46526746497697123895.0 / 94116917485856474137.0,  -10297879244026594958.0 / 49199457603717988219.0,  0,                                                 0,                                               0},
217:         {23364788935845982499.0 / 87425311444725389446.0,  -79205144337496116638.0 / 148994349441340815519.0, 40051189859317443782.0 / 36487615018004984309.0,   0,                                               0},
218:         {42089522664062539205.0 / 124911313006412840286.0, -15074384760342762939.0 / 137927286865289746282.0, -62274678522253371016.0 / 125918573676298591413.0, 13755475729852471739.0 / 79257927066651693390.0, 0}
219:     },
220:                     B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0, -55810892792806293355.0 / 206957624151308356511.0, 24061048952676379087.0 / 158739347956038723465.0, 3577972206874351339.0 / 7599733370677197135.0, -59449832954780563947.0 / 137360038685338563670.0}, {-9738262186984159168.0 / 99299082461487742983.0, -32797097931948613195.0 / 61521565616362163366.0, 42895514606418420631.0 / 71714201188501437336.0, 22608567633166065068.0 / 55371917805607957003.0, 94655809487476459565.0 / 151517167160302729021.0}}, U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0, 10043614439674808267.0 / 80863923579509469826.0}, {161694774978034105510.0 / 106187653640211060371.0, -55507121337823045139.0 / 106187653640211060371.0}, {78486094644566264568.0 / 88171030896733822981.0, 9684936252167558413.0 / 88171030896733822981.0}, {65394922146334854435.0 / 84570853840405479554.0, 19175931694070625119.0 / 84570853840405479554.0}, {8607282770183754108.0 / 108658046436496925911.0, 100050763666313171803.0 / 108658046436496925911.0}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0};
221:     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);
222:   }
223:   {
224: #undef GAMMA
225: #define GAMMA 0.25
226:     /* y-eps form */
227:     const PetscInt  p = 2, s = 6, r = 2;
228:     const PetscReal A[6][6] =
229:       {
230:         {0, 0, 0,    0,    0,   0},
231:         {1, 0, 0,    0,    0,   0},
232:         {0, 0, 0,    0,    0,   0},
233:         {0, 0, 0.5,  0,    0,   0},
234:         {0, 0, 0.25, 0.25, 0,   0},
235:         {0, 0, 0.25, 0.25, 0.5, 0}
236:     },
237:                     B[2][6] = {{0.5, 0.5, 0, 0, 0, 0}, {-2.0 / 3.0, -2.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}}, U[6][2] = {{1, 0}, {1, 0}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
238:     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);
239:   }
240:   {
241: #undef GAMMA
242: #define GAMMA 0.0
243:     /* y-eps form */
244:     const PetscInt  p = 3, s = 8, r = 2;
245:     const PetscReal A[8][8] =
246:       {
247:         {0,           0,          0,          0,          0,         0,         0,         0},
248:         {0.5,         0,          0,          0,          0,         0,         0,         0},
249:         {-1,          2,          0,          0,          0,         0,         0,         0},
250:         {1.0 / 6.0,   2.0 / 3.0,  1.0 / 6.0,  0,          0,         0,         0,         0},
251:         {0,           0,          0,          0,          0,         0,         0,         0},
252:         {-7.0 / 24.0, 1.0 / 3.0,  1.0 / 12.0, -1.0 / 8.0, 0.5,       0,         0,         0},
253:         {7.0 / 6.0,   -4.0 / 3.0, -1.0 / 3.0, 0.5,        -1.0,      2.0,       0,         0},
254:         {0,           0,          0,          0,          1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0}
255:     },
256:                     B[2][8] = {{1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0, 0, 0, 0, 0}, {-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}}, U[8][2] = {{1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 1}, {1, 1}, {1, 1}, {1, 1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
257:     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);
258:   }
259:   {
260: #undef GAMMA
261: #define GAMMA 0.25
262:     /* y-eps form */
263:     const PetscInt  p = 2, s = 9, r = 2;
264:     const PetscReal A[9][9] =
265:       {
266:         {0,                    0,                    0, 0,                     0,                    0,                    0,                     0,                    0},
267:         {0.585786437626904966, 0,                    0, 0,                     0,                    0,                    0,                     0,                    0},
268:         {0.149999999999999994, 0.849999999999999978, 0, 0,                     0,                    0,                    0,                     0,                    0},
269:         {0,                    0,                    0, 0,                     0,                    0,                    0,                     0,                    0},
270:         {0,                    0,                    0, 0.292893218813452483,  0,                    0,                    0,                     0,                    0},
271:         {0,                    0,                    0, 0.0749999999999999972, 0.424999999999999989, 0,                    0,                     0,                    0},
272:         {0,                    0,                    0, 0.176776695296636893,  0.176776695296636893, 0.146446609406726241, 0,                     0,                    0},
273:         {0,                    0,                    0, 0.176776695296636893,  0.176776695296636893, 0.146446609406726241, 0.292893218813452483,  0,                    0},
274:         {0,                    0,                    0, 0.176776695296636893,  0.176776695296636893, 0.146446609406726241, 0.0749999999999999972, 0.424999999999999989, 0}
275:     },
276:                     B[2][9] = {{0.353553390593273786, 0.353553390593273786, 0.292893218813452483, 0, 0, 0, 0, 0, 0}, {-0.471404520791031678, -0.471404520791031678, -0.390524291751269959, 0.235702260395515839, 0.235702260395515839, 0.195262145875634979, 0.235702260395515839, 0.235702260395515839, 0.195262145875634979}}, 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}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
277:     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);
278:   }

280:   return 0;
281: }

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

286:    Not Collective

288:    Level: advanced

290: .seealso: `TSGLEERegister()`, `TSGLEERegisterAll()`
291: @*/
292: PetscErrorCode TSGLEERegisterDestroy(void)
293: {
294:   GLEETableauLink link;

296:   while ((link = GLEETableauList)) {
297:     GLEETableau t   = &link->tab;
298:     GLEETableauList = link->next;
299:     PetscFree5(t->A, t->B, t->U, t->V, t->c);
300:     PetscFree2(t->S, t->F);
301:     PetscFree(t->Fembed);
302:     PetscFree(t->Ferror);
303:     PetscFree(t->Serror);
304:     PetscFree(t->binterp);
305:     PetscFree(t->name);
306:     PetscFree(link);
307:   }
308:   TSGLEERegisterAllCalled = PETSC_FALSE;
309:   return 0;
310: }

312: /*@C
313:   TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called
314:   from TSInitializePackage().

316:   Level: developer

318: .seealso: `PetscInitialize()`
319: @*/
320: PetscErrorCode TSGLEEInitializePackage(void)
321: {
322:   if (TSGLEEPackageInitialized) return 0;
323:   TSGLEEPackageInitialized = PETSC_TRUE;
324:   TSGLEERegisterAll();
325:   PetscObjectComposedDataRegister(&explicit_stage_time_id);
326:   PetscRegisterFinalize(TSGLEEFinalizePackage);
327:   return 0;
328: }

330: /*@C
331:   TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is
332:   called from PetscFinalize().

334:   Level: developer

336: .seealso: `PetscFinalize()`
337: @*/
338: PetscErrorCode TSGLEEFinalizePackage(void)
339: {
340:   TSGLEEPackageInitialized = PETSC_FALSE;
341:   TSGLEERegisterDestroy();
342:   return 0;
343: }

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

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

350:    Input Parameters:
351: +  name   - identifier for method
352: .  order  - order of method
353: .  s - number of stages
354: .  r - number of steps
355: .  gamma - LTE ratio
356: .  A - stage coefficients (dimension s*s, row-major)
357: .  B - step completion coefficients (dimension r*s, row-major)
358: .  U - method coefficients (dimension s*r, row-major)
359: .  V - method coefficients (dimension r*r, row-major)
360: .  S - starting coefficients
361: .  F - finishing coefficients
362: .  c - abscissa (dimension s; NULL to use row sums of A)
363: .  Fembed - step completion coefficients for embedded method
364: .  Ferror - error computation coefficients
365: .  Serror - error initialization coefficients
366: .  pinterp - order of interpolation (0 if unavailable)
367: -  binterp - array of interpolation coefficients (NULL if unavailable)

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

372:    Level: advanced

374: .seealso: `TSGLEE`
375: @*/
376: PetscErrorCode TSGLEERegister(TSGLEEType name, PetscInt order, PetscInt s, PetscInt r, PetscReal gamma, const PetscReal A[], const PetscReal B[], const PetscReal U[], const PetscReal V[], const PetscReal S[], const PetscReal F[], const PetscReal c[], const PetscReal Fembed[], const PetscReal Ferror[], const PetscReal Serror[], PetscInt pinterp, const PetscReal binterp[])
377: {
378:   GLEETableauLink link;
379:   GLEETableau     t;
380:   PetscInt        i, j;

382:   TSGLEEInitializePackage();
383:   PetscNew(&link);
384:   t = &link->tab;
385:   PetscStrallocpy(name, &t->name);
386:   t->order = order;
387:   t->s     = s;
388:   t->r     = r;
389:   t->gamma = gamma;
390:   PetscMalloc5(s * s, &t->A, r * r, &t->V, s, &t->c, r * s, &t->B, s * r, &t->U);
391:   PetscMalloc2(r, &t->S, r, &t->F);
392:   PetscArraycpy(t->A, A, s * s);
393:   PetscArraycpy(t->B, B, r * s);
394:   PetscArraycpy(t->U, U, s * r);
395:   PetscArraycpy(t->V, V, r * r);
396:   PetscArraycpy(t->S, S, r);
397:   PetscArraycpy(t->F, F, r);
398:   if (c) {
399:     PetscArraycpy(t->c, c, s);
400:   } else {
401:     for (i = 0; i < s; i++)
402:       for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j];
403:   }
404:   PetscMalloc1(r, &t->Fembed);
405:   PetscMalloc1(r, &t->Ferror);
406:   PetscMalloc1(r, &t->Serror);
407:   PetscArraycpy(t->Fembed, Fembed, r);
408:   PetscArraycpy(t->Ferror, Ferror, r);
409:   PetscArraycpy(t->Serror, Serror, r);
410:   t->pinterp = pinterp;
411:   PetscMalloc1(s * pinterp, &t->binterp);
412:   PetscArraycpy(t->binterp, binterp, s * pinterp);

414:   link->next      = GLEETableauList;
415:   GLEETableauList = link;
416:   return 0;
417: }

419: static PetscErrorCode TSEvaluateStep_GLEE(TS ts, PetscInt order, Vec X, PetscBool *done)
420: {
421:   TS_GLEE     *glee = (TS_GLEE *)ts->data;
422:   GLEETableau  tab  = glee->tableau;
423:   PetscReal    h, *B = tab->B, *V = tab->V, *F = tab->F, *Fembed = tab->Fembed;
424:   PetscInt     s = tab->s, r = tab->r, i, j;
425:   Vec         *Y = glee->Y, *YdotStage = glee->YdotStage;
426:   PetscScalar *ws = glee->swork, *wr = glee->rwork;


429:   switch (glee->status) {
430:   case TS_STEP_INCOMPLETE:
431:   case TS_STEP_PENDING:
432:     h = ts->time_step;
433:     break;
434:   case TS_STEP_COMPLETE:
435:     h = ts->ptime - ts->ptime_prev;
436:     break;
437:   default:
438:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
439:   }

441:   if (order == tab->order) {
442:     /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
443:              or TS_STEP_COMPLETE, glee->X has the solution at the
444:              beginning of the time step. So no need to roll-back.
445:     */
446:     if (glee->status == TS_STEP_INCOMPLETE) {
447:       for (i = 0; i < r; i++) {
448:         VecZeroEntries(Y[i]);
449:         for (j = 0; j < r; j++) wr[j] = V[i * r + j];
450:         VecMAXPY(Y[i], r, wr, glee->X);
451:         for (j = 0; j < s; j++) ws[j] = h * B[i * s + j];
452:         VecMAXPY(Y[i], s, ws, YdotStage);
453:       }
454:       VecZeroEntries(X);
455:       for (j = 0; j < r; j++) wr[j] = F[j];
456:       VecMAXPY(X, r, wr, Y);
457:     } else VecCopy(ts->vec_sol, X);
458:     return 0;

460:   } else if (order == tab->order - 1) {
461:     /* Complete with the embedded method (Fembed) */
462:     for (i = 0; i < r; i++) {
463:       VecZeroEntries(Y[i]);
464:       for (j = 0; j < r; j++) wr[j] = V[i * r + j];
465:       VecMAXPY(Y[i], r, wr, glee->X);
466:       for (j = 0; j < s; j++) ws[j] = h * B[i * s + j];
467:       VecMAXPY(Y[i], s, ws, YdotStage);
468:     }
469:     VecZeroEntries(X);
470:     for (j = 0; j < r; j++) wr[j] = Fembed[j];
471:     VecMAXPY(X, r, wr, Y);

473:     if (done) *done = PETSC_TRUE;
474:     return 0;
475:   }
476:   if (done) *done = PETSC_FALSE;
477:   else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "GLEE '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT, tab->name, tab->order, order);
478:   return 0;
479: }

481: static PetscErrorCode TSStep_GLEE(TS ts)
482: {
483:   TS_GLEE       *glee = (TS_GLEE *)ts->data;
484:   GLEETableau    tab  = glee->tableau;
485:   const PetscInt s = tab->s, r = tab->r;
486:   PetscReal     *A = tab->A, *U = tab->U, *F = tab->F, *c = tab->c;
487:   Vec           *Y = glee->Y, *X = glee->X, *YStage = glee->YStage, *YdotStage = glee->YdotStage, W = glee->W;
488:   SNES           snes;
489:   PetscScalar   *ws = glee->swork, *wr = glee->rwork;
490:   TSAdapt        adapt;
491:   PetscInt       i, j, reject, next_scheme, its, lits;
492:   PetscReal      next_time_step;
493:   PetscReal      t;
494:   PetscBool      accept;

496:   PetscCitationsRegister(citation, &cited);

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

500:   TSGetSNES(ts, &snes);
501:   next_time_step = ts->time_step;
502:   t              = ts->ptime;
503:   accept         = PETSC_TRUE;
504:   glee->status   = TS_STEP_INCOMPLETE;

506:   for (reject = 0; reject < ts->max_reject && !ts->reason; reject++, ts->reject++) {
507:     PetscReal h = ts->time_step;
508:     TSPreStep(ts);

510:     for (i = 0; i < s; i++) {
511:       glee->stage_time = t + h * c[i];
512:       TSPreStage(ts, glee->stage_time);

514:       if (A[i * s + i] == 0) { /* Explicit stage */
515:         VecZeroEntries(YStage[i]);
516:         for (j = 0; j < r; j++) wr[j] = U[i * r + j];
517:         VecMAXPY(YStage[i], r, wr, X);
518:         for (j = 0; j < i; j++) ws[j] = h * A[i * s + j];
519:         VecMAXPY(YStage[i], i, ws, YdotStage);
520:       } else { /* Implicit stage */
521:         glee->scoeff = 1.0 / A[i * s + i];
522:         /* compute right-hand-side */
523:         VecZeroEntries(W);
524:         for (j = 0; j < r; j++) wr[j] = U[i * r + j];
525:         VecMAXPY(W, r, wr, X);
526:         for (j = 0; j < i; j++) ws[j] = h * A[i * s + j];
527:         VecMAXPY(W, i, ws, YdotStage);
528:         VecScale(W, glee->scoeff / h);
529:         /* set initial guess */
530:         VecCopy(i > 0 ? YStage[i - 1] : ts->vec_sol, YStage[i]);
531:         /* solve for this stage */
532:         SNESSolve(snes, W, YStage[i]);
533:         SNESGetIterationNumber(snes, &its);
534:         SNESGetLinearSolveIterations(snes, &lits);
535:         ts->snes_its += its;
536:         ts->ksp_its += lits;
537:       }
538:       TSGetAdapt(ts, &adapt);
539:       TSAdaptCheckStage(adapt, ts, glee->stage_time, YStage[i], &accept);
540:       if (!accept) goto reject_step;
541:       TSPostStage(ts, glee->stage_time, i, YStage);
542:       TSComputeRHSFunction(ts, t + h * c[i], YStage[i], YdotStage[i]);
543:     }
544:     TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL);
545:     glee->status = TS_STEP_PENDING;

547:     /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
548:     TSGetAdapt(ts, &adapt);
549:     TSAdaptCandidatesClear(adapt);
550:     TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE);
551:     TSAdaptChoose(adapt, ts, ts->time_step, &next_scheme, &next_time_step, &accept);
552:     if (accept) {
553:       /* ignore next_scheme for now */
554:       ts->ptime += ts->time_step;
555:       ts->time_step = next_time_step;
556:       glee->status  = TS_STEP_COMPLETE;
557:       /* compute and store the global error */
558:       /* Note: this is not needed if TSAdaptGLEE is not used */
559:       TSGetTimeError(ts, 0, &(glee->yGErr));
560:       PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol, explicit_stage_time_id, ts->ptime);
561:       break;
562:     } else { /* Roll back the current step */
563:       for (j = 0; j < r; j++) wr[j] = F[j];
564:       VecMAXPY(ts->vec_sol, r, wr, X);
565:       ts->time_step = next_time_step;
566:       glee->status  = TS_STEP_INCOMPLETE;
567:     }
568:   reject_step:
569:     continue;
570:   }
571:   if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
572:   return 0;
573: }

575: static PetscErrorCode TSInterpolate_GLEE(TS ts, PetscReal itime, Vec X)
576: {
577:   TS_GLEE         *glee = (TS_GLEE *)ts->data;
578:   PetscInt         s = glee->tableau->s, pinterp = glee->tableau->pinterp, i, j;
579:   PetscReal        h, tt, t;
580:   PetscScalar     *b;
581:   const PetscReal *B = glee->tableau->binterp;

584:   switch (glee->status) {
585:   case TS_STEP_INCOMPLETE:
586:   case TS_STEP_PENDING:
587:     h = ts->time_step;
588:     t = (itime - ts->ptime) / h;
589:     break;
590:   case TS_STEP_COMPLETE:
591:     h = ts->ptime - ts->ptime_prev;
592:     t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
593:     break;
594:   default:
595:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
596:   }
597:   PetscMalloc1(s, &b);
598:   for (i = 0; i < s; i++) b[i] = 0;
599:   for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
600:     for (i = 0; i < s; i++) b[i] += h * B[i * pinterp + j] * tt;
601:   }
602:   VecCopy(glee->YStage[0], X);
603:   VecMAXPY(X, s, b, glee->YdotStage);
604:   PetscFree(b);
605:   return 0;
606: }

608: /*------------------------------------------------------------*/
609: static PetscErrorCode TSReset_GLEE(TS ts)
610: {
611:   TS_GLEE *glee = (TS_GLEE *)ts->data;
612:   PetscInt s, r;

614:   if (!glee->tableau) return 0;
615:   s = glee->tableau->s;
616:   r = glee->tableau->r;
617:   VecDestroyVecs(r, &glee->Y);
618:   VecDestroyVecs(r, &glee->X);
619:   VecDestroyVecs(s, &glee->YStage);
620:   VecDestroyVecs(s, &glee->YdotStage);
621:   VecDestroy(&glee->Ydot);
622:   VecDestroy(&glee->yGErr);
623:   VecDestroy(&glee->W);
624:   PetscFree2(glee->swork, glee->rwork);
625:   return 0;
626: }

628: static PetscErrorCode TSGLEEGetVecs(TS ts, DM dm, Vec *Ydot)
629: {
630:   TS_GLEE *glee = (TS_GLEE *)ts->data;

632:   if (Ydot) {
633:     if (dm && dm != ts->dm) {
634:       DMGetNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot);
635:     } else *Ydot = glee->Ydot;
636:   }
637:   return 0;
638: }

640: static PetscErrorCode TSGLEERestoreVecs(TS ts, DM dm, Vec *Ydot)
641: {
642:   if (Ydot) {
643:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot);
644:   }
645:   return 0;
646: }

648: /*
649:   This defines the nonlinear equation that is to be solved with SNES
650: */
651: static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes, Vec X, Vec F, TS ts)
652: {
653:   TS_GLEE  *glee = (TS_GLEE *)ts->data;
654:   DM        dm, dmsave;
655:   Vec       Ydot;
656:   PetscReal shift = glee->scoeff / ts->time_step;

658:   SNESGetDM(snes, &dm);
659:   TSGLEEGetVecs(ts, dm, &Ydot);
660:   /* Set Ydot = shift*X */
661:   VecCopy(X, Ydot);
662:   VecScale(Ydot, shift);
663:   dmsave = ts->dm;
664:   ts->dm = dm;

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

668:   ts->dm = dmsave;
669:   TSGLEERestoreVecs(ts, dm, &Ydot);
670:   return 0;
671: }

673: static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes, Vec X, Mat A, Mat B, TS ts)
674: {
675:   TS_GLEE  *glee = (TS_GLEE *)ts->data;
676:   DM        dm, dmsave;
677:   Vec       Ydot;
678:   PetscReal shift = glee->scoeff / ts->time_step;

680:   SNESGetDM(snes, &dm);
681:   TSGLEEGetVecs(ts, dm, &Ydot);
682:   /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
683:   dmsave = ts->dm;
684:   ts->dm = dm;

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

688:   ts->dm = dmsave;
689:   TSGLEERestoreVecs(ts, dm, &Ydot);
690:   return 0;
691: }

693: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine, DM coarse, void *ctx)
694: {
695:   return 0;
696: }

698: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
699: {
700:   return 0;
701: }

703: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm, DM subdm, void *ctx)
704: {
705:   return 0;
706: }

708: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
709: {
710:   return 0;
711: }

713: static PetscErrorCode TSSetUp_GLEE(TS ts)
714: {
715:   TS_GLEE    *glee = (TS_GLEE *)ts->data;
716:   GLEETableau tab;
717:   PetscInt    s, r;
718:   DM          dm;

720:   if (!glee->tableau) TSGLEESetType(ts, TSGLEEDefaultType);
721:   tab = glee->tableau;
722:   s   = tab->s;
723:   r   = tab->r;
724:   VecDuplicateVecs(ts->vec_sol, r, &glee->Y);
725:   VecDuplicateVecs(ts->vec_sol, r, &glee->X);
726:   VecDuplicateVecs(ts->vec_sol, s, &glee->YStage);
727:   VecDuplicateVecs(ts->vec_sol, s, &glee->YdotStage);
728:   VecDuplicate(ts->vec_sol, &glee->Ydot);
729:   VecDuplicate(ts->vec_sol, &glee->yGErr);
730:   VecZeroEntries(glee->yGErr);
731:   VecDuplicate(ts->vec_sol, &glee->W);
732:   PetscMalloc2(s, &glee->swork, r, &glee->rwork);
733:   TSGetDM(ts, &dm);
734:   DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts);
735:   DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts);
736:   return 0;
737: }

739: PetscErrorCode TSStartingMethod_GLEE(TS ts)
740: {
741:   TS_GLEE    *glee = (TS_GLEE *)ts->data;
742:   GLEETableau tab  = glee->tableau;
743:   PetscInt    r    = tab->r, i;
744:   PetscReal  *S    = tab->S;

746:   for (i = 0; i < r; i++) {
747:     VecZeroEntries(glee->Y[i]);
748:     VecAXPY(glee->Y[i], S[i], ts->vec_sol);
749:   }

751:   return 0;
752: }

754: /*------------------------------------------------------------*/

756: static PetscErrorCode TSSetFromOptions_GLEE(TS ts, PetscOptionItems *PetscOptionsObject)
757: {
758:   char gleetype[256];

760:   PetscOptionsHeadBegin(PetscOptionsObject, "GLEE ODE solver options");
761:   {
762:     GLEETableauLink link;
763:     PetscInt        count, choice;
764:     PetscBool       flg;
765:     const char    **namelist;

767:     PetscStrncpy(gleetype, TSGLEEDefaultType, sizeof(gleetype));
768:     for (link = GLEETableauList, count = 0; link; link = link->next, count++)
769:       ;
770:     PetscMalloc1(count, (char ***)&namelist);
771:     for (link = GLEETableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
772:     PetscOptionsEList("-ts_glee_type", "Family of GLEE method", "TSGLEESetType", (const char *const *)namelist, count, gleetype, &choice, &flg);
773:     TSGLEESetType(ts, flg ? namelist[choice] : gleetype);
774:     PetscFree(namelist);
775:   }
776:   PetscOptionsHeadEnd();
777:   return 0;
778: }

780: static PetscErrorCode TSView_GLEE(TS ts, PetscViewer viewer)
781: {
782:   TS_GLEE    *glee = (TS_GLEE *)ts->data;
783:   GLEETableau tab  = glee->tableau;
784:   PetscBool   iascii;

786:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
787:   if (iascii) {
788:     TSGLEEType gleetype;
789:     char       buf[512];
790:     TSGLEEGetType(ts, &gleetype);
791:     PetscViewerASCIIPrintf(viewer, "  GLEE type %s\n", gleetype);
792:     PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c);
793:     PetscViewerASCIIPrintf(viewer, "  Abscissa     c = %s\n", buf);
794:     /* Note: print out r as well */
795:   }
796:   return 0;
797: }

799: static PetscErrorCode TSLoad_GLEE(TS ts, PetscViewer viewer)
800: {
801:   SNES    snes;
802:   TSAdapt tsadapt;

804:   TSGetAdapt(ts, &tsadapt);
805:   TSAdaptLoad(tsadapt, viewer);
806:   TSGetSNES(ts, &snes);
807:   SNESLoad(snes, viewer);
808:   /* function and Jacobian context for SNES when used with TS is always ts object */
809:   SNESSetFunction(snes, NULL, NULL, ts);
810:   SNESSetJacobian(snes, NULL, NULL, NULL, ts);
811:   return 0;
812: }

814: /*@C
815:   TSGLEESetType - Set the type of GLEE scheme

817:   Logically collective

819:   Input Parameters:
820: +  ts - timestepping context
821: -  gleetype - type of GLEE-scheme

823:   Level: intermediate

825: .seealso: `TSGLEEGetType()`, `TSGLEE`
826: @*/
827: PetscErrorCode TSGLEESetType(TS ts, TSGLEEType gleetype)
828: {
831:   PetscTryMethod(ts, "TSGLEESetType_C", (TS, TSGLEEType), (ts, gleetype));
832:   return 0;
833: }

835: /*@C
836:   TSGLEEGetType - Get the type of GLEE scheme

838:   Logically collective

840:   Input Parameter:
841: .  ts - timestepping context

843:   Output Parameter:
844: .  gleetype - type of GLEE-scheme

846:   Level: intermediate

848: .seealso: `TSGLEESetType()`
849: @*/
850: PetscErrorCode TSGLEEGetType(TS ts, TSGLEEType *gleetype)
851: {
853:   PetscUseMethod(ts, "TSGLEEGetType_C", (TS, TSGLEEType *), (ts, gleetype));
854:   return 0;
855: }

857: PetscErrorCode TSGLEEGetType_GLEE(TS ts, TSGLEEType *gleetype)
858: {
859:   TS_GLEE *glee = (TS_GLEE *)ts->data;

861:   if (!glee->tableau) TSGLEESetType(ts, TSGLEEDefaultType);
862:   *gleetype = glee->tableau->name;
863:   return 0;
864: }
865: PetscErrorCode TSGLEESetType_GLEE(TS ts, TSGLEEType gleetype)
866: {
867:   TS_GLEE        *glee = (TS_GLEE *)ts->data;
868:   PetscBool       match;
869:   GLEETableauLink link;

871:   if (glee->tableau) {
872:     PetscStrcmp(glee->tableau->name, gleetype, &match);
873:     if (match) return 0;
874:   }
875:   for (link = GLEETableauList; link; link = link->next) {
876:     PetscStrcmp(link->tab.name, gleetype, &match);
877:     if (match) {
878:       TSReset_GLEE(ts);
879:       glee->tableau = &link->tab;
880:       return 0;
881:     }
882:   }
883:   SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", gleetype);
884: }

886: static PetscErrorCode TSGetStages_GLEE(TS ts, PetscInt *ns, Vec **Y)
887: {
888:   TS_GLEE *glee = (TS_GLEE *)ts->data;

890:   if (ns) *ns = glee->tableau->s;
891:   if (Y) *Y = glee->YStage;
892:   return 0;
893: }

895: PetscErrorCode TSGetSolutionComponents_GLEE(TS ts, PetscInt *n, Vec *Y)
896: {
897:   TS_GLEE    *glee = (TS_GLEE *)ts->data;
898:   GLEETableau tab  = glee->tableau;

900:   if (!Y) *n = tab->r;
901:   else {
902:     if ((*n >= 0) && (*n < tab->r)) {
903:       VecCopy(glee->Y[*n], *Y);
904:     } else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Second argument (%" PetscInt_FMT ") out of range[0,%" PetscInt_FMT "].", *n, tab->r - 1);
905:   }
906:   return 0;
907: }

909: PetscErrorCode TSGetAuxSolution_GLEE(TS ts, Vec *X)
910: {
911:   TS_GLEE     *glee = (TS_GLEE *)ts->data;
912:   GLEETableau  tab  = glee->tableau;
913:   PetscReal   *F    = tab->Fembed;
914:   PetscInt     r    = tab->r;
915:   Vec         *Y    = glee->Y;
916:   PetscScalar *wr   = glee->rwork;
917:   PetscInt     i;

919:   VecZeroEntries(*X);
920:   for (i = 0; i < r; i++) wr[i] = F[i];
921:   VecMAXPY((*X), r, wr, Y);
922:   return 0;
923: }

925: PetscErrorCode TSGetTimeError_GLEE(TS ts, PetscInt n, Vec *X)
926: {
927:   TS_GLEE     *glee = (TS_GLEE *)ts->data;
928:   GLEETableau  tab  = glee->tableau;
929:   PetscReal   *F    = tab->Ferror;
930:   PetscInt     r    = tab->r;
931:   Vec         *Y    = glee->Y;
932:   PetscScalar *wr   = glee->rwork;
933:   PetscInt     i;

935:   VecZeroEntries(*X);
936:   if (n == 0) {
937:     for (i = 0; i < r; i++) wr[i] = F[i];
938:     VecMAXPY((*X), r, wr, Y);
939:   } else if (n == -1) {
940:     *X = glee->yGErr;
941:   }
942:   return 0;
943: }

945: PetscErrorCode TSSetTimeError_GLEE(TS ts, Vec X)
946: {
947:   TS_GLEE    *glee = (TS_GLEE *)ts->data;
948:   GLEETableau tab  = glee->tableau;
949:   PetscReal  *S    = tab->Serror;
950:   PetscInt    r    = tab->r, i;
951:   Vec        *Y    = glee->Y;

954:   for (i = 1; i < r; i++) {
955:     VecCopy(ts->vec_sol, Y[i]);
956:     VecAXPBY(Y[i], S[0], S[1], X);
957:     VecCopy(X, glee->yGErr);
958:   }
959:   return 0;
960: }

962: static PetscErrorCode TSDestroy_GLEE(TS ts)
963: {
964:   TSReset_GLEE(ts);
965:   if (ts->dm) {
966:     DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts);
967:     DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts);
968:   }
969:   PetscFree(ts->data);
970:   PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", NULL);
971:   PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", NULL);
972:   return 0;
973: }

975: /* ------------------------------------------------------------ */
976: /*MC
977:       TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes

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

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

985:   Level: beginner

987: .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSGLEESetType()`, `TSGLEEGetType()`,
988:           `TSGLEE23`, `TTSGLEE24`, `TSGLEE35`, `TSGLEE25I`, `TSGLEEEXRK2A`,
989:           `TSGLEERK32G1`, `TSGLEERK285EX`, `TSGLEEType`, `TSGLEERegister()`

991: M*/
992: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
993: {
994:   TS_GLEE *th;

996:   TSGLEEInitializePackage();

998:   ts->ops->reset                 = TSReset_GLEE;
999:   ts->ops->destroy               = TSDestroy_GLEE;
1000:   ts->ops->view                  = TSView_GLEE;
1001:   ts->ops->load                  = TSLoad_GLEE;
1002:   ts->ops->setup                 = TSSetUp_GLEE;
1003:   ts->ops->step                  = TSStep_GLEE;
1004:   ts->ops->interpolate           = TSInterpolate_GLEE;
1005:   ts->ops->evaluatestep          = TSEvaluateStep_GLEE;
1006:   ts->ops->setfromoptions        = TSSetFromOptions_GLEE;
1007:   ts->ops->getstages             = TSGetStages_GLEE;
1008:   ts->ops->snesfunction          = SNESTSFormFunction_GLEE;
1009:   ts->ops->snesjacobian          = SNESTSFormJacobian_GLEE;
1010:   ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE;
1011:   ts->ops->getauxsolution        = TSGetAuxSolution_GLEE;
1012:   ts->ops->gettimeerror          = TSGetTimeError_GLEE;
1013:   ts->ops->settimeerror          = TSSetTimeError_GLEE;
1014:   ts->ops->startingmethod        = TSStartingMethod_GLEE;
1015:   ts->default_adapt_type         = TSADAPTGLEE;

1017:   ts->usessnes = PETSC_TRUE;

1019:   PetscNew(&th);
1020:   ts->data = (void *)th;

1022:   PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", TSGLEEGetType_GLEE);
1023:   PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", TSGLEESetType_GLEE);
1024:   return 0;
1025: }