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