xref: /petsc/src/ts/impls/bdf/bdf.c (revision ba1e01c44f2f740c97c9026fe63e360558b7709b)
1 /*
2   Code for timestepping with BDF methods
3 */
4 #include <petsc/private/tsimpl.h>  /*I "petscts.h" I*/
5 
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";
14 
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;
22 
23   PetscInt     order;
24   PetscBool    adapt;
25   TSStepStatus status;
26 } TS_BDF;
27 
28 
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 }
37 
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 }
51 
52 #undef __FUNCT__
53 #define __FUNCT__ "TSBDF_Advance"
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];
59   PetscErrorCode ierr;
60 
61   PetscFunctionBegin;
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   ierr = VecCopy(X,tail);CHKERRQ(ierr);
70   PetscFunctionReturn(0);
71 }
72 
73 #undef __FUNCT__
74 #define __FUNCT__ "TSBDF_VecDot"
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];
82   PetscErrorCode ierr;
83 
84   PetscFunctionBegin;
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   ierr = VecZeroEntries(Xdot);CHKERRQ(ierr);
90   ierr = VecMAXPY(Xdot,n,alpha,vecs);CHKERRQ(ierr);
91   if (shift) *shift = PetscRealPart(alpha[0]);
92   PetscFunctionReturn(0);
93 }
94 
95 #undef __FUNCT__
96 #define __FUNCT__ "TSBDF_VecLTE"
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];
104   PetscErrorCode ierr;
105 
106   PetscFunctionBegin;
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   ierr = VecZeroEntries(lte);CHKERRQ(ierr);
111   ierr = VecMAXPY(lte,n+1,alpha,vecs);CHKERRQ(ierr);
112   PetscFunctionReturn(0);
113 }
114 
115 #undef __FUNCT__
116 #define __FUNCT__ "TSBDF_Extrapolate"
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];
124   PetscErrorCode ierr;
125 
126   PetscFunctionBegin;
127   n = PetscMin(n,bdf->n);
128   LagrangeBasisVals(n,t,time,alpha);
129   ierr = VecZeroEntries(X);CHKERRQ(ierr);
130   ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr);
131   PetscFunctionReturn(0);
132 }
133 
134 #undef __FUNCT__
135 #define __FUNCT__ "TSBDF_Interpolate"
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];
143   PetscErrorCode ierr;
144 
145   PetscFunctionBegin;
146   LagrangeBasisVals(n,t,time,alpha);
147   ierr = VecZeroEntries(X);CHKERRQ(ierr);
148   ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr);
149   PetscFunctionReturn(0);
150 }
151 
152 #undef __FUNCT__
153 #define __FUNCT__ "TS_SNESSolve"
154 static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
155 {
156   PetscInt       nits,lits;
157   PetscErrorCode ierr;
158 
159   PetscFunctionBegin;
160   ierr = SNESSolve(ts->snes,b,x);CHKERRQ(ierr);
161   ierr = SNESGetIterationNumber(ts->snes,&nits);CHKERRQ(ierr);
162   ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr);
163   ts->snes_its += nits; ts->ksp_its += lits;
164   PetscFunctionReturn(0);
165 }
166 
167 #undef __FUNCT__
168 #define __FUNCT__ "TSBDF_Restart"
169 static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
170 {
171   TS_BDF         *bdf = (TS_BDF*)ts->data;
172   PetscErrorCode ierr;
173 
174   PetscFunctionBegin;
175   bdf->k = 1; bdf->n = 0;
176   ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr);
177 
178   bdf->time[0] = ts->ptime + ts->time_step/2;
179   ierr = VecCopy(bdf->work[1],bdf->work[0]);CHKERRQ(ierr);
180   ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr);
181   ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr);
182   ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr);
183   ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);CHKERRQ(ierr);
184   if (!*accept) PetscFunctionReturn(0);
185 
186   bdf->k = PetscMin(2,bdf->order); bdf->n++;
187   ierr = VecCopy(bdf->work[0],bdf->work[2]);CHKERRQ(ierr);
188   bdf->time[2] = bdf->time[0];
189   PetscFunctionReturn(0);
190 }
191 
192 static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
193 
194 #undef __FUNCT__
195 #define __FUNCT__ "TSStep_BDF"
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;
202   PetscErrorCode ierr;
203 
204   PetscFunctionBegin;
205   ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
206 
207   if (!ts->steprollback && !ts->steprestart) {
208     bdf->k = PetscMin(bdf->k+1,bdf->order);
209     ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr);
210   }
211 
212   bdf->status = TS_STEP_INCOMPLETE;
213   while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
214 
215     if (ts->steprestart) {
216       ierr = TSBDF_Restart(ts,&stageok);CHKERRQ(ierr);
217       if (!stageok) goto reject_step;
218     }
219 
220     bdf->time[0] = ts->ptime + ts->time_step;
221     ierr = TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);CHKERRQ(ierr);
222     ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr);
223     ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr);
224     ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr);
225     ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);CHKERRQ(ierr);
226     if (!stageok) goto reject_step;
227 
228     bdf->status = TS_STEP_PENDING;
229     ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr);
230     ierr = TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);CHKERRQ(ierr);
231     ierr = TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr);
232     bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
233     if (!accept) { ts->time_step = next_time_step; goto reject_step; }
234 
235     ierr = VecCopy(bdf->work[0],ts->vec_sol);CHKERRQ(ierr);
236     ts->ptime += ts->time_step;
237     ts->time_step = next_time_step;
238     break;
239 
240   reject_step:
241     ts->reject++; accept = PETSC_FALSE;
242     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
243       ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr);
244       ts->reason = TS_DIVERGED_STEP_REJECTED;
245     }
246   }
247   PetscFunctionReturn(0);
248 }
249 
250 #undef __FUNCT__
251 #define __FUNCT__ "TSInterpolate_BDF"
252 static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
253 {
254   TS_BDF         *bdf = (TS_BDF*)ts->data;
255   PetscErrorCode ierr;
256 
257   PetscFunctionBegin;
258   ierr = TSBDF_Interpolate(ts,bdf->k,t,X);CHKERRQ(ierr);
259   PetscFunctionReturn(0);
260 }
261 
262 #undef __FUNCT__
263 #define __FUNCT__ "TSEvaluateWLTE_BDF"
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   PetscReal      wltea,wlter;
269   Vec            X = bdf->work[0], Y = bdf->vec_lte;
270   PetscErrorCode ierr;
271 
272   PetscFunctionBegin;
273   k = PetscMin(k,bdf->n-1);
274   ierr = TSBDF_VecLTE(ts,k,Y);CHKERRQ(ierr);
275   ierr = VecAXPY(Y,1,X);CHKERRQ(ierr);
276   ierr = TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);CHKERRQ(ierr);
277   if (order) *order = k + 1;
278   PetscFunctionReturn(0);
279 }
280 
281 #undef __FUNCT__
282 #define __FUNCT__ "TSRollBack_BDF"
283 static PetscErrorCode TSRollBack_BDF(TS ts)
284 {
285   TS_BDF         *bdf = (TS_BDF*)ts->data;
286   PetscErrorCode ierr;
287 
288   PetscFunctionBegin;
289   ierr = VecCopy(bdf->work[1],ts->vec_sol);CHKERRQ(ierr);
290   PetscFunctionReturn(0);
291 }
292 
293 #undef __FUNCT__
294 #define __FUNCT__ "SNESTSFormFunction_BDF"
295 static PetscErrorCode SNESTSFormFunction_BDF(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
296 {
297   TS_BDF         *bdf = (TS_BDF*)ts->data;
298   PetscInt       k = bdf->k;
299   PetscReal      t = bdf->time[0];
300   Vec            V = bdf->vec_dot;
301   PetscErrorCode ierr;
302 
303   PetscFunctionBegin;
304   ierr = TSBDF_VecDot(ts,k,t,X,V,&bdf->shift);CHKERRQ(ierr);
305   /* F = Function(t,X,V) */
306   ierr = TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);CHKERRQ(ierr);
307   PetscFunctionReturn(0);
308 }
309 
310 #undef __FUNCT__
311 #define __FUNCT__ "SNESTSFormJacobian_BDF"
312 static PetscErrorCode SNESTSFormJacobian_BDF(PETSC_UNUSED SNES snes,
313                                              PETSC_UNUSED Vec X,
314                                              Mat J,Mat P,
315                                              TS ts)
316 {
317   TS_BDF         *bdf = (TS_BDF*)ts->data;
318   PetscReal      t = bdf->time[0];
319   Vec            V = bdf->vec_dot;
320   PetscReal      dVdX = bdf->shift;
321   PetscErrorCode ierr;
322 
323   PetscFunctionBegin;
324   /* J,P = Jacobian(t,X,V) */
325   ierr = TSComputeIJacobian(ts,t,X,V,dVdX,J,P,PETSC_FALSE);CHKERRQ(ierr);
326   PetscFunctionReturn(0);
327 }
328 
329 #undef __FUNCT__
330 #define __FUNCT__ "TSReset_BDF"
331 static PetscErrorCode TSReset_BDF(TS ts)
332 {
333   TS_BDF         *bdf = (TS_BDF*)ts->data;
334   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);
335   PetscErrorCode ierr;
336 
337   PetscFunctionBegin;
338   for (i=0; i<n; i++) {ierr = VecDestroy(&bdf->work[i]);CHKERRQ(ierr);}
339   ierr = VecDestroy(&bdf->vec_dot);CHKERRQ(ierr);
340   ierr = VecDestroy(&bdf->vec_lte);CHKERRQ(ierr);
341   PetscFunctionReturn(0);
342 }
343 
344 #undef __FUNCT__
345 #define __FUNCT__ "TSDestroy_BDF"
346 static PetscErrorCode TSDestroy_BDF(TS ts)
347 {
348   PetscErrorCode ierr;
349 
350   PetscFunctionBegin;
351   ierr = TSReset_BDF(ts);CHKERRQ(ierr);
352   ierr = PetscFree(ts->data);CHKERRQ(ierr);
353   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);CHKERRQ(ierr);
354   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);CHKERRQ(ierr);
355   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",NULL);CHKERRQ(ierr);
356   PetscFunctionReturn(0);
357 }
358 
359 #undef __FUNCT__
360 #define __FUNCT__ "TSSetUp_BDF"
361 static PetscErrorCode TSSetUp_BDF(TS ts)
362 {
363   TS_BDF         *bdf = (TS_BDF*)ts->data;
364   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);
365   PetscErrorCode ierr;
366 
367   PetscFunctionBegin;
368   bdf->k = bdf->n = 0;
369   for (i=0; i<n; i++) {ierr = VecDuplicate(ts->vec_sol,&bdf->work[i]);CHKERRQ(ierr);}
370   ierr = VecDuplicate(ts->vec_sol,&bdf->vec_dot);CHKERRQ(ierr);
371   ierr = VecDuplicate(ts->vec_sol,&bdf->vec_lte);CHKERRQ(ierr);
372 
373   ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr);
374   ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr);
375   if (!bdf->adapt) {
376     ierr = TSAdaptSetType(ts->adapt,TSADAPTNONE);CHKERRQ(ierr);
377   } else {
378     PetscReal low,high;
379     ierr = TSAdaptBasicGetClip(ts->adapt,&low,&high);CHKERRQ(ierr);
380     high = PetscMin(high,2.0);
381     ierr = TSAdaptBasicSetClip(ts->adapt,low,high);CHKERRQ(ierr);
382     if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
383       ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
384   }
385 
386   ierr = TSGetSNES(ts,&ts->snes);CHKERRQ(ierr);
387   PetscFunctionReturn(0);
388 }
389 
390 #undef __FUNCT__
391 #define __FUNCT__ "TSSetFromOptions_BDF"
392 static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
393 {
394   TS_BDF         *bdf = (TS_BDF*)ts->data;
395   PetscErrorCode ierr;
396 
397   PetscFunctionBegin;
398   ierr = PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");CHKERRQ(ierr);
399   {
400     PetscBool flg;
401     PetscInt  order = bdf->order;
402     PetscBool adapt = bdf->adapt;
403     ierr = PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);CHKERRQ(ierr);
404     if (flg) {ierr = TSBDFSetOrder(ts,order);CHKERRQ(ierr);}
405     ierr = PetscOptionsBool("-ts_bdf_adapt","Use time-step adaptivity with the BDF method","TSBDFUseAdapt",adapt,&adapt,&flg);CHKERRQ(ierr);
406     if (flg) {ierr = TSBDFUseAdapt(ts,adapt);CHKERRQ(ierr);}
407   }
408   ierr = PetscOptionsTail();CHKERRQ(ierr);
409   PetscFunctionReturn(0);
410 }
411 
412 #undef __FUNCT__
413 #define __FUNCT__ "TSView_BDF"
414 static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
415 {
416   TS_BDF         *bdf = (TS_BDF*)ts->data;
417   PetscBool      iascii;
418   PetscErrorCode ierr;
419 
420   PetscFunctionBegin;
421   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
422   if (iascii)    {
423     ierr = PetscViewerASCIIPrintf(viewer,"  Order=%D\n",bdf->order);CHKERRQ(ierr);
424   }
425   if (ts->adapt) {ierr = TSAdaptView(ts->adapt,viewer);CHKERRQ(ierr);}
426   if (ts->snes)  {ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr);}
427   PetscFunctionReturn(0);
428 }
429 
430 /* ------------------------------------------------------------ */
431 
432 #undef __FUNCT__
433 #define __FUNCT__ "TSBDFSetOrder_BDF"
434 static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
435 {
436   TS_BDF *bdf = (TS_BDF*)ts->data;
437 
438   PetscFunctionBegin;
439   if (order == bdf->order) PetscFunctionReturn(0);
440   if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
441   bdf->order = order;
442   PetscFunctionReturn(0);
443 }
444 
445 #undef __FUNCT__
446 #define __FUNCT__ "TSBDFGetOrder_BDF"
447 static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
448 {
449   TS_BDF *bdf = (TS_BDF*)ts->data;
450 
451   PetscFunctionBegin;
452   *order = bdf->order;
453   PetscFunctionReturn(0);
454 }
455 
456 #undef __FUNCT__
457 #define __FUNCT__ "TSBDFUseAdapt_BDF"
458 static PetscErrorCode TSBDFUseAdapt_BDF(TS ts,PetscBool use)
459 {
460   TS_BDF *bdf = (TS_BDF*)ts->data;
461 
462   PetscFunctionBegin;
463   if (use == bdf->adapt) PetscFunctionReturn(0);
464   if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()");
465   bdf->adapt = use;
466   PetscFunctionReturn(0);
467 }
468 
469 /* ------------------------------------------------------------ */
470 
471 /*MC
472       TSBDF - DAE solver using BDF methods
473 
474   Level: beginner
475 
476 .seealso:  TS, TSCreate(), TSSetType()
477 M*/
478 #undef __FUNCT__
479 #define __FUNCT__ "TSCreate_BDF"
480 PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
481 {
482   TS_BDF         *bdf;
483   PetscErrorCode ierr;
484 
485   PetscFunctionBegin;
486   ts->ops->reset          = TSReset_BDF;
487   ts->ops->destroy        = TSDestroy_BDF;
488   ts->ops->view           = TSView_BDF;
489   ts->ops->setup          = TSSetUp_BDF;
490   ts->ops->setfromoptions = TSSetFromOptions_BDF;
491   ts->ops->step           = TSStep_BDF;
492   ts->ops->evaluatewlte   = TSEvaluateWLTE_BDF;
493   ts->ops->rollback       = TSRollBack_BDF;
494   ts->ops->interpolate    = TSInterpolate_BDF;
495   ts->ops->snesfunction   = SNESTSFormFunction_BDF;
496   ts->ops->snesjacobian   = SNESTSFormJacobian_BDF;
497 
498   ierr = PetscNewLog(ts,&bdf);CHKERRQ(ierr);
499   ts->data = (void*)bdf;
500 
501   bdf->order  = 2;
502   bdf->adapt  = PETSC_FALSE;
503   bdf->status = TS_STEP_COMPLETE;
504 
505   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);CHKERRQ(ierr);
506   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);CHKERRQ(ierr);
507   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",TSBDFUseAdapt_BDF);CHKERRQ(ierr);
508   PetscFunctionReturn(0);
509 }
510 
511 /* ------------------------------------------------------------ */
512 
513 #undef __FUNCT__
514 #define __FUNCT__ "TSBDFSetOrder"
515 /*@
516   TSBDFSetOrder - Set the order of the BDF method
517 
518   Logically Collective on TS
519 
520   Input Parameter:
521 +  ts - timestepping context
522 -  order - order of the method
523 
524   Options Database:
525 .  -ts_bdf_order <order>
526 
527   Level: intermediate
528 
529 @*/
530 PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
531 {
532   PetscErrorCode ierr;
533 
534   PetscFunctionBegin;
535   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
536   PetscValidLogicalCollectiveInt(ts,order,2);
537   ierr = PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));CHKERRQ(ierr);
538   PetscFunctionReturn(0);
539 }
540 
541 #undef __FUNCT__
542 #define __FUNCT__ "TSBDFGetOrder"
543 /*@
544   TSBDFGetOrder - Get the order of the BDF method
545 
546   Not Collective
547 
548   Input Parameter:
549 .  ts - timestepping context
550 
551   Output Parameter:
552 .  order - order of the method
553 
554   Level: intermediate
555 
556 @*/
557 PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
558 {
559   PetscErrorCode ierr;
560 
561   PetscFunctionBegin;
562   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
563   PetscValidIntPointer(order,2);
564   ierr = PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));CHKERRQ(ierr);
565   PetscFunctionReturn(0);
566 }
567 
568 #undef __FUNCT__
569 #define __FUNCT__ "TSBDFUseAdapt"
570 /*@
571   TSBDFUseAdapt - Use time-step adaptivity with the BDF method
572 
573   Logically Collective on TS
574 
575   Input Parameter:
576 +  ts - timestepping context
577 -  use - flag to use adaptivity
578 
579   Options Database:
580 .  -ts_bdf_adapt
581 
582   Level: intermediate
583 
584 .seealso: TSAdapt
585 @*/
586 PetscErrorCode TSBDFUseAdapt(TS ts,PetscBool use)
587 {
588   PetscErrorCode ierr;
589 
590   PetscFunctionBegin;
591   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
592   PetscValidLogicalCollectiveBool(ts,use,2);
593   ierr = PetscTryMethod(ts,"TSBDFUseAdapt_C",(TS,PetscBool),(ts,use));CHKERRQ(ierr);
594   PetscFunctionReturn(0);
595 }
596