Actual source code: bdf.c
petsc-3.7.5 2017-01-01
1: /*
2: Code for timestepping with BDF methods
3: */
4: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
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: PetscBool adapt;
25: TSStepStatus status;
26: } TS_BDF;
29: PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
30: {
31: PetscInt k,j;
32: for (k=0; k<n; k++)
33: for (L[k]=1, j=0; j<n; j++)
34: if (j != k)
35: L[k] *= (t - T[j])/(T[k] - T[j]);
36: }
38: PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])
39: {
40: PetscInt k,j,i;
41: for (k=0; k<n; k++)
42: for (dL[k]=0, j=0; j<n; j++)
43: if (j != k) {
44: PetscReal L = 1/(T[k] - T[j]);
45: for (i=0; i<n; i++)
46: if (i != j && i != k)
47: L *= (t - T[i])/(T[k] - T[i]);
48: dL[k] += L;
49: }
50: }
54: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
55: {
56: TS_BDF *bdf = (TS_BDF*)ts->data;
57: PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
58: Vec tail = bdf->work[n-1];
62: for (i=n-1; i>=2; i--) {
63: bdf->time[i] = bdf->time[i-1];
64: bdf->work[i] = bdf->work[i-1];
65: }
66: bdf->n = PetscMin(bdf->n+1,n-1);
67: bdf->time[1] = t;
68: bdf->work[1] = tail;
69: VecCopy(X,tail);
70: return(0);
71: }
75: static PetscErrorCode TSBDF_VecDot(TS ts,PetscInt order,PetscReal t,Vec X,Vec Xdot,PetscReal *shift)
76: {
77: TS_BDF *bdf = (TS_BDF*)ts->data;
78: PetscInt i,n = order+1;
79: PetscReal time[7];
80: Vec vecs[7];
81: PetscScalar alpha[7];
85: n = PetscMax(n,2);
86: time[0] = t; for (i=1; i<n; i++) time[i] = bdf->time[i];
87: vecs[0] = X; for (i=1; i<n; i++) vecs[i] = bdf->work[i];
88: LagrangeBasisDers(n,t,time,alpha);
89: VecZeroEntries(Xdot);
90: VecMAXPY(Xdot,n,alpha,vecs);
91: if (shift) *shift = PetscRealPart(alpha[0]);
92: return(0);
93: }
97: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
98: {
99: TS_BDF *bdf = (TS_BDF*)ts->data;
100: PetscInt i,n = order+1;
101: PetscReal *time = bdf->time;
102: Vec *vecs = bdf->work;
103: PetscScalar a[8],b[8],alpha[8];
107: LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
108: LagrangeBasisDers(n+1,time[0],time,b);
109: for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
110: VecZeroEntries(lte);
111: VecMAXPY(lte,n+1,alpha,vecs);
112: return(0);
113: }
117: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
118: {
119: TS_BDF *bdf = (TS_BDF*)ts->data;
120: PetscInt n = order+1;
121: PetscReal *time = bdf->time+1;
122: Vec *vecs = bdf->work+1;
123: PetscScalar alpha[7];
127: n = PetscMin(n,bdf->n);
128: LagrangeBasisVals(n,t,time,alpha);
129: VecZeroEntries(X);
130: VecMAXPY(X,n,alpha,vecs);
131: return(0);
132: }
136: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
137: {
138: TS_BDF *bdf = (TS_BDF*)ts->data;
139: PetscInt n = order+1;
140: PetscReal *time = bdf->time;
141: Vec *vecs = bdf->work;
142: PetscScalar alpha[7];
146: LagrangeBasisVals(n,t,time,alpha);
147: VecZeroEntries(X);
148: VecMAXPY(X,n,alpha,vecs);
149: return(0);
150: }
154: static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
155: {
156: PetscInt nits,lits;
160: SNESSolve(ts->snes,b,x);
161: SNESGetIterationNumber(ts->snes,&nits);
162: SNESGetLinearSolveIterations(ts->snes,&lits);
163: ts->snes_its += nits; ts->ksp_its += lits;
164: return(0);
165: }
169: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
170: {
171: TS_BDF *bdf = (TS_BDF*)ts->data;
175: bdf->k = 1; bdf->n = 0;
176: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
178: bdf->time[0] = ts->ptime + ts->time_step/2;
179: VecCopy(bdf->work[1],bdf->work[0]);
180: TSPreStage(ts,bdf->time[0]);
181: TS_SNESSolve(ts,NULL,bdf->work[0]);
182: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
183: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
184: if (!*accept) return(0);
186: bdf->k = PetscMin(2,bdf->order); bdf->n++;
187: VecCopy(bdf->work[0],bdf->work[2]);
188: bdf->time[2] = bdf->time[0];
189: return(0);
190: }
192: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
196: static PetscErrorCode TSStep_BDF(TS ts)
197: {
198: TS_BDF *bdf = (TS_BDF*)ts->data;
199: PetscInt rejections = 0;
200: PetscBool stageok,accept = PETSC_TRUE;
201: PetscReal next_time_step = ts->time_step;
205: PetscCitationsRegister(citation,&cited);
207: if (!ts->steprollback && !ts->steprestart) {
208: bdf->k = PetscMin(bdf->k+1,bdf->order);
209: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
210: }
212: bdf->status = TS_STEP_INCOMPLETE;
213: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
215: if (ts->steprestart) {
216: TSBDF_Restart(ts,&stageok);
217: if (!stageok) goto reject_step;
218: }
220: bdf->time[0] = ts->ptime + ts->time_step;
221: TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
222: TSPreStage(ts,bdf->time[0]);
223: TS_SNESSolve(ts,NULL,bdf->work[0]);
224: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
225: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
226: if (!stageok) goto reject_step;
228: bdf->status = TS_STEP_PENDING;
229: TSAdaptCandidatesClear(ts->adapt);
230: TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
231: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
232: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
233: if (!accept) { ts->time_step = next_time_step; goto reject_step; }
235: VecCopy(bdf->work[0],ts->vec_sol);
236: ts->ptime += ts->time_step;
237: ts->time_step = next_time_step;
238: break;
240: reject_step:
241: ts->reject++; accept = PETSC_FALSE;
242: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
243: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
244: ts->reason = TS_DIVERGED_STEP_REJECTED;
245: }
246: }
247: return(0);
248: }
252: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
253: {
254: TS_BDF *bdf = (TS_BDF*)ts->data;
258: TSBDF_Interpolate(ts,bdf->k,t,X);
259: return(0);
260: }
264: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
265: {
266: TS_BDF *bdf = (TS_BDF*)ts->data;
267: PetscInt k = bdf->k;
268: Vec X = bdf->work[0], Y = bdf->vec_lte;
272: k = PetscMin(k,bdf->n-1);
273: TSBDF_VecLTE(ts,k,Y);
274: VecAXPY(Y,1,X);
275: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte);
276: if (order) *order = k + 1;
277: return(0);
278: }
282: static PetscErrorCode TSRollBack_BDF(TS ts)
283: {
284: TS_BDF *bdf = (TS_BDF*)ts->data;
288: VecCopy(bdf->work[1],ts->vec_sol);
289: return(0);
290: }
294: static PetscErrorCode SNESTSFormFunction_BDF(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
295: {
296: TS_BDF *bdf = (TS_BDF*)ts->data;
297: PetscInt k = bdf->k;
298: PetscReal t = bdf->time[0];
299: Vec V = bdf->vec_dot;
303: TSBDF_VecDot(ts,k,t,X,V,&bdf->shift);
304: /* F = Function(t,X,V) */
305: TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
306: return(0);
307: }
311: static PetscErrorCode SNESTSFormJacobian_BDF(PETSC_UNUSED SNES snes,
312: PETSC_UNUSED Vec X,
313: Mat J,Mat P,
314: TS ts)
315: {
316: TS_BDF *bdf = (TS_BDF*)ts->data;
317: PetscReal t = bdf->time[0];
318: Vec V = bdf->vec_dot;
319: PetscReal dVdX = bdf->shift;
323: /* J,P = Jacobian(t,X,V) */
324: TSComputeIJacobian(ts,t,X,V,dVdX,J,P,PETSC_FALSE);
325: return(0);
326: }
330: static PetscErrorCode TSReset_BDF(TS ts)
331: {
332: TS_BDF *bdf = (TS_BDF*)ts->data;
333: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
337: for (i=0; i<n; i++) {VecDestroy(&bdf->work[i]);}
338: VecDestroy(&bdf->vec_dot);
339: VecDestroy(&bdf->vec_lte);
340: return(0);
341: }
345: static PetscErrorCode TSDestroy_BDF(TS ts)
346: {
350: TSReset_BDF(ts);
351: PetscFree(ts->data);
352: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
353: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
354: PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",NULL);
355: return(0);
356: }
360: static PetscErrorCode TSSetUp_BDF(TS ts)
361: {
362: TS_BDF *bdf = (TS_BDF*)ts->data;
363: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
367: bdf->k = bdf->n = 0;
368: for (i=0; i<n; i++) {VecDuplicate(ts->vec_sol,&bdf->work[i]);}
369: VecDuplicate(ts->vec_sol,&bdf->vec_dot);
370: VecDuplicate(ts->vec_sol,&bdf->vec_lte);
372: TSGetAdapt(ts,&ts->adapt);
373: TSAdaptCandidatesClear(ts->adapt);
374: if (!bdf->adapt) {
375: TSAdaptSetType(ts->adapt,TSADAPTNONE);
376: } else {
377: PetscReal low,high;
378: TSAdaptBasicGetClip(ts->adapt,&low,&high);
379: high = PetscMin(high,2.0);
380: TSAdaptBasicSetClip(ts->adapt,low,high);
381: if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
382: ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
383: }
385: TSGetSNES(ts,&ts->snes);
386: return(0);
387: }
391: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
392: {
393: TS_BDF *bdf = (TS_BDF*)ts->data;
397: PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
398: {
399: PetscBool flg;
400: PetscInt order = bdf->order;
401: PetscBool adapt = bdf->adapt;
402: PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
403: if (flg) {TSBDFSetOrder(ts,order);}
404: PetscOptionsBool("-ts_bdf_adapt","Use time-step adaptivity with the BDF method","TSBDFUseAdapt",adapt,&adapt,&flg);
405: if (flg) {TSBDFUseAdapt(ts,adapt);}
406: }
407: PetscOptionsTail();
408: return(0);
409: }
413: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
414: {
415: TS_BDF *bdf = (TS_BDF*)ts->data;
416: PetscBool iascii;
420: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
421: if (iascii) {
422: PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);
423: }
424: if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
425: if (ts->snes) {SNESView(ts->snes,viewer);}
426: return(0);
427: }
429: /* ------------------------------------------------------------ */
433: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
434: {
435: TS_BDF *bdf = (TS_BDF*)ts->data;
438: if (order == bdf->order) return(0);
439: if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
440: bdf->order = order;
441: return(0);
442: }
446: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
447: {
448: TS_BDF *bdf = (TS_BDF*)ts->data;
451: *order = bdf->order;
452: return(0);
453: }
457: static PetscErrorCode TSBDFUseAdapt_BDF(TS ts,PetscBool use)
458: {
459: TS_BDF *bdf = (TS_BDF*)ts->data;
462: if (use == bdf->adapt) return(0);
463: if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()");
464: bdf->adapt = use;
465: return(0);
466: }
468: /* ------------------------------------------------------------ */
470: /*MC
471: TSBDF - DAE solver using BDF methods
473: Level: beginner
475: .seealso: TS, TSCreate(), TSSetType()
476: M*/
479: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
480: {
481: TS_BDF *bdf;
485: ts->ops->reset = TSReset_BDF;
486: ts->ops->destroy = TSDestroy_BDF;
487: ts->ops->view = TSView_BDF;
488: ts->ops->setup = TSSetUp_BDF;
489: ts->ops->setfromoptions = TSSetFromOptions_BDF;
490: ts->ops->step = TSStep_BDF;
491: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
492: ts->ops->rollback = TSRollBack_BDF;
493: ts->ops->interpolate = TSInterpolate_BDF;
494: ts->ops->snesfunction = SNESTSFormFunction_BDF;
495: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
497: PetscNewLog(ts,&bdf);
498: ts->data = (void*)bdf;
500: bdf->order = 2;
501: bdf->adapt = PETSC_FALSE;
502: bdf->status = TS_STEP_COMPLETE;
504: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
505: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
506: PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",TSBDFUseAdapt_BDF);
507: return(0);
508: }
510: /* ------------------------------------------------------------ */
514: /*@
515: TSBDFSetOrder - Set the order of the BDF method
517: Logically Collective on TS
519: Input Parameter:
520: + ts - timestepping context
521: - order - order of the method
523: Options Database:
524: . -ts_bdf_order <order>
526: Level: intermediate
528: @*/
529: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
530: {
536: PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
537: return(0);
538: }
542: /*@
543: TSBDFGetOrder - Get the order of the BDF method
545: Not Collective
547: Input Parameter:
548: . ts - timestepping context
550: Output Parameter:
551: . order - order of the method
553: Level: intermediate
555: @*/
556: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
557: {
563: PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
564: return(0);
565: }
569: /*@
570: TSBDFUseAdapt - Use time-step adaptivity with the BDF method
572: Logically Collective on TS
574: Input Parameter:
575: + ts - timestepping context
576: - use - flag to use adaptivity
578: Options Database:
579: . -ts_bdf_adapt
581: Level: intermediate
583: .seealso: TSAdapt
584: @*/
585: PetscErrorCode TSBDFUseAdapt(TS ts,PetscBool use)
586: {
592: PetscTryMethod(ts,"TSBDFUseAdapt_C",(TS,PetscBool),(ts,use));
593: return(0);
594: }