Actual source code: eimex.c

petsc-master 2020-10-20
Report Typos and Errors

  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;
 39:   VecCopy(ext->T[Map(ext->row_ind,ext->col_ind,ns)],X);
 40:   return(0);
 41: }

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

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

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

 74: static PetscErrorCode TSStep_EIMEX(TS ts)
 75: {
 76:   TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
 77:   const PetscInt  ns = ext->nstages;
 78:   Vec             *T=ext->T, Y=ext->Y;

 80:   SNES            snes;
 81:   PetscInt        i,j;
 82:   PetscBool       accept = PETSC_FALSE;
 83:   PetscErrorCode  ierr;
 84:   PetscReal       alpha,local_error,local_error_a,local_error_r;

 87:   TSGetSNES(ts,&snes);
 88:   SNESSetType(snes,"ksponly");
 89:   ext->status = TS_STEP_INCOMPLETE;

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

 93:   /* Apply n_j steps of the base method to obtain solutions of T(j,1),1<=j<=s */
 94:   for (j=0; j<ns; j++){
 95:         TSStage_EIMEX(ts,j);
 96:         VecCopy(Y,T[j]);
 97:   }

 99:   for (i=1;i<ns;i++){
100:     for (j=i;j<ns;j++){
101:       alpha = -(PetscReal)ext->N[j]/ext->N[j-i];
102:       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] */
103:       alpha = 1.0/(1.0 + alpha);
104:       VecScale(T[Map(j,i,ns)],alpha);
105:     }
106:   }

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

110:   if (ext->ord_adapt && ext->nstages < ext->max_rows){
111:         accept = PETSC_FALSE;
112:         while (!accept && ext->nstages < ext->max_rows){
113:           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);
114:           accept = (local_error < 1.0)? PETSC_TRUE : PETSC_FALSE;

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

142:         if (ext->nstages == ext->max_rows){
143:           PetscInfo(ts,"Max number of rows has been used\n");
144:         }
145:   }/*end if ext->ord_adapt*/
146:   ts->ptime += ts->time_step;
147:   ext->status = TS_STEP_COMPLETE;

149:   if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
150:   return(0);
151: }

153: /* cubic Hermit spline */
154: static PetscErrorCode TSInterpolate_EIMEX(TS ts,PetscReal itime,Vec X)
155: {
156:   TS_EIMEX       *ext = (TS_EIMEX*)ts->data;
157:   PetscReal      t,a,b;
158:   Vec            Y0=ext->VecSolPrev,Y1=ext->Y,Ydot=ext->Ydot,YdotI=ext->YdotI;
159:   const PetscReal h = ts->ptime - ts->ptime_prev;
162:   t = (itime -ts->ptime + h)/h;
163:   /* YdotI = -f(x)-g(x) */

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

168:   a    = 2.0*t*t*t - 3.0*t*t + 1.0;
169:   b    = -(t*t*t - 2.0*t*t + t)*h;
170:   VecAXPBYPCZ(X,a,b,0.0,Y0,YdotI);

172:   TSComputeIFunction(ts,ts->ptime,Y1,Ydot,YdotI,PETSC_FALSE);
173:   a    = -2.0*t*t*t+3.0*t*t;
174:   b    = -(t*t*t - t*t)*h;
175:   VecAXPBYPCZ(X,a,b,1.0,Y1,YdotI);

177:   return(0);
178: }

180: static PetscErrorCode TSReset_EIMEX(TS ts)
181: {
182:   TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
183:   PetscInt        ns;
184:   PetscErrorCode  ierr;

187:   ns = ext->nstages;
188:   VecDestroyVecs((1+ns)*ns/2,&ext->T);
189:   VecDestroy(&ext->Y);
190:   VecDestroy(&ext->Z);
191:   VecDestroy(&ext->YdotRHS);
192:   VecDestroy(&ext->YdotI);
193:   VecDestroy(&ext->Ydot);
194:   VecDestroy(&ext->VecSolPrev);
195:   PetscFree(ext->N);
196:   return(0);
197: }

