Actual source code: eimex.c


  2: #include <petsc/private/tsimpl.h>
  3: #include <petscdm.h>

  5: static const PetscInt TSEIMEXDefault = 3;

  7: typedef struct {
  8:   PetscInt     row_ind;    /* Return the term T[row_ind][col_ind] */
  9:   PetscInt     col_ind;    /* Return the term T[row_ind][col_ind] */
 10:   PetscInt     nstages;    /* Numbers of stages in current scheme */
 11:   PetscInt     max_rows;   /* Maximum number of rows */
 12:   PetscInt    *N;          /* Harmonic sequence N[max_rows] */
 13:   Vec          Y;          /* States computed during the step, used to complete the step */
 14:   Vec          Z;          /* For shift*(Y-Z) */
 15:   Vec         *T;          /* Working table, size determined by nstages */
 16:   Vec          YdotRHS;    /* f(x) Work vector holding YdotRHS during residual evaluation */
 17:   Vec          YdotI;      /* xdot-g(x) Work vector holding YdotI = G(t,x,xdot) when xdot =0 */
 18:   Vec          Ydot;       /* f(x)+g(x) Work vector */
 19:   Vec          VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation) */
 20:   PetscReal    shift;
 21:   PetscReal    ctime;
 22:   PetscBool    recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */
 23:   PetscBool    ord_adapt;          /* order adapativity */
 24:   TSStepStatus status;
 25: } TS_EIMEX;

 27: /* This function is pure */
 28: static PetscInt Map(PetscInt i, PetscInt j, PetscInt s)
 29: {
 30:   return ((2 * s - j + 1) * j / 2 + i - j);
 31: }

 33: static PetscErrorCode TSEvaluateStep_EIMEX(TS ts, PetscInt order, Vec X, PetscBool *done)
 34: {
 35:   TS_EIMEX      *ext = (TS_EIMEX *)ts->data;
 36:   const PetscInt ns  = ext->nstages;
 37:   VecCopy(ext->T[Map(ext->row_ind, ext->col_ind, ns)], X);
 38:   return 0;
 39: }

 41: static PetscErrorCode TSStage_EIMEX(TS ts, PetscInt istage)
 42: {
 43:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
 44:   PetscReal h;
 45:   Vec       Y = ext->Y, Z = ext->Z;
 46:   SNES      snes;
 47:   TSAdapt   adapt;
 48:   PetscInt  i, its, lits;
 49:   PetscBool accept;

 51:   TSGetSNES(ts, &snes);
 52:   h          = ts->time_step / ext->N[istage]; /* step size for the istage-th stage */
 53:   ext->shift = 1. / h;
 54:   SNESSetLagJacobian(snes, -2); /* Recompute the Jacobian on this solve, but not again */
 55:   VecCopy(ext->VecSolPrev, Y);  /* Take the previous solution as initial step */

 57:   for (i = 0; i < ext->N[istage]; i++) {
 58:     ext->ctime = ts->ptime + h * i;
 59:     VecCopy(Y, Z); /* Save the solution of the previous substep */
 60:     SNESSolve(snes, NULL, Y);
 61:     SNESGetIterationNumber(snes, &its);
 62:     SNESGetLinearSolveIterations(snes, &lits);
 63:     ts->snes_its += its;
 64:     ts->ksp_its += lits;
 65:     TSGetAdapt(ts, &adapt);
 66:     TSAdaptCheckStage(adapt, ts, ext->ctime, Y, &accept);
 67:   }
 68:   return 0;
 69: }

 71: static PetscErrorCode TSStep_EIMEX(TS ts)
 72: {
 73:   TS_EIMEX      *ext = (TS_EIMEX *)ts->data;
 74:   const PetscInt ns  = ext->nstages;
 75:   Vec           *T = ext->T, Y = ext->Y;
 76:   SNES           snes;
 77:   PetscInt       i, j;
 78:   PetscBool      accept = PETSC_FALSE;
 79:   PetscReal      alpha, local_error, local_error_a, local_error_r;

 81:   TSGetSNES(ts, &snes);
 82:   SNESSetType(snes, "ksponly");
 83:   ext->status = TS_STEP_INCOMPLETE;

 85:   VecCopy(ts->vec_sol, ext->VecSolPrev);

 87:   /* Apply n_j steps of the base method to obtain solutions of T(j,1),1<=j<=s */
 88:   for (j = 0; j < ns; j++) {
 89:     TSStage_EIMEX(ts, j);
 90:     VecCopy(Y, T[j]);
 91:   }

 93:   for (i = 1; i < ns; i++) {
 94:     for (j = i; j < ns; j++) {
 95:       alpha = -(PetscReal)ext->N[j] / ext->N[j - i];
 96:       VecAXPBYPCZ(T[Map(j, i, ns)], alpha, 1.0, 0, T[Map(j, i - 1, ns)], T[Map(j - 1, i - 1, ns)]); /* T[j][i]=alpha*T[j][i-1]+T[j-1][i-1] */
 97:       alpha = 1.0 / (1.0 + alpha);
 98:       VecScale(T[Map(j, i, ns)], alpha);
 99:     }
100:   }

102:   TSEvaluateStep(ts, ns, ts->vec_sol, NULL); /*update ts solution */

104:   if (ext->ord_adapt && ext->nstages < ext->max_rows) {
105:     accept = PETSC_FALSE;
106:     while (!accept && ext->nstages < ext->max_rows) {
107:       TSErrorWeightedNorm(ts, ts->vec_sol, T[Map(ext->nstages - 1, ext->nstages - 2, ext->nstages)], ts->adapt->wnormtype, &local_error, &local_error_a, &local_error_r);
108:       accept = (local_error < 1.0) ? PETSC_TRUE : PETSC_FALSE;

110:       if (!accept) { /* add one more stage*/
111:         TSStage_EIMEX(ts, ext->nstages);
112:         ext->nstages++;
113:         ext->row_ind++;
114:         ext->col_ind++;
115:         /*T table need to be recycled*/
116:         VecDuplicateVecs(ts->vec_sol, (1 + ext->nstages) * ext->nstages / 2, &ext->T);
117:         for (i = 0; i < ext->nstages - 1; i++) {
118:           for (j = 0; j <= i; j++) VecCopy(T[Map(i, j, ext->nstages - 1)], ext->T[Map(i, j, ext->nstages)]);
119:         }
120:         VecDestroyVecs(ext->nstages * (ext->nstages - 1) / 2, &T);
121:         T = ext->T; /*reset the pointer*/
122:         /*recycling finished, store the new solution*/
123:         VecCopy(Y, T[ext->nstages - 1]);
124:         /*extrapolation for the newly added stage*/
125:         for (i = 1; i < ext->nstages; i++) {
126:           alpha = -(PetscReal)ext->N[ext->nstages - 1] / ext->N[ext->nstages - 1 - i];
127:           VecAXPBYPCZ(T[Map(ext->nstages - 1, i, ext->nstages)], alpha, 1.0, 0, T[Map(ext->nstages - 1, i - 1, ext->nstages)], T[Map(ext->nstages - 1 - 1, i - 1, ext->nstages)]); /*T[ext->nstages-1][i]=alpha*T[ext->nstages-1][i-1]+T[ext->nstages-1-1][i-1]*/
128:           alpha = 1.0 / (1.0 + alpha);
129:           VecScale(T[Map(ext->nstages - 1, i, ext->nstages)], alpha);
130:         }
131:         /*update ts solution */
132:         TSEvaluateStep(ts, ext->nstages, ts->vec_sol, NULL);
133:       } /*end if !accept*/
134:     }   /*end while*/

136:     if (ext->nstages == ext->max_rows) PetscInfo(ts, "Max number of rows has been used\n");
137:   } /*end if ext->ord_adapt*/
138:   ts->ptime += ts->time_step;
139:   ext->status = TS_STEP_COMPLETE;

141:   if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
142:   return 0;
143: }

