Actual source code: eimex.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  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:           TSErrorNormWRMS(ts,T[Map(ext->nstages-1,ext->nstages-2,ext->nstages)],&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 */

192:   ts->ptime += ts->time_step;
193:   ts->steps++;
194:   ext->status = TS_STEP_COMPLETE;

196:   if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
197:   return(0);
198: }


201: /* cubic Hermit spline */
204: static PetscErrorCode TSInterpolate_EIMEX(TS ts,PetscReal itime,Vec X)
205: {
206:   TS_EIMEX       *ext = (TS_EIMEX*)ts->data;
207:   PetscReal      t,a,b;
208:   Vec            Y0=ext->VecSolPrev,Y1=ext->Y,
209:                          Ydot=ext->Ydot,YdotI=ext->YdotI;
210:   const PetscReal h = ts->time_step_prev;
213:   t = (itime -ts->ptime + h)/h;
214:   /* YdotI = -f(x)-g(x) */

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

219:   a    = 2.0*t*t*t - 3.0*t*t + 1.0;
220:   b    = -(t*t*t - 2.0*t*t + t)*h;
221:   VecAXPBYPCZ(X,a,b,0.0,Y0,YdotI);

223:   TSComputeIFunction(ts,ts->ptime,Y1,Ydot,YdotI,PETSC_FALSE);
224:   a    = -2.0*t*t*t+3.0*t*t;
225:   b    = -(t*t*t - t*t)*h;
226:   VecAXPBYPCZ(X,a,b,1.0,Y1,YdotI);

228:   return(0);
229: }


234: static PetscErrorCode TSReset_EIMEX(TS ts)
235: {
236:   TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
237:   PetscInt        ns;
238:   PetscErrorCode  ierr;

241:   ns = ext->nstages;
242:   VecDestroyVecs((1+ns)*ns/2,&ext->T);
243:   VecDestroy(&ext->Y);
244:   VecDestroy(&ext->Z);
245:   VecDestroy(&ext->YdotRHS);
246:   VecDestroy(&ext->YdotI);
247:   VecDestroy(&ext->Ydot);
248:   VecDestroy(&ext->VecSolPrev);
249:   PetscFree(ext->N);
250:   return(0);
251: }

255: static PetscErrorCode TSDestroy_EIMEX(TS ts)
256: {
257:   PetscErrorCode  ierr;

260:   TSReset_EIMEX(ts);
261:   PetscFree(ts->data);
262:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C",NULL);
263:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C",NULL);
264:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",NULL);

266:   return(0);
267: }


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

278:   if (Z) {
279:     if (dm && dm != ts->dm) {
280:       DMGetNamedGlobalVector(dm,"TSEIMEX_Z",Z);
281:     } else *Z = ext->Z;
282:   }
283:   if (Ydot) {
284:     if (dm && dm != ts->dm) {
285:       DMGetNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
286:     } else *Ydot = ext->Ydot;
287:   }
288:   if (YdotI) {
289:     if (dm && dm != ts->dm) {
290:       DMGetNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
291:     } else *YdotI = ext->YdotI;
292:   }
293:   if (YdotRHS) {
294:     if (dm && dm != ts->dm) {
295:       DMGetNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
296:     } else *YdotRHS = ext->YdotRHS;
297:   }
298:   return(0);
299: }


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

309:   if (Z) {
310:     if (dm && dm != ts->dm) {
311:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_Z",Z);
312:     }
313:   }
314:   if (Ydot) {
315:     if (dm && dm != ts->dm) {
316:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
317:     }
318:   }
319:   if (YdotI) {
320:     if (dm && dm != ts->dm) {
321:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
322:     }
323:   }
324:   if (YdotRHS) {
325:     if (dm && dm != ts->dm) {
326:       DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
327:     }
328:   }
329:   return(0);
330: }


333: /*
334:   This defines the nonlinear equation that is to be solved with SNES
335:   Fn[t0+Theta*dt, U, (U-U0)*shift] = 0
336:   In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U))
337:   Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h
338: */
341: static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes,Vec X,Vec G,TS ts)
342: {
343:   TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
344:   PetscErrorCode  ierr;
345:   Vec             Ydot,Z;
346:   DM              dm,dmsave;

349:   VecZeroEntries(G);

351:   SNESGetDM(snes,&dm);
352:   TSEIMEXGetVecs(ts,dm,&Z,&Ydot,NULL,NULL);
353:   VecZeroEntries(Ydot);
354:   dmsave = ts->dm;
355:   ts->dm = dm;
356:   TSComputeIFunction(ts,ext->ctime,X,Ydot,G,PETSC_FALSE);
357:   /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function.  */
358:   VecCopy(G,Ydot);
359:   ts->dm = dmsave;
360:   TSEIMEXRestoreVecs(ts,dm,&Z,&Ydot,NULL,NULL);

362:   return(0);
363: }