199: static PetscErrorCode TSDestroy_EIMEX(TS ts)
200: {
201:   PetscErrorCode  ierr;

204:   TSReset_EIMEX(ts);
205:   PetscFree(ts->data);
206:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C",NULL);
207:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C",NULL);
208:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",NULL);
209:   return(0);
210: }

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

218:   if (Z) {
219:     if (dm && dm != ts->dm) {
220:       DMGetNamedGlobalVector(dm,"TSEIMEX_Z",Z);
221:     } else *Z = ext->Z;
222:   }
223:   if (Ydot) {
224:     if (dm && dm != ts->dm) {
225:       DMGetNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
226:     } else *Ydot = ext->Ydot;
227:   }
228:   if (YdotI) {
229:     if (dm && dm != ts->dm) {
230:       DMGetNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
231:     } else *YdotI = ext->YdotI;
232:   }
233:   if (YdotRHS) {
234:     if (dm && dm != ts->dm) {
235:       DMGetNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
236:     } else *YdotRHS = ext->YdotRHS;
237:   }
238:   return(0);
239: }

241: static PetscErrorCode TSEIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI,Vec *YdotRHS)
242: {

246:   if (Z) {
247:     if (dm && dm != ts->dm) {
248:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_Z",Z);
249:     }
250:   }
251:   if (Ydot) {
252:     if (dm && dm != ts->dm) {
253:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
254:     }
255:   }
256:   if (YdotI) {
257:     if (dm && dm != ts->dm) {
258:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
259:     }
260:   }
261:   if (YdotRHS) {
262:     if (dm && dm != ts->dm) {
263:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
264:     }
265:   }
266:   return(0);
267: }

269: /*
270:   This defines the nonlinear equation that is to be solved with SNES
271:   Fn[t0+Theta*dt, U, (U-U0)*shift] = 0
272:   In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U))
273:   Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h
274: */
275: static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes,Vec X,Vec G,TS ts)
276: {
277:   TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
278:   PetscErrorCode  ierr;
279:   Vec             Ydot,Z;
280:   DM              dm,dmsave;

283:   VecZeroEntries(G);

285:   SNESGetDM(snes,&dm);
286:   TSEIMEXGetVecs(ts,dm,&Z,&Ydot,NULL,NULL);
287:   VecZeroEntries(Ydot);
288:   dmsave = ts->dm;
289:   ts->dm = dm;
290:   TSComputeIFunction(ts,ext->ctime,X,Ydot,G,PETSC_FALSE);
291:   /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function.  */
292:   VecCopy(G,Ydot);
293:   ts->dm = dmsave;
294:   TSEIMEXRestoreVecs(ts,dm,&Z,&Ydot,NULL,NULL);

296:   return(0);
297: }

299: /*
300:  This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y))
301:  */
302: static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)
303: {
304:   TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
305:   Vec             Ydot;
306:   PetscErrorCode  ierr;
307:   DM              dm,dmsave;
309:   SNESGetDM(snes,&dm);
310:   TSEIMEXGetVecs(ts,dm,NULL,&Ydot,NULL,NULL);
311:   /*  VecZeroEntries(Ydot); */
312:   /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */
313:   dmsave = ts->dm;
314:   ts->dm = dm;
315:   TSComputeIJacobian(ts,ts->ptime,X,Ydot,ext->shift,A,B,PETSC_TRUE);
316:   ts->dm = dmsave;
317:   TSEIMEXRestoreVecs(ts,dm,NULL,&Ydot,NULL,NULL);
318:   return(0);
319: }

321: static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine,DM coarse,void *ctx)
322: {

325:   return(0);
326: }

328: static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
329: {
330:   TS ts = (TS)ctx;
332:   Vec Z,Z_c;

335:   TSEIMEXGetVecs(ts,fine,&Z,NULL,NULL,NULL);
336:   TSEIMEXGetVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
337:   MatRestrict(restrct,Z,Z_c);
338:   VecPointwiseMult(Z_c,rscale,Z_c);
339:   TSEIMEXRestoreVecs(ts,fine,&Z,NULL,NULL,NULL);
340:   TSEIMEXRestoreVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
341:   return(0);
342: }