145: /* cubic Hermit spline */
146: static PetscErrorCode TSInterpolate_EIMEX(TS ts, PetscReal itime, Vec X)
147: {
148:   TS_EIMEX       *ext = (TS_EIMEX *)ts->data;
149:   PetscReal       t, a, b;
150:   Vec             Y0 = ext->VecSolPrev, Y1 = ext->Y, Ydot = ext->Ydot, YdotI = ext->YdotI;
151:   const PetscReal h = ts->ptime - ts->ptime_prev;
152:   t = (itime - ts->ptime + h) / h;
153:   /* YdotI = -f(x)-g(x) */

155:   VecZeroEntries(Ydot);
156:   TSComputeIFunction(ts, ts->ptime - h, Y0, Ydot, YdotI, PETSC_FALSE);

158:   a = 2.0 * t * t * t - 3.0 * t * t + 1.0;
159:   b = -(t * t * t - 2.0 * t * t + t) * h;
160:   VecAXPBYPCZ(X, a, b, 0.0, Y0, YdotI);

162:   TSComputeIFunction(ts, ts->ptime, Y1, Ydot, YdotI, PETSC_FALSE);
163:   a = -2.0 * t * t * t + 3.0 * t * t;
164:   b = -(t * t * t - t * t) * h;
165:   VecAXPBYPCZ(X, a, b, 1.0, Y1, YdotI);

167:   return 0;
168: }