365: /*
366:  This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y))
367:  */
370: static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)
371: {
372:   TS_EIMEX        *ext = (TS_EIMEX*)ts->data;
373:   Vec             Ydot;
374:   PetscErrorCode  ierr;
375:   DM              dm,dmsave;
377:   SNESGetDM(snes,&dm);
378:   TSEIMEXGetVecs(ts,dm,NULL,&Ydot,NULL,NULL);
379:   /*  VecZeroEntries(Ydot); */
380:   /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */
381:   dmsave = ts->dm;
382:   ts->dm = dm;
383:   TSComputeIJacobian(ts,ts->ptime,X,Ydot,ext->shift,A,B,PETSC_TRUE);
384:   ts->dm = dmsave;
385:   TSEIMEXRestoreVecs(ts,dm,NULL,&Ydot,NULL,NULL);
386:   return(0);
387: }

391: static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine,DM coarse,void *ctx)
392: {

395:   return(0);
396: }

400: static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
401: {
402:   TS ts = (TS)ctx;
404:   Vec Z,Z_c;

407:   TSEIMEXGetVecs(ts,fine,&Z,NULL,NULL,NULL);
408:   TSEIMEXGetVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
409:   MatRestrict(restrct,Z,Z_c);
410:   VecPointwiseMult(Z_c,rscale,Z_c);
411:   TSEIMEXRestoreVecs(ts,fine,&Z,NULL,NULL,NULL);
412:   TSEIMEXRestoreVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
413:   return(0);
414: }


419: static PetscErrorCode TSSetUp_EIMEX(TS ts)
420: {
421:   TS_EIMEX       *ext = (TS_EIMEX*)ts->data;
423:   DM             dm;

426:   if (!ext->N){ /* ext->max_rows not set */
427:     TSEIMEXSetMaxRows(ts,TSEIMEXDefault);
428:   }
429:   if(-1 == ext->row_ind && -1 == ext->col_ind){
430:         TSEIMEXSetRowCol(ts,ext->max_rows,ext->max_rows);
431:   }else{/* ext->row_ind and col_ind already set */
432:     if(ext->ord_adapt){
433:       PetscInfo(ts,"Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");
434:         }
435:   }

437:   if(ext->ord_adapt){
438:     ext->nstages = 2; /* Start with the 2-stage scheme */
439:         TSEIMEXSetRowCol(ts,ext->nstages,ext->nstages);
440:   }else{
441:         ext->nstages = ext->max_rows; /* by default nstages is the same as max_rows, this can be changed by setting order adaptivity */
442:   }

444:   VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);/* full T table */
445:   VecDuplicate(ts->vec_sol,&ext->YdotI);
446:   VecDuplicate(ts->vec_sol,&ext->YdotRHS);
447:   VecDuplicate(ts->vec_sol,&ext->Ydot);
448:   VecDuplicate(ts->vec_sol,&ext->VecSolPrev);
449:   VecDuplicate(ts->vec_sol,&ext->Y);
450:   VecDuplicate(ts->vec_sol,&ext->Z);
451:   TSGetDM(ts,&dm);
452:   if (dm) {
453:     DMCoarsenHookAdd(dm,DMCoarsenHook_TSEIMEX,DMRestrictHook_TSEIMEX,ts);
454:   }
455:   return(0);
456: }

460: static PetscErrorCode TSSetFromOptions_EIMEX(TS ts)
461: {
462:   TS_EIMEX       *ext = (TS_EIMEX*)ts->data;
464:   PetscInt       tindex[2];
465:   PetscInt       np = 2, nrows=TSEIMEXDefault;

468:   tindex[0] = TSEIMEXDefault;
469:   tindex[1] = TSEIMEXDefault;
470:   PetscOptionsHead("EIMEX ODE solver options");
471:   {
472:     PetscBool flg;
473:     flg  = PETSC_FALSE;
474:     PetscOptionsInt("-ts_eimex_max_rows","Define the maximum number of rows used","TSEIMEXSetMaxRows",nrows,&nrows,&flg); /* default value 3 */
475:     if(flg){
476:       TSEIMEXSetMaxRows(ts,nrows);
477:     }
478:     PetscOptionsIntArray("-ts_eimex_row_col","Return the specific term in the T table","TSEIMEXSetRowCol",tindex,&np,&flg);
479:     if(flg){
480:       TSEIMEXSetRowCol(ts,tindex[0],tindex[1]);
481:     }
482:     PetscOptionsBool("-ts_eimex_order_adapt","Solve the problem with adaptive order","TSEIMEXSetOrdAdapt",ext->ord_adapt,&ext->ord_adapt,NULL);
483:     SNESSetFromOptions(ts->snes);
484:   }
485:   PetscOptionsTail();
486:   return(0);
487: }

