Actual source code: eimex.c
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;
37: VecCopy(ext->T[Map(ext->row_ind,ext->col_ind,ns)],X);
38: return 0;
39: }
41: static PetscErrorCode TSStage_EIMEX(TS ts,PetscInt istage)
42: {
43: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
44: PetscReal h;
45: Vec Y=ext->Y, Z=ext->Z;
46: SNES snes;
47: TSAdapt adapt;
48: PetscInt i,its,lits;
49: PetscBool accept;
51: TSGetSNES(ts,&snes);
52: h = ts->time_step/ext->N[istage];/* step size for the istage-th stage */
53: ext->shift = 1./h;
54: SNESSetLagJacobian(snes,-2); /* Recompute the Jacobian on this solve, but not again */
55: VecCopy(ext->VecSolPrev,Y); /* Take the previous solution as initial step */
57: for (i=0; i<ext->N[istage]; i++) {
58: ext->ctime = ts->ptime + h*i;
59: VecCopy(Y,Z);/* Save the solution of the previous substep */
60: SNESSolve(snes,NULL,Y);
61: SNESGetIterationNumber(snes,&its);
62: SNESGetLinearSolveIterations(snes,&lits);
63: ts->snes_its += its; ts->ksp_its += lits;
64: TSGetAdapt(ts,&adapt);
65: TSAdaptCheckStage(adapt,ts,ext->ctime,Y,&accept);
66: }
67: return 0;
68: }
70: static PetscErrorCode TSStep_EIMEX(TS ts)
71: {
72: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
73: const PetscInt ns = ext->nstages;
74: Vec *T=ext->T, Y=ext->Y;
76: SNES snes;
77: PetscInt i,j;
78: PetscBool accept = PETSC_FALSE;
79: PetscErrorCode ierr;
80: PetscReal alpha,local_error,local_error_a,local_error_r;
82: TSGetSNES(ts,&snes);
83: SNESSetType(snes,"ksponly");
84: ext->status = TS_STEP_INCOMPLETE;
86: VecCopy(ts->vec_sol,ext->VecSolPrev);
88: /* Apply n_j steps of the base method to obtain solutions of T(j,1),1<=j<=s */
89: for (j=0; j<ns; j++) {
90: TSStage_EIMEX(ts,j);
91: VecCopy(Y,T[j]);
92: }
94: for (i=1;i<ns;i++) {
95: for (j=i;j<ns;j++) {
96: alpha = -(PetscReal)ext->N[j]/ext->N[j-i];
97: 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] */ierr;
98: alpha = 1.0/(1.0 + alpha);
99: VecScale(T[Map(j,i,ns)],alpha);
100: }
101: }
103: TSEvaluateStep(ts,ns,ts->vec_sol,NULL);/*update ts solution */
105: if (ext->ord_adapt && ext->nstages < ext->max_rows) {
106: accept = PETSC_FALSE;
107: while (!accept && ext->nstages < ext->max_rows) {
108: 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);
109: accept = (local_error < 1.0)? PETSC_TRUE : PETSC_FALSE;
111: if (!accept) {/* add one more stage*/
112: TSStage_EIMEX(ts,ext->nstages);
113: ext->nstages++; ext->row_ind++; ext->col_ind++;
114: /*T table need to be recycled*/
115: VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);
116: for (i=0; i<ext->nstages-1; i++) {
117: for (j=0; j<=i; j++) {
118: VecCopy(T[Map(i,j,ext->nstages-1)],ext->T[Map(i,j,ext->nstages)]);
119: }
120: }
121: VecDestroyVecs(ext->nstages*(ext->nstages-1)/2,&T);
122: T = ext->T; /*reset the pointer*/
123: /*recycling finished, store the new solution*/
124: VecCopy(Y,T[ext->nstages-1]);
125: /*extrapolation for the newly added stage*/
126: for (i=1;i<ext->nstages;i++) {
127: alpha = -(PetscReal)ext->N[ext->nstages-1]/ext->N[ext->nstages-1-i];
128: 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]*/ierr;
129: alpha = 1.0/(1.0 + alpha);
130: VecScale(T[Map(ext->nstages-1,i,ext->nstages)],alpha);
131: }
132: /*update ts solution */
133: TSEvaluateStep(ts,ext->nstages,ts->vec_sol,NULL);
134: }/*end if !accept*/
135: }/*end while*/
137: if (ext->nstages == ext->max_rows) {
138: PetscInfo(ts,"Max number of rows has been used\n");
139: }
140: }/*end if ext->ord_adapt*/
141: ts->ptime += ts->time_step;
142: ext->status = TS_STEP_COMPLETE;
144: if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
145: return 0;
146: }
148: /* cubic Hermit spline */
149: static PetscErrorCode TSInterpolate_EIMEX(TS ts,PetscReal itime,Vec X)
150: {
151: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
152: PetscReal t,a,b;
153: Vec Y0=ext->VecSolPrev,Y1=ext->Y,Ydot=ext->Ydot,YdotI=ext->YdotI;
154: const PetscReal h = ts->ptime - ts->ptime_prev;
155: t = (itime -ts->ptime + h)/h;
156: /* YdotI = -f(x)-g(x) */
158: VecZeroEntries(Ydot);
159: TSComputeIFunction(ts,ts->ptime-h,Y0,Ydot,YdotI,PETSC_FALSE);
161: a = 2.0*t*t*t - 3.0*t*t + 1.0;
162: b = -(t*t*t - 2.0*t*t + t)*h;
163: VecAXPBYPCZ(X,a,b,0.0,Y0,YdotI);
165: TSComputeIFunction(ts,ts->ptime,Y1,Ydot,YdotI,PETSC_FALSE);
166: a = -2.0*t*t*t+3.0*t*t;
167: b = -(t*t*t - t*t)*h;
168: VecAXPBYPCZ(X,a,b,1.0,Y1,YdotI);
170: return 0;
171: }
173: static PetscErrorCode TSReset_EIMEX(TS ts)
174: {
175: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
176: PetscInt ns;
178: ns = ext->nstages;
179: VecDestroyVecs((1+ns)*ns/2,&ext->T);
180: VecDestroy(&ext->Y);
181: VecDestroy(&ext->Z);
182: VecDestroy(&ext->YdotRHS);
183: VecDestroy(&ext->YdotI);
184: VecDestroy(&ext->Ydot);
185: VecDestroy(&ext->VecSolPrev);
186: PetscFree(ext->N);
187: return 0;
188: }
190: static PetscErrorCode TSDestroy_EIMEX(TS ts)
191: {
192: TSReset_EIMEX(ts);
193: PetscFree(ts->data);
194: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C",NULL);
195: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C",NULL);
196: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",NULL);
197: return 0;
198: }
200: static PetscErrorCode TSEIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI, Vec *YdotRHS)
201: {
202: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
204: if (Z) {
205: if (dm && dm != ts->dm) {
206: DMGetNamedGlobalVector(dm,"TSEIMEX_Z",Z);
207: } else *Z = ext->Z;
208: }
209: if (Ydot) {
210: if (dm && dm != ts->dm) {
211: DMGetNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
212: } else *Ydot = ext->Ydot;
213: }
214: if (YdotI) {
215: if (dm && dm != ts->dm) {
216: DMGetNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
217: } else *YdotI = ext->YdotI;
218: }
219: if (YdotRHS) {
220: if (dm && dm != ts->dm) {
221: DMGetNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
222: } else *YdotRHS = ext->YdotRHS;
223: }
224: return 0;
225: }
227: static PetscErrorCode TSEIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI,Vec *YdotRHS)
228: {
229: if (Z) {
230: if (dm && dm != ts->dm) {
231: DMRestoreNamedGlobalVector(dm,"TSEIMEX_Z",Z);
232: }
233: }
234: if (Ydot) {
235: if (dm && dm != ts->dm) {
236: DMRestoreNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);
237: }
238: }
239: if (YdotI) {
240: if (dm && dm != ts->dm) {
241: DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);
242: }
243: }
244: if (YdotRHS) {
245: if (dm && dm != ts->dm) {
246: DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);
247: }
248: }
249: return 0;
250: }
252: /*
253: This defines the nonlinear equation that is to be solved with SNES
254: Fn[t0+Theta*dt, U, (U-U0)*shift] = 0
255: In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U))
256: Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h
257: */
258: static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes,Vec X,Vec G,TS ts)
259: {
260: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
261: Vec Ydot,Z;
262: DM dm,dmsave;
264: VecZeroEntries(G);
266: SNESGetDM(snes,&dm);
267: TSEIMEXGetVecs(ts,dm,&Z,&Ydot,NULL,NULL);
268: VecZeroEntries(Ydot);
269: dmsave = ts->dm;
270: ts->dm = dm;
271: TSComputeIFunction(ts,ext->ctime,X,Ydot,G,PETSC_FALSE);
272: /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function. */
273: VecCopy(G,Ydot);
274: ts->dm = dmsave;
275: TSEIMEXRestoreVecs(ts,dm,&Z,&Ydot,NULL,NULL);
277: return 0;
278: }
280: /*
281: This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y))
282: */
283: static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)
284: {
285: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
286: Vec Ydot;
287: DM dm,dmsave;
288: SNESGetDM(snes,&dm);
289: TSEIMEXGetVecs(ts,dm,NULL,&Ydot,NULL,NULL);
290: /* VecZeroEntries(Ydot); */
291: /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */
292: dmsave = ts->dm;
293: ts->dm = dm;
294: TSComputeIJacobian(ts,ts->ptime,X,Ydot,ext->shift,A,B,PETSC_TRUE);
295: ts->dm = dmsave;
296: TSEIMEXRestoreVecs(ts,dm,NULL,&Ydot,NULL,NULL);
297: return 0;
298: }
300: static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine,DM coarse,void *ctx)
301: {
302: return 0;
303: }
305: static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
306: {
307: TS ts = (TS)ctx;
308: Vec Z,Z_c;
310: TSEIMEXGetVecs(ts,fine,&Z,NULL,NULL,NULL);
311: TSEIMEXGetVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
312: MatRestrict(restrct,Z,Z_c);
313: VecPointwiseMult(Z_c,rscale,Z_c);
314: TSEIMEXRestoreVecs(ts,fine,&Z,NULL,NULL,NULL);
315: TSEIMEXRestoreVecs(ts,coarse,&Z_c,NULL,NULL,NULL);
316: return 0;
317: }
319: static PetscErrorCode TSSetUp_EIMEX(TS ts)
320: {
321: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
322: DM dm;
324: if (!ext->N) { /* ext->max_rows not set */
325: TSEIMEXSetMaxRows(ts,TSEIMEXDefault);
326: }
327: if (-1 == ext->row_ind && -1 == ext->col_ind) {
328: TSEIMEXSetRowCol(ts,ext->max_rows,ext->max_rows);
329: } else{/* ext->row_ind and col_ind already set */
330: if (ext->ord_adapt) {
331: PetscInfo(ts,"Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");
332: }
333: }
335: if (ext->ord_adapt) {
336: ext->nstages = 2; /* Start with the 2-stage scheme */
337: TSEIMEXSetRowCol(ts,ext->nstages,ext->nstages);
338: } else{
339: ext->nstages = ext->max_rows; /* by default nstages is the same as max_rows, this can be changed by setting order adaptivity */
340: }
342: TSGetAdapt(ts,&ts->adapt);
344: VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);/* full T table */
345: VecDuplicate(ts->vec_sol,&ext->YdotI);
346: VecDuplicate(ts->vec_sol,&ext->YdotRHS);
347: VecDuplicate(ts->vec_sol,&ext->Ydot);
348: VecDuplicate(ts->vec_sol,&ext->VecSolPrev);
349: VecDuplicate(ts->vec_sol,&ext->Y);
350: VecDuplicate(ts->vec_sol,&ext->Z);
351: TSGetDM(ts,&dm);
352: if (dm) {
353: DMCoarsenHookAdd(dm,DMCoarsenHook_TSEIMEX,DMRestrictHook_TSEIMEX,ts);
354: }
355: return 0;
356: }
358: static PetscErrorCode TSSetFromOptions_EIMEX(PetscOptionItems *PetscOptionsObject,TS ts)
359: {
360: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
361: PetscInt tindex[2];
362: PetscInt np = 2, nrows=TSEIMEXDefault;
364: tindex[0] = TSEIMEXDefault;
365: tindex[1] = TSEIMEXDefault;
366: PetscOptionsHead(PetscOptionsObject,"EIMEX ODE solver options");
367: {
368: PetscBool flg;
369: PetscOptionsInt("-ts_eimex_max_rows","Define the maximum number of rows used","TSEIMEXSetMaxRows",nrows,&nrows,&flg); /* default value 3 */
370: if (flg) {
371: TSEIMEXSetMaxRows(ts,nrows);
372: }
373: PetscOptionsIntArray("-ts_eimex_row_col","Return the specific term in the T table","TSEIMEXSetRowCol",tindex,&np,&flg);
374: if (flg) {
375: TSEIMEXSetRowCol(ts,tindex[0],tindex[1]);
376: }
377: PetscOptionsBool("-ts_eimex_order_adapt","Solve the problem with adaptive order","TSEIMEXSetOrdAdapt",ext->ord_adapt,&ext->ord_adapt,NULL);
378: }
379: PetscOptionsTail();
380: return 0;
381: }
383: static PetscErrorCode TSView_EIMEX(TS ts,PetscViewer viewer)
384: {
385: return 0;
386: }
388: /*@C
389: TSEIMEXSetMaxRows - Set the maximum number of rows for EIMEX schemes
391: Logically collective
393: Input Parameters:
394: + ts - timestepping context
395: - nrows - maximum number of rows
397: Level: intermediate
399: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
400: @*/
401: PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows)
402: {
404: PetscTryMethod(ts,"TSEIMEXSetMaxRows_C",(TS,PetscInt),(ts,nrows));
405: return 0;
406: }
408: /*@C
409: TSEIMEXSetRowCol - Set the type index in the T table for the return value
411: Logically collective
413: Input Parameters:
414: + ts - timestepping context
415: - tindex - index in the T table
417: Level: intermediate
419: .seealso: TSEIMEXSetMaxRows(), TSEIMEXSetOrdAdapt(), TSEIMEX
420: @*/
421: PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col)
422: {
424: PetscTryMethod(ts,"TSEIMEXSetRowCol_C",(TS,PetscInt, PetscInt),(ts,row,col));
425: return 0;
426: }
428: /*@C
429: TSEIMEXSetOrdAdapt - Set the order adaptativity
431: Logically collective
433: Input Parameters:
434: + ts - timestepping context
435: - tindex - index in the T table
437: Level: intermediate
439: .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX
440: @*/
441: PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg)
442: {
444: PetscTryMethod(ts,"TSEIMEXSetOrdAdapt_C",(TS,PetscBool),(ts,flg));
445: return 0;
446: }
448: static PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts,PetscInt nrows)
449: {
450: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
451: PetscInt i;
454: PetscFree(ext->N);
455: ext->max_rows = nrows;
456: PetscMalloc1(nrows,&ext->N);
457: for (i=0;i<nrows;i++) ext->N[i]=i+1;
458: return 0;
459: }
461: static PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts,PetscInt row,PetscInt col)
462: {
463: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
469: ext->row_ind = row - 1;
470: ext->col_ind = col - 1; /* Array index in C starts from 0 */
471: return 0;
472: }
474: static PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts,PetscBool flg)
475: {
476: TS_EIMEX *ext = (TS_EIMEX*)ts->data;
477: ext->ord_adapt = flg;
478: return 0;
479: }
481: /*MC
482: TSEIMEX - Time stepping with Extrapolated IMEX methods.
484: These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly nonlinear such that it
485: is expensive to solve with a fully implicit method. The user should provide the stiff part of the equation using TSSetIFunction() and the
486: non-stiff part with TSSetRHSFunction().
488: Notes:
489: The default is a 3-stage scheme, it can be changed with TSEIMEXSetMaxRows() or -ts_eimex_max_rows
491: This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).
493: The general system is written as
495: G(t,X,Xdot) = F(t,X)
497: where G represents the stiff part and F represents the non-stiff part. The user should provide the stiff part
498: of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction().
499: This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian.
501: Another common form for the system is
503: y'=f(x)+g(x)
505: The relationship between F,G and f,g is
507: G = y'-g(x), F = f(x)
509: References
510: E. Constantinescu and A. Sandu, Extrapolated implicit-explicit time stepping, SIAM Journal on Scientific
511: Computing, 31 (2010), pp. 4452-4477.
513: Level: beginner
515: .seealso: TSCreate(), TS, TSSetType(), TSEIMEXSetMaxRows(), TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt()
517: M*/
518: PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS ts)
519: {
520: TS_EIMEX *ext;
523: ts->ops->reset = TSReset_EIMEX;
524: ts->ops->destroy = TSDestroy_EIMEX;
525: ts->ops->view = TSView_EIMEX;
526: ts->ops->setup = TSSetUp_EIMEX;
527: ts->ops->step = TSStep_EIMEX;
528: ts->ops->interpolate = TSInterpolate_EIMEX;
529: ts->ops->evaluatestep = TSEvaluateStep_EIMEX;
530: ts->ops->setfromoptions = TSSetFromOptions_EIMEX;
531: ts->ops->snesfunction = SNESTSFormFunction_EIMEX;
532: ts->ops->snesjacobian = SNESTSFormJacobian_EIMEX;
533: ts->default_adapt_type = TSADAPTNONE;
535: ts->usessnes = PETSC_TRUE;
537: PetscNewLog(ts,&ext);
538: ts->data = (void*)ext;
540: ext->ord_adapt = PETSC_FALSE; /* By default, no order adapativity */
541: ext->row_ind = -1;
542: ext->col_ind = -1;
543: ext->max_rows = TSEIMEXDefault;
544: ext->nstages = TSEIMEXDefault;
546: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetMaxRows_C", TSEIMEXSetMaxRows_EIMEX);
547: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetRowCol_C", TSEIMEXSetRowCol_EIMEX);
548: PetscObjectComposeFunction((PetscObject)ts,"TSEIMEXSetOrdAdapt_C",TSEIMEXSetOrdAdapt_EIMEX);
549: return 0;
550: }