170: static PetscErrorCode TSReset_EIMEX(TS ts)
171: {
172:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
173:   PetscInt  ns;

175:   ns = ext->nstages;
176:   VecDestroyVecs((1 + ns) * ns / 2, &ext->T);
177:   VecDestroy(&ext->Y);
178:   VecDestroy(&ext->Z);
179:   VecDestroy(&ext->YdotRHS);
180:   VecDestroy(&ext->YdotI);
181:   VecDestroy(&ext->Ydot);
182:   VecDestroy(&ext->VecSolPrev);
183:   PetscFree(ext->N);
184:   return 0;
185: }

187: static PetscErrorCode TSDestroy_EIMEX(TS ts)
188: {
189:   TSReset_EIMEX(ts);
190:   PetscFree(ts->data);
191:   PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetMaxRows_C", NULL);
192:   PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetRowCol_C", NULL);
193:   PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetOrdAdapt_C", NULL);
194:   return 0;
195: }

197: static PetscErrorCode TSEIMEXGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydot, Vec *YdotI, Vec *YdotRHS)
198: {
199:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;

201:   if (Z) {
202:     if (dm && dm != ts->dm) {
203:       DMGetNamedGlobalVector(dm, "TSEIMEX_Z", Z);
204:     } else *Z = ext->Z;
205:   }
206:   if (Ydot) {
207:     if (dm && dm != ts->dm) {
208:       DMGetNamedGlobalVector(dm, "TSEIMEX_Ydot", Ydot);
209:     } else *Ydot = ext->Ydot;
210:   }
211:   if (YdotI) {
212:     if (dm && dm != ts->dm) {
213:       DMGetNamedGlobalVector(dm, "TSEIMEX_YdotI", YdotI);
214:     } else *YdotI = ext->YdotI;
215:   }
216:   if (YdotRHS) {
217:     if (dm && dm != ts->dm) {
218:       DMGetNamedGlobalVector(dm, "TSEIMEX_YdotRHS", YdotRHS);
219:     } else *YdotRHS = ext->YdotRHS;
220:   }
221:   return 0;
222: }

224: static PetscErrorCode TSEIMEXRestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydot, Vec *YdotI, Vec *YdotRHS)
225: {
226:   if (Z) {
227:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSEIMEX_Z", Z);
228:   }
229:   if (Ydot) {
230:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSEIMEX_Ydot", Ydot);
231:   }
232:   if (YdotI) {
233:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSEIMEX_YdotI", YdotI);
234:   }
235:   if (YdotRHS) {
236:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSEIMEX_YdotRHS", YdotRHS);
237:   }
238:   return 0;
239: }

