Actual source code: eimex.c
petsc-3.6.1 2015-08-06
1: /*
2: * eimex.c
3: *
4: * Created on: Jun 21, 2012
5: * Written by Hong Zhang (zhang@vt.edu), Virginia Tech
6: * Emil Constantinescu (emconsta@mcs.anl.gov), Argonne National Laboratory
7: */
8: /*MC
9: EIMEX - Time stepping with Extrapolated IMEX methods.
11: Notes:
12: The general system is written as
14: G(t,X,Xdot) = F(t,X)
16: where G represents the stiff part and F represents the non-stiff part. The user should provide the stiff part
17: of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction().
18: This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian.
20: Another common form for the system is
22: y'=f(x)+g(x)
24: The relationship between F,G and f,g is
26: G = y'-g(x), F = f(x)
28: References
29: E. Constantinescu and A. Sandu, Extrapolated implicit-explicit time stepping, SIAM Journal on Scientific
30: Computing, 31 (2010), pp. 4452-4477.
32: Level: beginner
34: .seealso: TSCreate(), TS, TSSetType(), TSEIMEXSetMaxRows(), TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt()
36: M*/
37: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
38: #include <petscdm.h>
40: static const PetscInt TSEIMEXDefault = 3;
42: typedef struct {
43: PetscInt row_ind; /* Return the term T[row_ind][col_ind] */
44: PetscInt col_ind; /* Return the term T[row_ind][col_ind] */
45: PetscInt nstages; /* Numbers of stages in current scheme */
46: PetscInt max_rows; /* Maximum number of rows */
47: PetscInt *N; /* Harmonic sequence N[max_rows] */
48: Vec Y; /* States computed during the step, used to complete the step */
49: Vec Z; /* For shift*(Y-Z) */
50: Vec *T; /* Working table, size determined by nstages */
51: Vec YdotRHS; /* f(x) Work vector holding YdotRHS during residual evaluation */
52: Vec YdotI; /* xdot-g(x) Work vector holding YdotI = G(t,x,xdot) when xdot =0 */
53: Vec Ydot; /* f(x)+g(x) Work vector */
54: Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation) */
55: PetscReal shift;
56: PetscReal ctime;
57: PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */
58: PetscBool ord_adapt; /* order adapativity */
59: TSStepStatus status;
60: } TS_EIMEX;
62: /* This function is pure */
63: static PetscInt Map(PetscInt i, PetscInt j, PetscInt s)
64: {
65: return ((2*s-j+1)*j/2+i-j);
66: }
71: static PetscErrorCode TSEvaluateStep_EIMEX(TS ts,PetscInt order,Vec X,PetscBool *done)
72: {
73: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
74: const PetscInt ns = ext->nstages;
77: VecCopy(ext->T[Map(ext->row_ind,ext->col_ind,ns)],X);
78: return(0);
79: }
84: static PetscErrorCode TSStage_EIMEX(TS ts,PetscInt istage)
85: {
86: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
87: PetscReal h;
88: Vec Y=ext->Y, Z=ext->Z;
89: SNES snes;
90: TSAdapt adapt;
91: PetscInt i,its,lits;
92: PetscBool accept;
93: PetscErrorCode ierr;
96: TSGetSNES(ts,&snes);
97: h = ts->time_step/ext->N[istage];/* step size for the istage-th stage */
98: ext->shift = 1./h;
99: SNESSetLagJacobian(snes,-2); /* Recompute the Jacobian on this solve, but not again */
100: VecCopy(ext->VecSolPrev,Y); /* Take the previous solution as intial step */
102: for(i=0; i<ext->N[istage]; i++){
103: ext->ctime = ts->ptime + h*i;
104: VecCopy(Y,Z);/* Save the solution of the previous substep */
105: SNESSolve(snes,NULL,Y);
106: SNESGetIterationNumber(snes,&its);
107: SNESGetLinearSolveIterations(snes,&lits);
108: ts->snes_its += its; ts->ksp_its += lits;
109: TSGetAdapt(ts,&adapt);
110: TSAdaptCheckStage(adapt,ts,&accept);
111: }
113: return(0);
114: }
119: static PetscErrorCode TSStep_EIMEX(TS ts)
120: {
121: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
122: const PetscInt ns = ext->nstages;
123: Vec *T=ext->T, Y=ext->Y;
125: SNES snes;
126: PetscInt i,j;
127: PetscBool accept = PETSC_FALSE;
128: PetscErrorCode ierr;
129: PetscReal alpha,local_error;
132: TSGetSNES(ts,&snes);
133: SNESSetType(snes,"ksponly");
134: ext->status = TS_STEP_INCOMPLETE;
136: VecCopy(ts->vec_sol,ext->VecSolPrev);
138: /* Apply n_j steps of the base method to obtain solutions of T(j,1),1<=j<=s */
139: for(j=0; j<ns; j++){
140: TSStage_EIMEX(ts,j);
141: VecCopy(Y,T[j]);
142: }
144: for(i=1;i<ns;i++){
145: for(j=i;j<ns;j++){
146: alpha = -(PetscReal)ext->N[j]/ext->N[j-i];
147: 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] */
148: alpha = 1.0/(1.0 + alpha);
149: VecScale(T[Map(j,i,ns)],alpha);
150: }
151: }
153: TSEvaluateStep(ts,ns,ts->vec_sol,NULL);/*update ts solution */
155: if(ext->ord_adapt && ext->nstages < ext->max_rows){
156: accept = PETSC_FALSE;
157: while(!accept && ext->nstages < ext->max_rows){
158: TSErrorWeightedNorm(ts,ts->vec_sol,T[Map(ext->nstages-1,ext->nstages-2,ext->nstages)],ts->adapt->wnormtype,&local_error);
159: accept = (local_error < 1.0)? PETSC_TRUE : PETSC_FALSE;
161: if(!accept){/* add one more stage*/
162: TSStage_EIMEX(ts,ext->nstages);
163: ext->nstages++; ext->row_ind++; ext->col_ind++;
164: /*T table need to be recycled*/
165: VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);
166: for(i=0; i<ext->nstages-1; i++){
167: for(j=0; j<=i; j++){
168: VecCopy(T[Map(i,j,ext->nstages-1)],ext->T[Map(i,j,ext->nstages)]);
169: }
170: }
171: VecDestroyVecs(ext->nstages*(ext->nstages-1)/2,&T);
172: T = ext->T; /*reset the pointer*/
173: /*recycling finished, store the new solution*/
174: VecCopy(Y,T[ext->nstages-1]);
175: /*extrapolation for the newly added stage*/
176: for(i=1;i<ext->nstages;i++){
177: alpha = -(PetscReal)ext->N[ext->nstages-1]/ext->N[ext->nstages-1-i];
178: 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]*/
179: alpha = 1.0/(1.0 + alpha);
180: VecScale(T[Map(ext->nstages-1,i,ext->nstages)],alpha);
181: }
182: /*update ts solution */
183: TSEvaluateStep(ts,ext->nstages,ts->vec_sol,NULL);
184: }/*end if !accept*/
185: }/*end while*/
187: if(ext->nstages == ext->max_rows){
188: PetscInfo(ts,"Max number of rows has been used\n");
189: }
190: }/*end if ext->ord_adapt*/
191: ts->ptime += ts->time_step;
192: ts->steps++;
193: ext->status = TS_STEP_COMPLETE;
195: if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
196: return(0);
197: }
199: /* cubic Hermit spline */
202: static PetscErrorCode TSInterpolate_EIMEX(TS ts,PetscReal itime,Vec X)
203: {
204: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
205: PetscReal t,a,b;
206: Vec Y0=ext->VecSolPrev,Y1=ext->Y,
207: Ydot=ext->Ydot,YdotI=ext->YdotI;
208: const PetscReal h = ts->time_step_prev;
211: t = (itime -ts->ptime + h)/h;
212: /* YdotI = -f(x)-g(x) */
214: VecZeroEntries(Ydot);
215: TSComputeIFunction(ts,ts->ptime-h,Y0,Ydot,YdotI,PETSC_FALSE);
217: a = 2.0*t*t*t - 3.0*t*t + 1.0;
218: b = -(t*t*t - 2.0*t*t + t)*h;
219: VecAXPBYPCZ(X,a,b,0.0,Y0,YdotI);
221: TSComputeIFunction(ts,ts->ptime,Y1,Ydot,YdotI,PETSC_FALSE);
222: a = -2.0*t*t*t+3.0*t*t;
223: b = -(t*t*t - t*t)*h;
224: VecAXPBYPCZ(X,a,b,1.0,Y1,YdotI);
226: return(0);
227: }
232: static PetscErrorCode TSReset_EIMEX(TS ts)
233: {
234: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
235: PetscInt ns;
236: PetscErrorCode ierr;
239: ns = ext->nstages;
240: VecDestroyVecs((1+ns)*ns/2,&ext->T);
241: VecDestroy(&ext->Y);
242: VecDestroy(&ext->Z);
243: VecDestroy(&ext->YdotRHS);
244: VecDestroy(&ext->YdotI);
245: VecDestroy(&ext->Ydot);
246: VecDestroy(&ext->VecSolPrev);
247: PetscFree(ext->N);
248: return(0);
249: }
253: static PetscErrorCode TSDestroy_EIMEX(TS ts)
254: {
255: PetscErrorCode ierr;
258: TSReset_EIMEX(ts);
259: PetscFree(ts->data);
260: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C",NULL);
261: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C",NULL);
262: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",NULL);
264: return(0);
265: }
270: static PetscErrorCode TSEIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI, Vec *YdotRHS)
271: {
272: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
276: if (Z) {
277: if (dm && dm != ts->dm) {
278: DMGetNamedGlobalVector(dm,"TSEIMEX_Z",Z);
279: } else *Z = ext->Z;
280: }
281: if (Ydot) {
282: if (dm && dm != ts->dm) {
283: DMGetNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
284: } else *Ydot = ext->Ydot;
285: }
286: if (YdotI) {
287: if (dm && dm != ts->dm) {
288: DMGetNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
289: } else *YdotI = ext->YdotI;
290: }
291: if (YdotRHS) {
292: if (dm && dm != ts->dm) {
293: DMGetNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
294: } else *YdotRHS = ext->YdotRHS;
295: }
296: return(0);
297: }
302: static PetscErrorCode TSEIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI,Vec *YdotRHS)
303: {
307: if (Z) {
308: if (dm && dm != ts->dm) {
309: DMRestoreNamedGlobalVector(dm,"TSEIMEX_Z",Z);
310: }
311: }
312: if (Ydot) {
313: if (dm && dm != ts->dm) {
314: DMRestoreNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
315: }
316: }
317: if (YdotI) {
318: if (dm && dm != ts->dm) {
319: DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
320: }
321: }
322: if (YdotRHS) {
323: if (dm && dm != ts->dm) {
324: DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
325: }
326: }
327: return(0);
328: }
331: /*
332: This defines the nonlinear equation that is to be solved with SNES
333: Fn[t0+Theta*dt, U, (U-U0)*shift] = 0
334: In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U))
335: Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h
336: */
339: static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes,Vec X,Vec G,TS ts)
340: {
341: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
342: PetscErrorCode ierr;
343: Vec Ydot,Z;
344: DM dm,dmsave;
347: VecZeroEntries(G);
349: SNESGetDM(snes,&dm);
350: TSEIMEXGetVecs(ts,dm,&Z,&Ydot,NULL,NULL);
351: VecZeroEntries(Ydot);
352: dmsave = ts->dm;
353: ts->dm = dm;
354: TSComputeIFunction(ts,ext->ctime,X,Ydot,G,PETSC_FALSE);
355: /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function. */
356: VecCopy(G,Ydot);
357: ts->dm = dmsave;
358: TSEIMEXRestoreVecs(ts,dm,&Z,&Ydot,NULL,NULL);
360: return(0);
361: }
363: /*
364: This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y))
365: */
368: static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)
369: {
370: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
371: Vec Ydot;
372: PetscErrorCode ierr;
373: DM dm,dmsave;
375: SNESGetDM(snes,&dm);
376: TSEIMEXGetVecs(ts,dm,NULL,&Ydot,NULL,NULL);
377: /* VecZeroEntries(Ydot); */
378: /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */
379: dmsave = ts->dm;
380: ts->dm = dm;
381: TSComputeIJacobian(ts,ts->ptime,X,Ydot,ext->shift,A,B,PETSC_TRUE);
382: ts->dm = dmsave;
383: TSEIMEXRestoreVecs(ts,dm,NULL,&Ydot,NULL,NULL);
384: return(0);
385: }
389: static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine,DM coarse,void *ctx)
390: {
393: return(0);
394: }
398: static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
399: {
400: TS ts = (TS)ctx;
402: Vec Z,Z_c;
405: TSEIMEXGetVecs(ts,fine,&Z,NULL,NULL,NULL);
406: TSEIMEXGetVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
407: MatRestrict(restrct,Z,Z_c);
408: VecPointwiseMult(Z_c,rscale,Z_c);
409: TSEIMEXRestoreVecs(ts,fine,&Z,NULL,NULL,NULL);
410: TSEIMEXRestoreVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
411: return(0);
412: }
417: static PetscErrorCode TSSetUp_EIMEX(TS ts)
418: {
419: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
421: DM dm;
424: if (!ext->N){ /* ext->max_rows not set */
425: TSEIMEXSetMaxRows(ts,TSEIMEXDefault);
426: }
427: if(-1 == ext->row_ind && -1 == ext->col_ind){
428: TSEIMEXSetRowCol(ts,ext->max_rows,ext->max_rows);
429: } else{/* ext->row_ind and col_ind already set */
430: if (ext->ord_adapt){
431: PetscInfo(ts,"Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");
432: }
433: }
435: if(ext->ord_adapt){
436: ext->nstages = 2; /* Start with the 2-stage scheme */
437: TSEIMEXSetRowCol(ts,ext->nstages,ext->nstages);
438: } else{
439: ext->nstages = ext->max_rows; /* by default nstages is the same as max_rows, this can be changed by setting order adaptivity */
440: }
442: VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);/* full T table */
443: VecDuplicate(ts->vec_sol,&ext->YdotI);
444: VecDuplicate(ts->vec_sol,&ext->YdotRHS);
445: VecDuplicate(ts->vec_sol,&ext->Ydot);
446: VecDuplicate(ts->vec_sol,&ext->VecSolPrev);
447: VecDuplicate(ts->vec_sol,&ext->Y);
448: VecDuplicate(ts->vec_sol,&ext->Z);
449: TSGetDM(ts,&dm);
450: if (dm) {
451: DMCoarsenHookAdd(dm,DMCoarsenHook_TSEIMEX,DMRestrictHook_TSEIMEX,ts);
452: }
453: return(0);
454: }
458: static PetscErrorCode TSSetFromOptions_EIMEX(PetscOptions *PetscOptionsObject,TS ts)
459: {
460: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
462: PetscInt tindex[2];
463: PetscInt np = 2, nrows=TSEIMEXDefault;
466: tindex[0] = TSEIMEXDefault;
467: tindex[1] = TSEIMEXDefault;
468: PetscOptionsHead(PetscOptionsObject,"EIMEX ODE solver options");
469: {
470: PetscBool flg;
471: PetscOptionsInt("-ts_eimex_max_rows","Define the maximum number of rows used","TSEIMEXSetMaxRows",nrows,&nrows,&flg); /* default value 3 */
472: if(flg){
473: TSEIMEXSetMaxRows(ts,nrows);
474: }
475: PetscOptionsIntArray("-ts_eimex_row_col","Return the specific term in the T table","TSEIMEXSetRowCol",tindex,&np,&flg);
476: if(flg){
477: TSEIMEXSetRowCol(ts,tindex[0],tindex[1]);
478: }
479: PetscOptionsBool("-ts_eimex_order_adapt","Solve the problem with adaptive order","TSEIMEXSetOrdAdapt",ext->ord_adapt,&ext->ord_adapt,NULL);
480: SNESSetFromOptions(ts->snes);
481: }
482: PetscOptionsTail();
483: return(0);
484: }
488: static PetscErrorCode TSView_EIMEX(TS ts,PetscViewer viewer)
489: {
490: /* TS_EIMEX *ext = (TS_EIMEX*)ts->data; */
491: PetscBool iascii;
492: PetscErrorCode ierr;
495: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
496: if (iascii) {
497: PetscViewerASCIIPrintf(viewer," EIMEX\n");
498: }
499: SNESView(ts->snes,viewer);
500: return(0);
501: }
506: /*@C
507: TSEIMEXSetMaxRows - Set the maximum number of rows for EIMEX schemes
509: Logically collective
511: Input Parameter:
512: + ts - timestepping context
513: - nrows - maximum number of rows
515: Level: intermediate
517: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
518: @*/
519: PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows)
520: {
524: PetscTryMethod(ts,"TSEIMEXSetMaxRows_C",(TS,PetscInt),(ts,nrows));
525: return(0);
526: }
531: /*@C
532: TSEIMEXSetRowCol - Set the type index in the T table for the return value
534: Logically collective
536: Input Parameter:
537: + ts - timestepping context
538: - tindex - index in the T table
540: Level: intermediate
542: .seealso: TSEIMEXSetMaxRows(), TSEIMEXSetOrdAdapt(), TSEIMEX
543: @*/
544: PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col)
545: {
549: PetscTryMethod(ts,"TSEIMEXSetRowCol_C",(TS,PetscInt, PetscInt),(ts,row,col));
550: return(0);
551: }
556: /*@C
557: TSEIMEXSetOrdAdapt - Set the order adaptativity
559: Logically collective
561: Input Parameter:
562: + ts - timestepping context
563: - tindex - index in the T table
565: Level: intermediate
567: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
568: @*/
569: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
570: {
574: PetscTryMethod(ts,"TSEIMEXSetOrdAdapt_C",(TS,PetscBool),(ts,flg));
575: return(0);
576: }
581: static PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts,PetscInt nrows)
582: {
583: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
585: PetscInt i;
588: 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);
589: PetscFree(ext->N);
590: ext->max_rows = nrows;
591: PetscMalloc1(nrows,&ext->N);
592: for(i=0;i<nrows;i++) ext->N[i]=i+1;
593: return(0);
594: }
598: static PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts,PetscInt row,PetscInt col)
599: {
600: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
603: 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);
604: 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);
605: if (col > row) SETERRQ2(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The column index (%d) exceeds the row index (%d)\n",col,row);
607: ext->row_ind = row - 1;
608: ext->col_ind = col - 1; /* Array index in C starts from 0 */
609: return(0);
610: }
614: static PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts,PetscBool flg)
615: {
616: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
618: ext->ord_adapt = flg;
619: return(0);
620: }
622: /* ------------------------------------------------------------ */
623: /*MC
624: TSEIMEX - ODE solver using extrapolated IMEX schemes
625: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction().
627: Notes:
628: The default is a 3-stage scheme, it can be changed with TSEIMEXSetMaxRows() or -ts_eimex_max_rows
630: This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).
632: Level: beginner
634: .seealso: TSCreate(), TS
635: M*/
638: PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS ts)
639: {
640: TS_EIMEX *ext;
645: ts->ops->reset = TSReset_EIMEX;
646: ts->ops->destroy = TSDestroy_EIMEX;
647: ts->ops->view = TSView_EIMEX;
648: ts->ops->setup = TSSetUp_EIMEX;
649: ts->ops->step = TSStep_EIMEX;
650: ts->ops->interpolate = TSInterpolate_EIMEX;
651: ts->ops->evaluatestep = TSEvaluateStep_EIMEX;
652: ts->ops->setfromoptions = TSSetFromOptions_EIMEX;
653: ts->ops->snesfunction = SNESTSFormFunction_EIMEX;
654: ts->ops->snesjacobian = SNESTSFormJacobian_EIMEX;
656: PetscNewLog(ts,&ext);
657: ts->data = (void*)ext;
659: ext->ord_adapt = PETSC_FALSE; /* By default, no order adapativity */
660: ext->row_ind = -1;
661: ext->col_ind = -1;
662: ext->max_rows = TSEIMEXDefault;
663: ext->nstages = TSEIMEXDefault;
665: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C", TSEIMEXSetMaxRows_EIMEX);
666: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C", TSEIMEXSetRowCol_EIMEX);
667: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",TSEIMEXSetOrdAdapt_EIMEX);
668: return(0);
669: }