491: static PetscErrorCode TSView_EIMEX(TS ts,PetscViewer viewer)
492: {
493:   /*  TS_EIMEX         *ext = (TS_EIMEX*)ts->data; */
494:   PetscBool        iascii;
495:   PetscErrorCode   ierr;

498:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
499:   if (iascii) {
500:     PetscViewerASCIIPrintf(viewer,"  EIMEX\n");
501:   }
502:   SNESView(ts->snes,viewer);
503:   return(0);
504: }


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

512:   Logically collective

514:   Input Parameter:
515: +  ts - timestepping context
516: -  nrows - maximum number of rows

518:   Level: intermediate

520: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
521: @*/
522: PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows)
523: {
527:   PetscTryMethod(ts,"TSEIMEXSetMaxRows_C",(TS,PetscInt),(ts,nrows));
528:   return(0);
529: }


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

537:   Logically collective

539:   Input Parameter:
540: +  ts - timestepping context
541: -  tindex - index in the T table

543:   Level: intermediate

545: .seealso: TSEIMEXSetMaxRows(), TSEIMEXSetOrdAdapt(), TSEIMEX
546: @*/
547: PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col)
548: {
552:   PetscTryMethod(ts,"TSEIMEXSetRowCol_C",(TS,PetscInt, PetscInt),(ts,row,col));
553:   return(0);
554: }


559: /*@C
560:   TSEIMEXSetOrdAdapt - Set the order adaptativity

562:   Logically collective

564:   Input Parameter:
565: +  ts - timestepping context
566: -  tindex - index in the T table

568:   Level: intermediate

570: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
571: @*/
572: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
573: {
577:   PetscTryMethod(ts,"TSEIMEXSetOrdAdapt_C",(TS,PetscBool),(ts,flg));
578:   return(0);
579: }


584: static PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts,PetscInt nrows)
585: {
586:   TS_EIMEX *ext = (TS_EIMEX*)ts->data;
588:   PetscInt       i;

591:   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);
592:   PetscFree(ext->N);
593:   ext->max_rows = nrows;
594:   PetscMalloc1(nrows,&ext->N);
595:   for(i=0;i<nrows;i++) ext->N[i]=i+1;
596:   return(0);
597: }

601: static PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts,PetscInt row,PetscInt col)
602: {
603:   TS_EIMEX *ext = (TS_EIMEX*)ts->data;

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

610:   ext->row_ind = row - 1;
611:   ext->col_ind = col - 1; /* Array index in C starts from 0 */
612:   return(0);
613: }

617: static PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts,PetscBool flg)
618: {
619:   TS_EIMEX *ext = (TS_EIMEX*)ts->data;
621:   ext->ord_adapt = flg;
622:   return(0);
623: }

625: /* ------------------------------------------------------------ */
626: /*MC
627:       TSEIMEX - ODE solver using extrapolated IMEX schemes
628:   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().

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

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

635:   Level: beginner

637: .seealso:  TSCreate(), TS
638: M*/
641: PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS ts)
642: {
643:   TS_EIMEX       *ext;


648:   ts->ops->reset          = TSReset_EIMEX;
649:   ts->ops->destroy        = TSDestroy_EIMEX;
650:   ts->ops->view           = TSView_EIMEX;
651:   ts->ops->setup          = TSSetUp_EIMEX;
652:   ts->ops->step           = TSStep_EIMEX;
653:   ts->ops->interpolate    = TSInterpolate_EIMEX;
654:   ts->ops->evaluatestep   = TSEvaluateStep_EIMEX;
655:   ts->ops->setfromoptions = TSSetFromOptions_EIMEX;
656:   ts->ops->snesfunction   = SNESTSFormFunction_EIMEX;
657:   ts->ops->snesjacobian   = SNESTSFormJacobian_EIMEX;

659:   PetscNewLog(ts,&ext);
660:   ts->data = (void*)ext;

662:   ext->ord_adapt = PETSC_FALSE; /* By default, no order adapativity */
663:   ext->row_ind   = -1;
664:   ext->col_ind   = -1;
665:   ext->max_rows  = TSEIMEXDefault;
666:   ext->nstages   = TSEIMEXDefault;

668:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C", TSEIMEXSetMaxRows_EIMEX);
669:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C",  TSEIMEXSetRowCol_EIMEX);
670:   PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",TSEIMEXSetOrdAdapt_EIMEX);
671:   return(0);
672: }