241: /*
242:   This defines the nonlinear equation that is to be solved with SNES
243:   Fn[t0+Theta*dt, U, (U-U0)*shift] = 0
244:   In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U))
245:   Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h
246: */
247: static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes, Vec X, Vec G, TS ts)
248: {
249:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
250:   Vec       Ydot, Z;
251:   DM        dm, dmsave;

253:   VecZeroEntries(G);

255:   SNESGetDM(snes, &dm);
256:   TSEIMEXGetVecs(ts, dm, &Z, &Ydot, NULL, NULL);
257:   VecZeroEntries(Ydot);
258:   dmsave = ts->dm;
259:   ts->dm = dm;
260:   TSComputeIFunction(ts, ext->ctime, X, Ydot, G, PETSC_FALSE);
261:   /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function.  */
262:   VecCopy(G, Ydot);
263:   ts->dm = dmsave;
264:   TSEIMEXRestoreVecs(ts, dm, &Z, &Ydot, NULL, NULL);

266:   return 0;
267: }

269: /*
270:  This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y))
271:  */
272: static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes, Vec X, Mat A, Mat B, TS ts)
273: {
274:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
275:   Vec       Ydot;
276:   DM        dm, dmsave;
277:   SNESGetDM(snes, &dm);
278:   TSEIMEXGetVecs(ts, dm, NULL, &Ydot, NULL, NULL);
279:   /*  VecZeroEntries(Ydot); */
280:   /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */
281:   dmsave = ts->dm;
282:   ts->dm = dm;
283:   TSComputeIJacobian(ts, ts->ptime, X, Ydot, ext->shift, A, B, PETSC_TRUE);
284:   ts->dm = dmsave;
285:   TSEIMEXRestoreVecs(ts, dm, NULL, &Ydot, NULL, NULL);
286:   return 0;
287: }

289: static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine, DM coarse, void *ctx)
290: {
291:   return 0;
292: }

294: static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
295: {
296:   TS  ts = (TS)ctx;
297:   Vec Z, Z_c;

299:   TSEIMEXGetVecs(ts, fine, &Z, NULL, NULL, NULL);
300:   TSEIMEXGetVecs(ts, coarse, &Z_c, NULL, NULL, NULL);
301:   MatRestrict(restrct, Z, Z_c);
302:   VecPointwiseMult(Z_c, rscale, Z_c);
303:   TSEIMEXRestoreVecs(ts, fine, &Z, NULL, NULL, NULL);
304:   TSEIMEXRestoreVecs(ts, coarse, &Z_c, NULL, NULL, NULL);
305:   return 0;
306: }

308: static PetscErrorCode TSSetUp_EIMEX(TS ts)
309: {
310:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
311:   DM        dm;

313:   if (!ext->N) { /* ext->max_rows not set */
314:     TSEIMEXSetMaxRows(ts, TSEIMEXDefault);
315:   }
316:   if (-1 == ext->row_ind && -1 == ext->col_ind) {
317:     TSEIMEXSetRowCol(ts, ext->max_rows, ext->max_rows);
318:   } else { /* ext->row_ind and col_ind already set */
319:     if (ext->ord_adapt) PetscInfo(ts, "Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");
320:   }

322:   if (ext->ord_adapt) {
323:     ext->nstages = 2; /* Start with the 2-stage scheme */
324:     TSEIMEXSetRowCol(ts, ext->nstages, ext->nstages);
325:   } else {
326:     ext->nstages = ext->max_rows; /* by default nstages is the same as max_rows, this can be changed by setting order adaptivity */
327:   }

329:   TSGetAdapt(ts, &ts->adapt);

331:   VecDuplicateVecs(ts->vec_sol, (1 + ext->nstages) * ext->nstages / 2, &ext->T); /* full T table */
332:   VecDuplicate(ts->vec_sol, &ext->YdotI);
333:   VecDuplicate(ts->vec_sol, &ext->YdotRHS);
334:   VecDuplicate(ts->vec_sol, &ext->Ydot);
335:   VecDuplicate(ts->vec_sol, &ext->VecSolPrev);
336:   VecDuplicate(ts->vec_sol, &ext->Y);
337:   VecDuplicate(ts->vec_sol, &ext->Z);
338:   TSGetDM(ts, &dm);
339:   if (dm) DMCoarsenHookAdd(dm, DMCoarsenHook_TSEIMEX, DMRestrictHook_TSEIMEX, ts);
340:   return 0;
341: }

