Actual source code: eimex.c
petsc-3.9.4 2018-09-11
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: }
34: static PetscErrorCode TSEvaluateStep_EIMEX(TS ts,PetscInt order,Vec X,PetscBool *done)
35: {
36: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
37: const PetscInt ns = ext->nstages;
40: VecCopy(ext->T[Map(ext->row_ind,ext->col_ind,ns)],X);
41: return(0);
42: }
45: static PetscErrorCode TSStage_EIMEX(TS ts,PetscInt istage)
46: {
47: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
48: PetscReal h;
49: Vec Y=ext->Y, Z=ext->Z;
50: SNES snes;
51: TSAdapt adapt;
52: PetscInt i,its,lits;
53: PetscBool accept;
54: PetscErrorCode ierr;
57: TSGetSNES(ts,&snes);
58: h = ts->time_step/ext->N[istage];/* step size for the istage-th stage */
59: ext->shift = 1./h;
60: SNESSetLagJacobian(snes,-2); /* Recompute the Jacobian on this solve, but not again */
61: VecCopy(ext->VecSolPrev,Y); /* Take the previous solution as intial step */
63: for(i=0; i<ext->N[istage]; i++){
64: ext->ctime = ts->ptime + h*i;
65: VecCopy(Y,Z);/* Save the solution of the previous substep */
66: SNESSolve(snes,NULL,Y);
67: SNESGetIterationNumber(snes,&its);
68: SNESGetLinearSolveIterations(snes,&lits);
69: ts->snes_its += its; ts->ksp_its += lits;
70: TSGetAdapt(ts,&adapt);
71: TSAdaptCheckStage(adapt,ts,ext->ctime,Y,&accept);
72: }
74: return(0);
75: }
78: static PetscErrorCode TSStep_EIMEX(TS ts)
79: {
80: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
81: const PetscInt ns = ext->nstages;
82: Vec *T=ext->T, Y=ext->Y;
84: SNES snes;
85: PetscInt i,j;
86: PetscBool accept = PETSC_FALSE;
87: PetscErrorCode ierr;
88: PetscReal alpha,local_error,local_error_a,local_error_r;
91: TSGetSNES(ts,&snes);
92: SNESSetType(snes,"ksponly");
93: ext->status = TS_STEP_INCOMPLETE;
95: VecCopy(ts->vec_sol,ext->VecSolPrev);
97: /* Apply n_j steps of the base method to obtain solutions of T(j,1),1<=j<=s */
98: for(j=0; j<ns; j++){
99: TSStage_EIMEX(ts,j);
100: VecCopy(Y,T[j]);
101: }
103: for(i=1;i<ns;i++){
104: for(j=i;j<ns;j++){
105: alpha = -(PetscReal)ext->N[j]/ext->N[j-i];
106: 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] */
107: alpha = 1.0/(1.0 + alpha);
108: VecScale(T[Map(j,i,ns)],alpha);
109: }
110: }
112: TSEvaluateStep(ts,ns,ts->vec_sol,NULL);/*update ts solution */
114: if(ext->ord_adapt && ext->nstages < ext->max_rows){
115: accept = PETSC_FALSE;
116: while(!accept && ext->nstages < ext->max_rows){
117: 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);
118: accept = (local_error < 1.0)? PETSC_TRUE : PETSC_FALSE;
120: if(!accept){/* add one more stage*/
121: TSStage_EIMEX(ts,ext->nstages);
122: ext->nstages++; ext->row_ind++; ext->col_ind++;
123: /*T table need to be recycled*/
124: VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);
125: for(i=0; i<ext->nstages-1; i++){
126: for(j=0; j<=i; j++){
127: VecCopy(T[Map(i,j,ext->nstages-1)],ext->T[Map(i,j,ext->nstages)]);
128: }
129: }
130: VecDestroyVecs(ext->nstages*(ext->nstages-1)/2,&T);
131: T = ext->T; /*reset the pointer*/
132: /*recycling finished, store the new solution*/
133: VecCopy(Y,T[ext->nstages-1]);
134: /*extrapolation for the newly added stage*/
135: for(i=1;i<ext->nstages;i++){
136: alpha = -(PetscReal)ext->N[ext->nstages-1]/ext->N[ext->nstages-1-i];
137: 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]*/
138: alpha = 1.0/(1.0 + alpha);
139: VecScale(T[Map(ext->nstages-1,i,ext->nstages)],alpha);
140: }
141: /*update ts solution */
142: TSEvaluateStep(ts,ext->nstages,ts->vec_sol,NULL);
143: }/*end if !accept*/
144: }/*end while*/
146: if(ext->nstages == ext->max_rows){
147: PetscInfo(ts,"Max number of rows has been used\n");
148: }
149: }/*end if ext->ord_adapt*/
150: ts->ptime += ts->time_step;
151: ext->status = TS_STEP_COMPLETE;
153: if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
154: return(0);
155: }
157: /* cubic Hermit spline */
158: static PetscErrorCode TSInterpolate_EIMEX(TS ts,PetscReal itime,Vec X)
159: {
160: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
161: PetscReal t,a,b;
162: Vec Y0=ext->VecSolPrev,Y1=ext->Y,Ydot=ext->Ydot,YdotI=ext->YdotI;
163: const PetscReal h = ts->ptime - ts->ptime_prev;
166: t = (itime -ts->ptime + h)/h;
167: /* YdotI = -f(x)-g(x) */
169: VecZeroEntries(Ydot);
170: TSComputeIFunction(ts,ts->ptime-h,Y0,Ydot,YdotI,PETSC_FALSE);
172: a = 2.0*t*t*t - 3.0*t*t + 1.0;
173: b = -(t*t*t - 2.0*t*t + t)*h;
174: VecAXPBYPCZ(X,a,b,0.0,Y0,YdotI);
176: TSComputeIFunction(ts,ts->ptime,Y1,Ydot,YdotI,PETSC_FALSE);
177: a = -2.0*t*t*t+3.0*t*t;
178: b = -(t*t*t - t*t)*h;
179: VecAXPBYPCZ(X,a,b,1.0,Y1,YdotI);
181: return(0);
182: }
185: static PetscErrorCode TSReset_EIMEX(TS ts)
186: {
187: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
188: PetscInt ns;
189: PetscErrorCode ierr;
192: ns = ext->nstages;
193: VecDestroyVecs((1+ns)*ns/2,&ext->T);
194: VecDestroy(&ext->Y);
195: VecDestroy(&ext->Z);
196: VecDestroy(&ext->YdotRHS);
197: VecDestroy(&ext->YdotI);
198: VecDestroy(&ext->Ydot);
199: VecDestroy(&ext->VecSolPrev);
200: PetscFree(ext->N);
201: return(0);
202: }
204: static PetscErrorCode TSDestroy_EIMEX(TS ts)
205: {
206: PetscErrorCode ierr;
209: TSReset_EIMEX(ts);
210: PetscFree(ts->data);
211: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C",NULL);
212: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C",NULL);
213: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",NULL);
215: return(0);
216: }
219: static PetscErrorCode TSEIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI, Vec *YdotRHS)
220: {
221: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
225: if (Z) {
226: if (dm && dm != ts->dm) {
227: DMGetNamedGlobalVector(dm,"TSEIMEX_Z",Z);
228: } else *Z = ext->Z;
229: }
230: if (Ydot) {
231: if (dm && dm != ts->dm) {
232: DMGetNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
233: } else *Ydot = ext->Ydot;
234: }
235: if (YdotI) {
236: if (dm && dm != ts->dm) {
237: DMGetNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
238: } else *YdotI = ext->YdotI;
239: }
240: if (YdotRHS) {
241: if (dm && dm != ts->dm) {
242: DMGetNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
243: } else *YdotRHS = ext->YdotRHS;
244: }
245: return(0);
246: }
249: static PetscErrorCode TSEIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI,Vec *YdotRHS)
250: {
254: if (Z) {
255: if (dm && dm != ts->dm) {
256: DMRestoreNamedGlobalVector(dm,"TSEIMEX_Z",Z);
257: }
258: }
259: if (Ydot) {
260: if (dm && dm != ts->dm) {
261: DMRestoreNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
262: }
263: }
264: if (YdotI) {
265: if (dm && dm != ts->dm) {
266: DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
267: }
268: }
269: if (YdotRHS) {
270: if (dm && dm != ts->dm) {
271: DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
272: }
273: }
274: return(0);
275: }
278: /*
279: This defines the nonlinear equation that is to be solved with SNES
280: Fn[t0+Theta*dt, U, (U-U0)*shift] = 0
281: In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U))
282: Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h
283: */
284: static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes,Vec X,Vec G,TS ts)
285: {
286: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
287: PetscErrorCode ierr;
288: Vec Ydot,Z;
289: DM dm,dmsave;
292: VecZeroEntries(G);
294: SNESGetDM(snes,&dm);
295: TSEIMEXGetVecs(ts,dm,&Z,&Ydot,NULL,NULL);
296: VecZeroEntries(Ydot);
297: dmsave = ts->dm;
298: ts->dm = dm;
299: TSComputeIFunction(ts,ext->ctime,X,Ydot,G,PETSC_FALSE);
300: /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function. */
301: VecCopy(G,Ydot);
302: ts->dm = dmsave;
303: TSEIMEXRestoreVecs(ts,dm,&Z,&Ydot,NULL,NULL);
305: return(0);
306: }
308: /*
309: This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y))
310: */
311: static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)
312: {
313: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
314: Vec Ydot;
315: PetscErrorCode ierr;
316: DM dm,dmsave;
318: SNESGetDM(snes,&dm);
319: TSEIMEXGetVecs(ts,dm,NULL,&Ydot,NULL,NULL);
320: /* VecZeroEntries(Ydot); */
321: /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */
322: dmsave = ts->dm;
323: ts->dm = dm;
324: TSComputeIJacobian(ts,ts->ptime,X,Ydot,ext->shift,A,B,PETSC_TRUE);
325: ts->dm = dmsave;
326: TSEIMEXRestoreVecs(ts,dm,NULL,&Ydot,NULL,NULL);
327: return(0);
328: }
330: static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine,DM coarse,void *ctx)
331: {
334: return(0);
335: }
337: static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
338: {
339: TS ts = (TS)ctx;
341: Vec Z,Z_c;
344: TSEIMEXGetVecs(ts,fine,&Z,NULL,NULL,NULL);
345: TSEIMEXGetVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
346: MatRestrict(restrct,Z,Z_c);
347: VecPointwiseMult(Z_c,rscale,Z_c);
348: TSEIMEXRestoreVecs(ts,fine,&Z,NULL,NULL,NULL);
349: TSEIMEXRestoreVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
350: return(0);
351: }
354: static PetscErrorCode TSSetUp_EIMEX(TS ts)
355: {
356: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
358: DM dm;
361: if (!ext->N){ /* ext->max_rows not set */
362: TSEIMEXSetMaxRows(ts,TSEIMEXDefault);
363: }
364: if(-1 == ext->row_ind && -1 == ext->col_ind){
365: TSEIMEXSetRowCol(ts,ext->max_rows,ext->max_rows);
366: } else{/* ext->row_ind and col_ind already set */
367: if (ext->ord_adapt){
368: PetscInfo(ts,"Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");
369: }
370: }
372: if(ext->ord_adapt){
373: ext->nstages = 2; /* Start with the 2-stage scheme */
374: TSEIMEXSetRowCol(ts,ext->nstages,ext->nstages);
375: } else{
376: ext->nstages = ext->max_rows; /* by default nstages is the same as max_rows, this can be changed by setting order adaptivity */
377: }
379: TSGetAdapt(ts,&ts->adapt);
381: VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);/* full T table */
382: VecDuplicate(ts->vec_sol,&ext->YdotI);
383: VecDuplicate(ts->vec_sol,&ext->YdotRHS);
384: VecDuplicate(ts->vec_sol,&ext->Ydot);
385: VecDuplicate(ts->vec_sol,&ext->VecSolPrev);
386: VecDuplicate(ts->vec_sol,&ext->Y);
387: VecDuplicate(ts->vec_sol,&ext->Z);
388: TSGetDM(ts,&dm);
389: if (dm) {
390: DMCoarsenHookAdd(dm,DMCoarsenHook_TSEIMEX,DMRestrictHook_TSEIMEX,ts);
391: }
392: return(0);
393: }
395: static PetscErrorCode TSSetFromOptions_EIMEX(PetscOptionItems *PetscOptionsObject,TS ts)
396: {
397: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
399: PetscInt tindex[2];
400: PetscInt np = 2, nrows=TSEIMEXDefault;
403: tindex[0] = TSEIMEXDefault;
404: tindex[1] = TSEIMEXDefault;
405: PetscOptionsHead(PetscOptionsObject,"EIMEX ODE solver options");
406: {
407: PetscBool flg;
408: PetscOptionsInt("-ts_eimex_max_rows","Define the maximum number of rows used","TSEIMEXSetMaxRows",nrows,&nrows,&flg); /* default value 3 */
409: if(flg){
410: TSEIMEXSetMaxRows(ts,nrows);
411: }
412: PetscOptionsIntArray("-ts_eimex_row_col","Return the specific term in the T table","TSEIMEXSetRowCol",tindex,&np,&flg);
413: if(flg){
414: TSEIMEXSetRowCol(ts,tindex[0],tindex[1]);
415: }
416: PetscOptionsBool("-ts_eimex_order_adapt","Solve the problem with adaptive order","TSEIMEXSetOrdAdapt",ext->ord_adapt,&ext->ord_adapt,NULL);
417: }
418: PetscOptionsTail();
419: return(0);
420: }
422: static PetscErrorCode TSView_EIMEX(TS ts,PetscViewer viewer)
423: {
425: return(0);
426: }
429: /*@C
430: TSEIMEXSetMaxRows - Set the maximum number of rows for EIMEX schemes
432: Logically collective
434: Input Parameter:
435: + ts - timestepping context
436: - nrows - maximum number of rows
438: Level: intermediate
440: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
441: @*/
442: PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows)
443: {
447: PetscTryMethod(ts,"TSEIMEXSetMaxRows_C",(TS,PetscInt),(ts,nrows));
448: return(0);
449: }
452: /*@C
453: TSEIMEXSetRowCol - Set the type index in the T table for the return value
455: Logically collective
457: Input Parameter:
458: + ts - timestepping context
459: - tindex - index in the T table
461: Level: intermediate
463: .seealso: TSEIMEXSetMaxRows(), TSEIMEXSetOrdAdapt(), TSEIMEX
464: @*/
465: PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col)
466: {
470: PetscTryMethod(ts,"TSEIMEXSetRowCol_C",(TS,PetscInt, PetscInt),(ts,row,col));
471: return(0);
472: }
475: /*@C
476: TSEIMEXSetOrdAdapt - Set the order adaptativity
478: Logically collective
480: Input Parameter:
481: + ts - timestepping context
482: - tindex - index in the T table
484: Level: intermediate
486: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
487: @*/
488: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
489: {
493: PetscTryMethod(ts,"TSEIMEXSetOrdAdapt_C",(TS,PetscBool),(ts,flg));
494: return(0);
495: }
498: static PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts,PetscInt nrows)
499: {
500: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
502: PetscInt i;
505: if (nrows < 0 || nrows > 100) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Max number of rows (current value %D) should be an integer number between 1 and 100\n",nrows);
506: PetscFree(ext->N);
507: ext->max_rows = nrows;
508: PetscMalloc1(nrows,&ext->N);
509: for(i=0;i<nrows;i++) ext->N[i]=i+1;
510: return(0);
511: }
513: static PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts,PetscInt row,PetscInt col)
514: {
515: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
518: if (row < 1 || col < 1) SETERRQ2(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The row or column index (current value %d,%d) should not be less than 1 \n",row,col);
519: if (row > ext->max_rows || col > ext->max_rows) SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The row or column index (current value %d,%d) exceeds the maximum number of rows %d\n",row,col,ext->max_rows);
520: if (col > row) SETERRQ2(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The column index (%d) exceeds the row index (%d)\n",col,row);
522: ext->row_ind = row - 1;
523: ext->col_ind = col - 1; /* Array index in C starts from 0 */
524: return(0);
525: }
527: static PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts,PetscBool flg)
528: {
529: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
531: ext->ord_adapt = flg;
532: return(0);
533: }
535: /*MC
536: TSEIMEX - Time stepping with Extrapolated IMEX methods.
538: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly nonlinear such that it
539: is expensive to solve with a fully implicit method. The user should provide the stiff part of the equation using TSSetIFunction() and the
540: non-stiff part with TSSetRHSFunction().
542: Notes:
543: The default is a 3-stage scheme, it can be changed with TSEIMEXSetMaxRows() or -ts_eimex_max_rows
545: This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).
547: The general system is written as
549: G(t,X,Xdot) = F(t,X)
551: where G represents the stiff part and F represents the non-stiff part. The user should provide the stiff part
552: of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction().
553: This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian.
555: Another common form for the system is
557: y'=f(x)+g(x)
559: The relationship between F,G and f,g is
561: G = y'-g(x), F = f(x)
563: References
564: E. Constantinescu and A. Sandu, Extrapolated implicit-explicit time stepping, SIAM Journal on Scientific
565: Computing, 31 (2010), pp. 4452-4477.
567: Level: beginner
569: .seealso: TSCreate(), TS, TSSetType(), TSEIMEXSetMaxRows(), TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt()
571: M*/
572: PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS ts)
573: {
574: TS_EIMEX *ext;
579: ts->ops->reset = TSReset_EIMEX;
580: ts->ops->destroy = TSDestroy_EIMEX;
581: ts->ops->view = TSView_EIMEX;
582: ts->ops->setup = TSSetUp_EIMEX;
583: ts->ops->step = TSStep_EIMEX;
584: ts->ops->interpolate = TSInterpolate_EIMEX;
585: ts->ops->evaluatestep = TSEvaluateStep_EIMEX;
586: ts->ops->setfromoptions = TSSetFromOptions_EIMEX;
587: ts->ops->snesfunction = SNESTSFormFunction_EIMEX;
588: ts->ops->snesjacobian = SNESTSFormJacobian_EIMEX;
589: ts->default_adapt_type = TSADAPTNONE;
591: ts->usessnes = PETSC_TRUE;
593: PetscNewLog(ts,&ext);
594: ts->data = (void*)ext;
596: ext->ord_adapt = PETSC_FALSE; /* By default, no order adapativity */
597: ext->row_ind = -1;
598: ext->col_ind = -1;
599: ext->max_rows = TSEIMEXDefault;
600: ext->nstages = TSEIMEXDefault;
602: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C", TSEIMEXSetMaxRows_EIMEX);
603: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C", TSEIMEXSetRowCol_EIMEX);
604: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",TSEIMEXSetOrdAdapt_EIMEX);
605: return(0);
606: }