Actual source code: eimex.c

petsc-3.9.2 2018-05-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 */
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;
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;
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){
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:         }
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);

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 */
368:       PetscInfo(ts,"Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");
369:     }
370:   }

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:   }

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) {
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;
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:     }
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

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

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

478:   Logically collective

480:   Input Parameter:
481: +  ts - timestepping context
482: -  tindex - index in the T table

484:   Level: intermediate

487: @*/
488: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
489: {
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;
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;

591:   ts->usessnes = PETSC_TRUE;

593:   PetscNewLog(ts,&ext);
594:   ts->data = (void*)ext;