Actual source code: bdf.c
petsc-3.11.4 2019-09-28
1: /*
2: Code for timestepping with BDF methods
3: */
4: #include <petsc/private/tsimpl.h>
6: static PetscBool cited = PETSC_FALSE;
7: static const char citation[] =
8: "@book{Brenan1995,\n"
9: " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n"
10: " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n"
11: " publisher = {Society for Industrial and Applied Mathematics},\n"
12: " year = {1995},\n"
13: " doi = {10.1137/1.9781611971224},\n}\n";
15: typedef struct {
16: PetscInt k,n;
17: PetscReal time[6+2];
18: Vec work[6+2];
19: PetscReal shift;
20: Vec vec_dot;
21: Vec vec_lte;
23: PetscInt order;
24: TSStepStatus status;
25: } TS_BDF;
28: PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
29: {
30: PetscInt k,j;
31: for (k=0; k<n; k++)
32: for (L[k]=1, j=0; j<n; j++)
33: if (j != k)
34: L[k] *= (t - T[j])/(T[k] - T[j]);
35: }
37: PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])
38: {
39: PetscInt k,j,i;
40: for (k=0; k<n; k++)
41: for (dL[k]=0, j=0; j<n; j++)
42: if (j != k) {
43: PetscReal L = 1/(T[k] - T[j]);
44: for (i=0; i<n; i++)
45: if (i != j && i != k)
46: L *= (t - T[i])/(T[k] - T[i]);
47: dL[k] += L;
48: }
49: }
51: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
52: {
53: TS_BDF *bdf = (TS_BDF*)ts->data;
54: PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
55: Vec tail = bdf->work[n-1];
59: for (i=n-1; i>=2; i--) {
60: bdf->time[i] = bdf->time[i-1];
61: bdf->work[i] = bdf->work[i-1];
62: }
63: bdf->n = PetscMin(bdf->n+1,n-1);
64: bdf->time[1] = t;
65: bdf->work[1] = tail;
66: VecCopy(X,tail);
67: return(0);
68: }
70: static PetscErrorCode TSBDF_VecDot(TS ts,PetscInt order,PetscReal t,Vec X,Vec Xdot,PetscReal *shift)
71: {
72: TS_BDF *bdf = (TS_BDF*)ts->data;
73: PetscInt i,n = order+1;
74: PetscReal time[7];
75: Vec vecs[7];
76: PetscScalar alpha[7];
80: n = PetscMax(n,2);
81: time[0] = t; for (i=1; i<n; i++) time[i] = bdf->time[i];
82: vecs[0] = X; for (i=1; i<n; i++) vecs[i] = bdf->work[i];
83: LagrangeBasisDers(n,t,time,alpha);
84: VecZeroEntries(Xdot);
85: VecMAXPY(Xdot,n,alpha,vecs);
86: if (shift) *shift = PetscRealPart(alpha[0]);
87: return(0);
88: }
90: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
91: {
92: TS_BDF *bdf = (TS_BDF*)ts->data;
93: PetscInt i,n = order+1;
94: PetscReal *time = bdf->time;
95: Vec *vecs = bdf->work;
96: PetscScalar a[8],b[8],alpha[8];
100: LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
101: LagrangeBasisDers(n+1,time[0],time,b);
102: for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
103: VecZeroEntries(lte);
104: VecMAXPY(lte,n+1,alpha,vecs);
105: return(0);
106: }
108: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
109: {
110: TS_BDF *bdf = (TS_BDF*)ts->data;
111: PetscInt n = order+1;
112: PetscReal *time = bdf->time+1;
113: Vec *vecs = bdf->work+1;
114: PetscScalar alpha[7];
118: n = PetscMin(n,bdf->n);
119: LagrangeBasisVals(n,t,time,alpha);
120: VecZeroEntries(X);
121: VecMAXPY(X,n,alpha,vecs);
122: return(0);
123: }
125: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
126: {
127: TS_BDF *bdf = (TS_BDF*)ts->data;
128: PetscInt n = order+1;
129: PetscReal *time = bdf->time;
130: Vec *vecs = bdf->work;
131: PetscScalar alpha[7];
135: LagrangeBasisVals(n,t,time,alpha);
136: VecZeroEntries(X);
137: VecMAXPY(X,n,alpha,vecs);
138: return(0);
139: }
141: static PetscErrorCode TSBDF_SNESSolve(TS ts,Vec b,Vec x)
142: {
143: PetscInt nits,lits;
147: SNESSolve(ts->snes,b,x);
148: SNESGetIterationNumber(ts->snes,&nits);
149: SNESGetLinearSolveIterations(ts->snes,&lits);
150: ts->snes_its += nits; ts->ksp_its += lits;
151: return(0);
152: }
154: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
155: {
156: TS_BDF *bdf = (TS_BDF*)ts->data;
160: bdf->k = 1; bdf->n = 0;
161: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
163: bdf->time[0] = ts->ptime + ts->time_step/2;
164: VecCopy(bdf->work[1],bdf->work[0]);
165: TSPreStage(ts,bdf->time[0]);
166: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
167: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
168: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
169: if (!*accept) return(0);
171: bdf->k = PetscMin(2,bdf->order); bdf->n++;
172: VecCopy(bdf->work[0],bdf->work[2]);
173: bdf->time[2] = bdf->time[0];
174: return(0);
175: }
177: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
179: static PetscErrorCode TSStep_BDF(TS ts)
180: {
181: TS_BDF *bdf = (TS_BDF*)ts->data;
182: PetscInt rejections = 0;
183: PetscBool stageok,accept = PETSC_TRUE;
184: PetscReal next_time_step = ts->time_step;
188: PetscCitationsRegister(citation,&cited);
190: if (!ts->steprollback && !ts->steprestart) {
191: bdf->k = PetscMin(bdf->k+1,bdf->order);
192: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
193: }
195: bdf->status = TS_STEP_INCOMPLETE;
196: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
198: if (ts->steprestart) {
199: TSBDF_Restart(ts,&stageok);
200: if (!stageok) goto reject_step;
201: }
203: bdf->time[0] = ts->ptime + ts->time_step;
204: TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
205: TSPreStage(ts,bdf->time[0]);
206: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
207: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
208: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
209: if (!stageok) goto reject_step;
211: bdf->status = TS_STEP_PENDING;
212: TSAdaptCandidatesClear(ts->adapt);
213: TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
214: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
215: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
216: if (!accept) { ts->time_step = next_time_step; goto reject_step; }
218: VecCopy(bdf->work[0],ts->vec_sol);
219: ts->ptime += ts->time_step;
220: ts->time_step = next_time_step;
221: break;
223: reject_step:
224: ts->reject++; accept = PETSC_FALSE;
225: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
226: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
227: ts->reason = TS_DIVERGED_STEP_REJECTED;
228: }
229: }
230: return(0);
231: }
233: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
234: {
235: TS_BDF *bdf = (TS_BDF*)ts->data;
239: TSBDF_Interpolate(ts,bdf->k,t,X);
240: return(0);
241: }
243: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
244: {
245: TS_BDF *bdf = (TS_BDF*)ts->data;
246: PetscInt k = bdf->k;
247: PetscReal wltea,wlter;
248: Vec X = bdf->work[0], Y = bdf->vec_lte;
252: k = PetscMin(k,bdf->n-1);
253: TSBDF_VecLTE(ts,k,Y);
254: VecAXPY(Y,1,X);
255: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
256: if (order) *order = k + 1;
257: return(0);
258: }
260: static PetscErrorCode TSRollBack_BDF(TS ts)
261: {
262: TS_BDF *bdf = (TS_BDF*)ts->data;
266: VecCopy(bdf->work[1],ts->vec_sol);
267: return(0);
268: }
270: static PetscErrorCode SNESTSFormFunction_BDF(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
271: {
272: TS_BDF *bdf = (TS_BDF*)ts->data;
273: PetscInt k = bdf->k;
274: PetscReal t = bdf->time[0];
275: Vec V = bdf->vec_dot;
279: TSBDF_VecDot(ts,k,t,X,V,&bdf->shift);
280: /* F = Function(t,X,V) */
281: TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
282: return(0);
283: }
285: static PetscErrorCode SNESTSFormJacobian_BDF(PETSC_UNUSED SNES snes,
286: PETSC_UNUSED Vec X,
287: Mat J,Mat P,
288: TS ts)
289: {
290: TS_BDF *bdf = (TS_BDF*)ts->data;
291: PetscReal t = bdf->time[0];
292: Vec V = bdf->vec_dot;
293: PetscReal dVdX = bdf->shift;
297: /* J,P = Jacobian(t,X,V) */
298: TSComputeIJacobian(ts,t,X,V,dVdX,J,P,PETSC_FALSE);
299: return(0);
300: }
302: static PetscErrorCode TSReset_BDF(TS ts)
303: {
304: TS_BDF *bdf = (TS_BDF*)ts->data;
305: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
309: for (i=0; i<n; i++) {VecDestroy(&bdf->work[i]);}
310: VecDestroy(&bdf->vec_dot);
311: VecDestroy(&bdf->vec_lte);
312: return(0);
313: }
315: static PetscErrorCode TSDestroy_BDF(TS ts)
316: {
320: TSReset_BDF(ts);
321: PetscFree(ts->data);
322: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
323: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
324: return(0);
325: }
327: static PetscErrorCode TSSetUp_BDF(TS ts)
328: {
329: TS_BDF *bdf = (TS_BDF*)ts->data;
330: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
331: PetscReal low,high,two = 2;
335: bdf->k = bdf->n = 0;
336: for (i=0; i<n; i++) {VecDuplicate(ts->vec_sol,&bdf->work[i]);}
337: VecDuplicate(ts->vec_sol,&bdf->vec_dot);
338: VecDuplicate(ts->vec_sol,&bdf->vec_lte);
340: TSGetAdapt(ts,&ts->adapt);
341: TSAdaptCandidatesClear(ts->adapt);
342: TSAdaptGetClip(ts->adapt,&low,&high);
343: TSAdaptSetClip(ts->adapt,low,PetscMin(high,two));
345: TSGetSNES(ts,&ts->snes);
346: return(0);
347: }
349: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
350: {
351: TS_BDF *bdf = (TS_BDF*)ts->data;
355: PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
356: {
357: PetscBool flg;
358: PetscInt order = bdf->order;
359: PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
360: if (flg) {TSBDFSetOrder(ts,order);}
361: }
362: PetscOptionsTail();
363: return(0);
364: }
366: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
367: {
368: TS_BDF *bdf = (TS_BDF*)ts->data;
369: PetscBool iascii;
373: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
374: if (iascii) {
375: PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);
376: }
377: return(0);
378: }
380: /* ------------------------------------------------------------ */
382: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
383: {
384: TS_BDF *bdf = (TS_BDF*)ts->data;
387: if (order == bdf->order) return(0);
388: if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
389: bdf->order = order;
390: return(0);
391: }
393: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
394: {
395: TS_BDF *bdf = (TS_BDF*)ts->data;
398: *order = bdf->order;
399: return(0);
400: }
402: /* ------------------------------------------------------------ */
404: /*MC
405: TSBDF - DAE solver using BDF methods
407: Level: beginner
409: .seealso: TS, TSCreate(), TSSetType()
410: M*/
411: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
412: {
413: TS_BDF *bdf;
417: ts->ops->reset = TSReset_BDF;
418: ts->ops->destroy = TSDestroy_BDF;
419: ts->ops->view = TSView_BDF;
420: ts->ops->setup = TSSetUp_BDF;
421: ts->ops->setfromoptions = TSSetFromOptions_BDF;
422: ts->ops->step = TSStep_BDF;
423: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
424: ts->ops->rollback = TSRollBack_BDF;
425: ts->ops->interpolate = TSInterpolate_BDF;
426: ts->ops->snesfunction = SNESTSFormFunction_BDF;
427: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
428: ts->default_adapt_type = TSADAPTBASIC;
430: ts->usessnes = PETSC_TRUE;
432: PetscNewLog(ts,&bdf);
433: ts->data = (void*)bdf;
435: bdf->order = 2;
436: bdf->status = TS_STEP_COMPLETE;
438: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
439: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
440: return(0);
441: }
443: /* ------------------------------------------------------------ */
445: /*@
446: TSBDFSetOrder - Set the order of the BDF method
448: Logically Collective on TS
450: Input Parameter:
451: + ts - timestepping context
452: - order - order of the method
454: Options Database:
455: . -ts_bdf_order <order>
457: Level: intermediate
459: @*/
460: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
461: {
467: PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
468: return(0);
469: }
471: /*@
472: TSBDFGetOrder - Get the order of the BDF method
474: Not Collective
476: Input Parameter:
477: . ts - timestepping context
479: Output Parameter:
480: . order - order of the method
482: Level: intermediate
484: @*/
485: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
486: {
492: PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
493: return(0);
494: }