343: static PetscErrorCode TSSetFromOptions_EIMEX(TS ts, PetscOptionItems *PetscOptionsObject)
344: {
345:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
346:   PetscInt  tindex[2];
347:   PetscInt  np = 2, nrows = TSEIMEXDefault;

349:   tindex[0] = TSEIMEXDefault;
350:   tindex[1] = TSEIMEXDefault;
351:   PetscOptionsHeadBegin(PetscOptionsObject, "EIMEX ODE solver options");
352:   {
353:     PetscBool flg;
354:     PetscOptionsInt("-ts_eimex_max_rows", "Define the maximum number of rows used", "TSEIMEXSetMaxRows", nrows, &nrows, &flg); /* default value 3 */
355:     if (flg) TSEIMEXSetMaxRows(ts, nrows);
356:     PetscOptionsIntArray("-ts_eimex_row_col", "Return the specific term in the T table", "TSEIMEXSetRowCol", tindex, &np, &flg);
357:     if (flg) TSEIMEXSetRowCol(ts, tindex[0], tindex[1]);
358:     PetscOptionsBool("-ts_eimex_order_adapt", "Solve the problem with adaptive order", "TSEIMEXSetOrdAdapt", ext->ord_adapt, &ext->ord_adapt, NULL);
359:   }
360:   PetscOptionsHeadEnd();
361:   return 0;
362: }

364: static PetscErrorCode TSView_EIMEX(TS ts, PetscViewer viewer)
365: {
366:   return 0;
367: }

369: /*@C
370:   TSEIMEXSetMaxRows - Set the maximum number of rows for `TSEIMEX` schemes

372:   Logically collective

374:   Input Parameters:
375: +  ts - timestepping context
376: -  nrows - maximum number of rows

378:   Level: intermediate

380: .seealso: [](chapter_ts), `TSEIMEXSetRowCol()`, `TSEIMEXSetOrdAdapt()`, `TSEIMEX`
381: @*/
382: PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows)
383: {
385:   PetscTryMethod(ts, "TSEIMEXSetMaxRows_C", (TS, PetscInt), (ts, nrows));
386:   return 0;
387: }

389: /*@C
390:   TSEIMEXSetRowCol - Set the type index in the T table for the return value for the `TSEIMEX` scheme

392:   Logically collective

394:   Input Parameters:
395: +  ts - timestepping context
396: -  tindex - index in the T table

398:   Level: intermediate

400: .seealso: [](chapter_ts), `TSEIMEXSetMaxRows()`, `TSEIMEXSetOrdAdapt()`, `TSEIMEX`
401: @*/
402: PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col)
403: {
405:   PetscTryMethod(ts, "TSEIMEXSetRowCol_C", (TS, PetscInt, PetscInt), (ts, row, col));
406:   return 0;
407: }

409: /*@C
410:   TSEIMEXSetOrdAdapt - Set the order adaptativity for the `TSEIMEX` schemes

412:   Logically collective

414:   Input Parameters:
415: +  ts - timestepping context
416: -  tindex - index in the T table

418:   Level: intermediate

420: .seealso: [](chapter_ts), `TSEIMEXSetRowCol()`, `TSEIMEXSetOrdAdapt()`, `TSEIMEX`
421: @*/
422: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
423: {
425:   PetscTryMethod(ts, "TSEIMEXSetOrdAdapt_C", (TS, PetscBool), (ts, flg));
426:   return 0;
427: }

