Actual source code: bdf.c
1: /*
2: Code for timestepping with BDF methods
3: */
4: #include <petsc/private/tsimpl.h>
5: #include <petscdm.h>
7: static PetscBool cited = PETSC_FALSE;
8: static const char citation[] =
9: "@book{Brenan1995,\n"
10: " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n"
11: " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n"
12: " publisher = {Society for Industrial and Applied Mathematics},\n"
13: " year = {1995},\n"
14: " doi = {10.1137/1.9781611971224},\n}\n";
16: typedef struct {
17: PetscInt k,n;
18: PetscReal time[6+2];
19: Vec work[6+2];
20: Vec tvwork[6+2];
21: PetscReal shift;
22: Vec vec_dot; /* Xdot when !transientvar, else Cdot where C(X) is the transient variable. */
23: Vec vec_wrk;
24: Vec vec_lte;
26: PetscBool transientvar;
27: PetscInt order;
28: TSStepStatus status;
29: } TS_BDF;
31: /* Compute Lagrange polynomials on T[:n] evaluated at t.
32: * If one has data (T[i], Y[i]), then the interpolation/extrapolation is f(t) = \sum_i L[i]*Y[i].
33: */
34: PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
35: {
36: PetscInt k,j;
37: for (k=0; k<n; k++)
38: for (L[k]=1, j=0; j<n; j++)
39: if (j != k)
40: L[k] *= (t - T[j])/(T[k] - T[j]);
41: }
43: PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])
44: {
45: PetscInt k,j,i;
46: for (k=0; k<n; k++)
47: for (dL[k]=0, j=0; j<n; j++)
48: if (j != k) {
49: PetscReal L = 1/(T[k] - T[j]);
50: for (i=0; i<n; i++)
51: if (i != j && i != k)
52: L *= (t - T[i])/(T[k] - T[i]);
53: dL[k] += L;
54: }
55: }
57: static PetscErrorCode TSBDF_GetVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
58: {
59: TS_BDF *bdf = (TS_BDF*)ts->data;
63: if (dm && dm != ts->dm) {
64: DMGetNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
65: DMGetNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
66: } else {
67: *Xdot = bdf->vec_dot;
68: *Ydot = bdf->vec_wrk;
69: }
70: return(0);
71: }
73: static PetscErrorCode TSBDF_RestoreVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
74: {
75: TS_BDF *bdf = (TS_BDF*)ts->data;
79: if (dm && dm != ts->dm) {
80: DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
81: DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
82: } else {
83: if (*Xdot != bdf->vec_dot) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache");
84: if (*Ydot != bdf->vec_wrk) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache");
85: *Xdot = NULL;
86: *Ydot = NULL;
87: }
88: return(0);
89: }
91: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine,DM coarse,void *ctx)
92: {
94: return(0);
95: }
97: static PetscErrorCode DMRestrictHook_TSBDF(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
98: {
99: TS ts = (TS)ctx;
100: Vec Ydot,Ydot_c;
101: Vec Xdot,Xdot_c;
105: TSBDF_GetVecs(ts,fine,&Xdot,&Ydot);
106: TSBDF_GetVecs(ts,coarse,&Xdot_c,&Ydot_c);
108: MatRestrict(restrct,Ydot,Ydot_c);
109: VecPointwiseMult(Ydot_c,rscale,Ydot_c);
111: TSBDF_RestoreVecs(ts,fine,&Xdot,&Ydot);
112: TSBDF_RestoreVecs(ts,coarse,&Xdot_c,&Ydot_c);
113: return(0);
114: }
116: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
117: {
118: TS_BDF *bdf = (TS_BDF*)ts->data;
119: PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
120: Vec tail = bdf->work[n-1],tvtail = bdf->tvwork[n-1];
124: for (i=n-1; i>=2; i--) {
125: bdf->time[i] = bdf->time[i-1];
126: bdf->work[i] = bdf->work[i-1];
127: bdf->tvwork[i] = bdf->tvwork[i-1];
128: }
129: bdf->n = PetscMin(bdf->n+1,n-1);
130: bdf->time[1] = t;
131: bdf->work[1] = tail;
132: bdf->tvwork[1] = tvtail;
133: VecCopy(X,tail);
134: TSComputeTransientVariable(ts,tail,tvtail);
135: return(0);
136: }
138: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
139: {
140: TS_BDF *bdf = (TS_BDF*)ts->data;
141: PetscInt i,n = order+1;
142: PetscReal *time = bdf->time;
143: Vec *vecs = bdf->work;
144: PetscScalar a[8],b[8],alpha[8];
148: LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
149: LagrangeBasisDers(n+1,time[0],time,b);
150: for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
151: VecZeroEntries(lte);
152: VecMAXPY(lte,n+1,alpha,vecs);
153: return(0);
154: }
156: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
157: {
158: TS_BDF *bdf = (TS_BDF*)ts->data;
159: PetscInt n = order+1;
160: PetscReal *time = bdf->time+1;
161: Vec *vecs = bdf->work+1;
162: PetscScalar alpha[7];
166: n = PetscMin(n,bdf->n);
167: LagrangeBasisVals(n,t,time,alpha);
168: VecZeroEntries(X);
169: VecMAXPY(X,n,alpha,vecs);
170: return(0);
171: }
173: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
174: {
175: TS_BDF *bdf = (TS_BDF*)ts->data;
176: PetscInt n = order+1;
177: PetscReal *time = bdf->time;
178: Vec *vecs = bdf->work;
179: PetscScalar alpha[7];
183: LagrangeBasisVals(n,t,time,alpha);
184: VecZeroEntries(X);
185: VecMAXPY(X,n,alpha,vecs);
186: return(0);
187: }
189: /* Compute the affine term V0 such that Xdot = shift*X + V0.
190: *
191: * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork.
192: */
193: static PetscErrorCode TSBDF_PreSolve(TS ts)
194: {
195: TS_BDF *bdf = (TS_BDF*)ts->data;
196: PetscInt i,n = PetscMax(bdf->k,1) + 1;
197: Vec V,V0;
198: Vec vecs[7];
199: PetscScalar alpha[7];
203: TSBDF_GetVecs(ts,NULL,&V,&V0);
204: LagrangeBasisDers(n,bdf->time[0],bdf->time,alpha);
205: for (i=1; i<n; i++) {
206: vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i];
207: }
208: VecZeroEntries(V0);
209: VecMAXPY(V0,n-1,alpha+1,vecs+1);
210: bdf->shift = PetscRealPart(alpha[0]);
211: TSBDF_RestoreVecs(ts,NULL,&V,&V0);
212: return(0);
213: }
215: static PetscErrorCode TSBDF_SNESSolve(TS ts,Vec b,Vec x)
216: {
217: PetscInt nits,lits;
221: TSBDF_PreSolve(ts);
222: SNESSolve(ts->snes,b,x);
223: SNESGetIterationNumber(ts->snes,&nits);
224: SNESGetLinearSolveIterations(ts->snes,&lits);
225: ts->snes_its += nits; ts->ksp_its += lits;
226: return(0);
227: }
229: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
230: {
231: TS_BDF *bdf = (TS_BDF*)ts->data;
235: bdf->k = 1; bdf->n = 0;
236: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
238: bdf->time[0] = ts->ptime + ts->time_step/2;
239: VecCopy(bdf->work[1],bdf->work[0]);
240: TSPreStage(ts,bdf->time[0]);
241: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
242: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
243: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
244: if (!*accept) return(0);
246: bdf->k = PetscMin(2,bdf->order); bdf->n++;
247: VecCopy(bdf->work[0],bdf->work[2]);
248: bdf->time[2] = bdf->time[0];
249: TSComputeTransientVariable(ts,bdf->work[2],bdf->tvwork[2]);
250: return(0);
251: }
253: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
255: static PetscErrorCode TSStep_BDF(TS ts)
256: {
257: TS_BDF *bdf = (TS_BDF*)ts->data;
258: PetscInt rejections = 0;
259: PetscBool stageok,accept = PETSC_TRUE;
260: PetscReal next_time_step = ts->time_step;
264: PetscCitationsRegister(citation,&cited);
266: if (!ts->steprollback && !ts->steprestart) {
267: bdf->k = PetscMin(bdf->k+1,bdf->order);
268: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
269: }
271: bdf->status = TS_STEP_INCOMPLETE;
272: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
274: if (ts->steprestart) {
275: TSBDF_Restart(ts,&stageok);
276: if (!stageok) goto reject_step;
277: }
279: bdf->time[0] = ts->ptime + ts->time_step;
280: TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
281: TSPreStage(ts,bdf->time[0]);
282: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
283: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
284: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
285: if (!stageok) goto reject_step;
287: bdf->status = TS_STEP_PENDING;
288: TSAdaptCandidatesClear(ts->adapt);
289: TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
290: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
291: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
292: if (!accept) { ts->time_step = next_time_step; goto reject_step; }
294: VecCopy(bdf->work[0],ts->vec_sol);
295: ts->ptime += ts->time_step;
296: ts->time_step = next_time_step;
297: break;
299: reject_step:
300: ts->reject++; accept = PETSC_FALSE;
301: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
302: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
303: ts->reason = TS_DIVERGED_STEP_REJECTED;
304: }
305: }
306: return(0);
307: }
309: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
310: {
311: TS_BDF *bdf = (TS_BDF*)ts->data;
315: TSBDF_Interpolate(ts,bdf->k,t,X);
316: return(0);
317: }
319: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
320: {
321: TS_BDF *bdf = (TS_BDF*)ts->data;
322: PetscInt k = bdf->k;
323: PetscReal wltea,wlter;
324: Vec X = bdf->work[0], Y = bdf->vec_lte;
328: k = PetscMin(k,bdf->n-1);
329: TSBDF_VecLTE(ts,k,Y);
330: VecAXPY(Y,1,X);
331: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
332: if (order) *order = k + 1;
333: return(0);
334: }
336: static PetscErrorCode TSRollBack_BDF(TS ts)
337: {
338: TS_BDF *bdf = (TS_BDF*)ts->data;
342: VecCopy(bdf->work[1],ts->vec_sol);
343: return(0);
344: }
346: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes,Vec X,Vec F,TS ts)
347: {
348: TS_BDF *bdf = (TS_BDF*)ts->data;
349: DM dm, dmsave = ts->dm;
350: PetscReal t = bdf->time[0];
351: PetscReal shift = bdf->shift;
352: Vec V,V0;
356: SNESGetDM(snes,&dm);
357: TSBDF_GetVecs(ts,dm,&V,&V0);
358: if (bdf->transientvar) { /* shift*C(X) + V0 */
359: TSComputeTransientVariable(ts,X,V);
360: VecAYPX(V,shift,V0);
361: } else { /* shift*X + V0 */
362: VecWAXPY(V,shift,X,V0);
363: }
365: /* F = Function(t,X,V) */
366: ts->dm = dm;
367: TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
368: ts->dm = dmsave;
370: TSBDF_RestoreVecs(ts,dm,&V,&V0);
371: return(0);
372: }
374: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes,Vec X,Mat J,Mat P,TS ts)
375: {
376: TS_BDF *bdf = (TS_BDF*)ts->data;
377: DM dm, dmsave = ts->dm;
378: PetscReal t = bdf->time[0];
379: PetscReal shift = bdf->shift;
380: Vec V,V0;
384: SNESGetDM(snes,&dm);
385: TSBDF_GetVecs(ts,dm,&V,&V0);
387: /* J,P = Jacobian(t,X,V) */
388: ts->dm = dm;
389: TSComputeIJacobian(ts,t,X,V,shift,J,P,PETSC_FALSE);
390: ts->dm = dmsave;
392: TSBDF_RestoreVecs(ts,dm,&V,&V0);
393: return(0);
394: }
396: static PetscErrorCode TSReset_BDF(TS ts)
397: {
398: TS_BDF *bdf = (TS_BDF*)ts->data;
399: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
403: bdf->k = bdf->n = 0;
404: for (i=0; i<n; i++) {
405: VecDestroy(&bdf->work[i]);
406: VecDestroy(&bdf->tvwork[i]);
407: }
408: VecDestroy(&bdf->vec_dot);
409: VecDestroy(&bdf->vec_wrk);
410: VecDestroy(&bdf->vec_lte);
411: if (ts->dm) {DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);}
412: return(0);
413: }
415: static PetscErrorCode TSDestroy_BDF(TS ts)
416: {
420: TSReset_BDF(ts);
421: PetscFree(ts->data);
422: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
423: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
424: return(0);
425: }
427: static PetscErrorCode TSSetUp_BDF(TS ts)
428: {
429: TS_BDF *bdf = (TS_BDF*)ts->data;
430: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
431: PetscReal low,high,two = 2;
435: TSHasTransientVariable(ts,&bdf->transientvar);
436: bdf->k = bdf->n = 0;
437: for (i=0; i<n; i++) {
438: VecDuplicate(ts->vec_sol,&bdf->work[i]);
439: if (i && bdf->transientvar) {
440: VecDuplicate(ts->vec_sol,&bdf->tvwork[i]);
441: }
442: }
443: VecDuplicate(ts->vec_sol,&bdf->vec_dot);
444: VecDuplicate(ts->vec_sol,&bdf->vec_wrk);
445: VecDuplicate(ts->vec_sol,&bdf->vec_lte);
446: TSGetDM(ts,&ts->dm);
447: DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);
449: TSGetAdapt(ts,&ts->adapt);
450: TSAdaptCandidatesClear(ts->adapt);
451: TSAdaptGetClip(ts->adapt,&low,&high);
452: TSAdaptSetClip(ts->adapt,low,PetscMin(high,two));
454: TSGetSNES(ts,&ts->snes);
455: return(0);
456: }
458: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
459: {
463: PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
464: {
465: PetscBool flg;
466: PetscInt order;
467: TSBDFGetOrder(ts,&order);
468: PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
469: if (flg) {TSBDFSetOrder(ts,order);}
470: }
471: PetscOptionsTail();
472: return(0);
473: }
475: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
476: {
477: TS_BDF *bdf = (TS_BDF*)ts->data;
478: PetscBool iascii;
482: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
483: if (iascii) {
484: PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);
485: }
486: return(0);
487: }
489: /* ------------------------------------------------------------ */
491: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
492: {
493: TS_BDF *bdf = (TS_BDF*)ts->data;
496: if (order == bdf->order) return(0);
497: if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
498: bdf->order = order;
499: return(0);
500: }
502: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
503: {
504: TS_BDF *bdf = (TS_BDF*)ts->data;
507: *order = bdf->order;
508: return(0);
509: }
511: /* ------------------------------------------------------------ */
513: /*MC
514: TSBDF - DAE solver using BDF methods
516: Level: beginner
518: .seealso: TS, TSCreate(), TSSetType()
519: M*/
520: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
521: {
522: TS_BDF *bdf;
526: ts->ops->reset = TSReset_BDF;
527: ts->ops->destroy = TSDestroy_BDF;
528: ts->ops->view = TSView_BDF;
529: ts->ops->setup = TSSetUp_BDF;
530: ts->ops->setfromoptions = TSSetFromOptions_BDF;
531: ts->ops->step = TSStep_BDF;
532: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
533: ts->ops->rollback = TSRollBack_BDF;
534: ts->ops->interpolate = TSInterpolate_BDF;
535: ts->ops->snesfunction = SNESTSFormFunction_BDF;
536: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
537: ts->default_adapt_type = TSADAPTBASIC;
539: ts->usessnes = PETSC_TRUE;
541: PetscNewLog(ts,&bdf);
542: ts->data = (void*)bdf;
544: bdf->status = TS_STEP_COMPLETE;
546: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
547: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
548: TSBDFSetOrder(ts,2);
549: return(0);
550: }
552: /* ------------------------------------------------------------ */
554: /*@
555: TSBDFSetOrder - Set the order of the BDF method
557: Logically Collective on TS
559: Input Parameters:
560: + ts - timestepping context
561: - order - order of the method
563: Options Database:
564: . -ts_bdf_order <order>
566: Level: intermediate
568: @*/
569: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
570: {
576: PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
577: return(0);
578: }
580: /*@
581: TSBDFGetOrder - Get the order of the BDF method
583: Not Collective
585: Input Parameter:
586: . ts - timestepping context
588: Output Parameter:
589: . order - order of the method
591: Level: intermediate
593: @*/
594: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
595: {
601: PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
602: return(0);
603: }