344: static PetscErrorCode TSSetUp_EIMEX(TS ts)
345: {
346:   TS_EIMEX       *ext = (TS_EIMEX*)ts->data;
348:   DM             dm;

351:   if (!ext->N){ /* ext->max_rows not set */
352:     TSEIMEXSetMaxRows(ts,TSEIMEXDefault);
353:   }
354:   if (-1 == ext->row_ind && -1 == ext->col_ind){
355:         TSEIMEXSetRowCol(ts,ext->max_rows,ext->max_rows);
356:   } else{/* ext->row_ind and col_ind already set */
357:     if (ext->ord_adapt){
358:       PetscInfo(ts,"Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");
359:     }
360:   }

362:   if (ext->ord_adapt){
363:     ext->nstages = 2; /* Start with the 2-stage scheme */
364:     TSEIMEXSetRowCol(ts,ext->nstages,ext->nstages);
365:   } else{
366:     ext->nstages = ext->max_rows; /* by default nstages is the same as max_rows, this can be changed by setting order adaptivity */
367:   }

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

371:   VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);/* full T table */
372:   VecDuplicate(ts->vec_sol,&ext->YdotI);
373:   VecDuplicate(ts->vec_sol,&ext->YdotRHS);
374:   VecDuplicate(ts->vec_sol,&ext->Ydot);
375:   VecDuplicate(ts->vec_sol,&ext->VecSolPrev);
376:   VecDuplicate(ts->vec_sol,&ext->Y);
377:   VecDuplicate(ts->vec_sol,&ext->Z);
378:   TSGetDM(ts,&dm);
379:   if (dm) {
380:     DMCoarsenHookAdd(dm,DMCoarsenHook_TSEIMEX,DMRestrictHook_TSEIMEX,ts);
381:   }
382:   return(0);
383: }

385: static PetscErrorCode TSSetFromOptions_EIMEX(PetscOptionItems *PetscOptionsObject,TS ts)
386: {
387:   TS_EIMEX       *ext = (TS_EIMEX*)ts->data;
389:   PetscInt       tindex[2];
390:   PetscInt       np = 2, nrows=TSEIMEXDefault;

393:   tindex[0] = TSEIMEXDefault;
394:   tindex[1] = TSEIMEXDefault;
395:   PetscOptionsHead(PetscOptionsObject,"EIMEX ODE solver options");
396:   {
397:     PetscBool flg;
398:     PetscOptionsInt("-ts_eimex_max_rows","Define the maximum number of rows used","TSEIMEXSetMaxRows",nrows,&nrows,&flg); /* default value 3 */
399:     if (flg){
400:       TSEIMEXSetMaxRows(ts,nrows);
401:     }
402:     PetscOptionsIntArray("-ts_eimex_row_col","Return the specific term in the T table","TSEIMEXSetRowCol",tindex,&np,&flg);
403:     if (flg){
404:       TSEIMEXSetRowCol(ts,tindex[0],tindex[1]);
405:     }
406:     PetscOptionsBool("-ts_eimex_order_adapt","Solve the problem with adaptive order","TSEIMEXSetOrdAdapt",ext->ord_adapt,&ext->ord_adapt,NULL);
407:   }
408:   PetscOptionsTail();
409:   return(0);
410: }

412: static PetscErrorCode TSView_EIMEX(TS ts,PetscViewer viewer)
413: {
415:   return(0);
416: }

418: /*@C
419:   TSEIMEXSetMaxRows - Set the maximum number of rows for EIMEX schemes

421:   Logically collective

423:   Input Parameter:
424: +  ts - timestepping context
425: -  nrows - maximum number of rows

427:   Level: intermediate

429: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
430: @*/
431: PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows)
432: {
436:   PetscTryMethod(ts,"TSEIMEXSetMaxRows_C",(TS,PetscInt),(ts,nrows));
437:   return(0);
438: }