429: static PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts, PetscInt nrows)
430: {
431:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
432:   PetscInt  i;

435:   PetscFree(ext->N);
436:   ext->max_rows = nrows;
437:   PetscMalloc1(nrows, &ext->N);
438:   for (i = 0; i < nrows; i++) ext->N[i] = i + 1;
439:   return 0;
440: }

442: static PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts, PetscInt row, PetscInt col)
443: {
444:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;

448:              ext->max_rows);

451:   ext->row_ind = row - 1;
452:   ext->col_ind = col - 1; /* Array index in C starts from 0 */
453:   return 0;
454: }

456: static PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts, PetscBool flg)
457: {
458:   TS_EIMEX *ext = (TS_EIMEX *)ts->data;
459:   ext->ord_adapt = flg;
460:   return 0;
461: }

463: /*MC
464:       TSEIMEX - Time stepping with Extrapolated IMEX methods.

466:    These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly nonlinear such that it
467:    is expensive to solve with a fully implicit method. The user should provide the stiff part of the equation using `TSSetIFunction()` and the
468:    non-stiff part with `TSSetRHSFunction()`.

470:       Level: beginner

472:   Notes:
473:   The default is a 3-stage scheme, it can be changed with `TSEIMEXSetMaxRows()` or -ts_eimex_max_rows

475:   This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).

477:   The general system is written as

479:   G(t,X,Xdot) = F(t,X)

481:   where G represents the stiff part and F represents the non-stiff part. The user should provide the stiff part
482:   of the equation using TSSetIFunction() and the non-stiff part with `TSSetRHSFunction()`.
483:   This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian.

485:   Another common form for the system is

487:   y'=f(x)+g(x)

489:   The relationship between F,G and f,g is

491:   G = y'-g(x), F = f(x)

493:  Reference:
494: . [1] -  E. Constantinescu and A. Sandu, Extrapolated implicit-explicit time stepping, SIAM Journal on Scientific Computing, 31 (2010), pp. 4452-4477.

496: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSEIMEXSetMaxRows()`, `TSEIMEXSetRowCol()`, `TSEIMEXSetOrdAdapt()`, `TSType`
497:  M*/
498: PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS ts)
499: {
500:   TS_EIMEX *ext;


503:   ts->ops->reset          = TSReset_EIMEX;
504:   ts->ops->destroy        = TSDestroy_EIMEX;
505:   ts->ops->view           = TSView_EIMEX;
506:   ts->ops->setup          = TSSetUp_EIMEX;
507:   ts->ops->step           = TSStep_EIMEX;
508:   ts->ops->interpolate    = TSInterpolate_EIMEX;
509:   ts->ops->evaluatestep   = TSEvaluateStep_EIMEX;
510:   ts->ops->setfromoptions = TSSetFromOptions_EIMEX;
511:   ts->ops->snesfunction   = SNESTSFormFunction_EIMEX;
512:   ts->ops->snesjacobian   = SNESTSFormJacobian_EIMEX;
513:   ts->default_adapt_type  = TSADAPTNONE;

515:   ts->usessnes = PETSC_TRUE;

517:   PetscNew(&ext);
518:   ts->data = (void *)ext;

520:   ext->ord_adapt = PETSC_FALSE; /* By default, no order adapativity */
521:   ext->row_ind   = -1;
522:   ext->col_ind   = -1;
523:   ext->max_rows  = TSEIMEXDefault;
524:   ext->nstages   = TSEIMEXDefault;

526:   PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetMaxRows_C", TSEIMEXSetMaxRows_EIMEX);
527:   PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetRowCol_C", TSEIMEXSetRowCol_EIMEX);
528:   PetscObjectComposeFunction((PetscObject)ts, "TSEIMEXSetOrdAdapt_C", TSEIMEXSetOrdAdapt_EIMEX);
529:   return 0;
530: }