440: /*@C
441:   TSEIMEXSetRowCol - Set the type index in the T table for the return value

443:   Logically collective

445:   Input Parameter:
446: +  ts - timestepping context
447: -  tindex - index in the T table

449:   Level: intermediate

451: .seealso: TSEIMEXSetMaxRows(), TSEIMEXSetOrdAdapt(), TSEIMEX
452: @*/
453: PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col)
454: {
458:   PetscTryMethod(ts,"TSEIMEXSetRowCol_C",(TS,PetscInt, PetscInt),(ts,row,col));
459:   return(0);
460: }

462: /*@C
463:   TSEIMEXSetOrdAdapt - Set the order adaptativity

465:   Logically collective

467:   Input Parameter:
468: +  ts - timestepping context
469: -  tindex - index in the T table

471:   Level: intermediate

473: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
474: @*/
475: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
476: {
480:   PetscTryMethod(ts,"TSEIMEXSetOrdAdapt_C",(TS,PetscBool),(ts,flg));
481:   return(0);
482: }

484: static PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts,PetscInt nrows)
485: {
486:   TS_EIMEX *ext = (TS_EIMEX*)ts->data;
488:   PetscInt       i;

491:   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);
492:   PetscFree(ext->N);
493:   ext->max_rows = nrows;
494:   PetscMalloc1(nrows,&ext->N);
495:   for (i=0;i<nrows;i++) ext->N[i]=i+1;
496:   return(0);
497: }

499: static PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts,PetscInt row,PetscInt col)
500: {
501:   TS_EIMEX *ext = (TS_EIMEX*)ts->data;

504:   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);
505:   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);
506:   if (col > row) SETERRQ2(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The column index (%d) exceeds the row index (%d)\n",col,row);

508:   ext->row_ind = row - 1;
509:   ext->col_ind = col - 1; /* Array index in C starts from 0 */
510:   return(0);
511: }

513: static PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts,PetscBool flg)
514: {
515:   TS_EIMEX *ext = (TS_EIMEX*)ts->data;
517:   ext->ord_adapt = flg;
518:   return(0);
519: }

521: /*MC
522:       TSEIMEX - Time stepping with Extrapolated IMEX methods.

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

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

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

533:   The general system is written as

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

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

541:   Another common form for the system is

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

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

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

549:  References
550:   E. Constantinescu and A. Sandu, Extrapolated implicit-explicit time stepping, SIAM Journal on Scientific
551: Computing, 31 (2010), pp. 4452-4477.

553:       Level: beginner

555: .seealso:  TSCreate(), TS, TSSetType(), TSEIMEXSetMaxRows(), TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt()

557:  M*/
558: PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS ts)
559: {
560:   TS_EIMEX       *ext;


565:   ts->ops->reset          = TSReset_EIMEX;
566:   ts->ops->destroy        = TSDestroy_EIMEX;
567:   ts->ops->view           = TSView_EIMEX;
568:   ts->ops->setup          = TSSetUp_EIMEX;
569:   ts->ops->step           = TSStep_EIMEX;
570:   ts->ops->interpolate    = TSInterpolate_EIMEX;
571:   ts->ops->evaluatestep   = TSEvaluateStep_EIMEX;
572:   ts->ops->setfromoptions = TSSetFromOptions_EIMEX;
573:   ts->ops->snesfunction   = SNESTSFormFunction_EIMEX;
574:   ts->ops->snesjacobian   = SNESTSFormJacobian_EIMEX;
575:   ts->default_adapt_type  = TSADAPTNONE;

577:   ts->usessnes = PETSC_TRUE;

579:   PetscNewLog(ts,&ext);
580:   ts->data = (void*)ext;

582:   ext->ord_adapt = PETSC_FALSE; /* By default, no order adapativity */
583:   ext->row_ind   = -1;
584:   ext->col_ind   = -1;
585:   ext->max_rows  = TSEIMEXDefault;
586:   ext->nstages   = TSEIMEXDefault;

588:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C", TSEIMEXSetMaxRows_EIMEX);
589:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C",  TSEIMEXSetRowCol_EIMEX);
590:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",TSEIMEXSetOrdAdapt_EIMEX);
591:   return(0);
592: }