xref: /petsc/src/ts/impls/bdf/bdf.c (revision 7453f77509fc006cbcc7de2dd772f60dc49feac5)
1211a84d6SLisandro Dalcin /*
2211a84d6SLisandro Dalcin   Code for timestepping with BDF methods
3211a84d6SLisandro Dalcin */
4211a84d6SLisandro Dalcin #include <petsc/private/tsimpl.h>  /*I "petscts.h" I*/
5211a84d6SLisandro Dalcin 
6211a84d6SLisandro Dalcin static PetscBool  cited = PETSC_FALSE;
7211a84d6SLisandro Dalcin static const char citation[] =
8211a84d6SLisandro Dalcin   "@book{Brenan1995,\n"
9211a84d6SLisandro Dalcin   "  title     = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n"
10211a84d6SLisandro Dalcin   "  author    = {Brenan, K. and Campbell, S. and Petzold, L.},\n"
11211a84d6SLisandro Dalcin   "  publisher = {Society for Industrial and Applied Mathematics},\n"
12211a84d6SLisandro Dalcin   "  year      = {1995},\n"
13211a84d6SLisandro Dalcin   "  doi       = {10.1137/1.9781611971224},\n}\n";
14211a84d6SLisandro Dalcin 
15211a84d6SLisandro Dalcin typedef struct {
16211a84d6SLisandro Dalcin   PetscInt  k,n;
17211a84d6SLisandro Dalcin   PetscReal time[6+2];
18211a84d6SLisandro Dalcin   Vec       work[6+2];
19211a84d6SLisandro Dalcin   PetscReal shift;
20211a84d6SLisandro Dalcin   Vec       vec_dot;
21211a84d6SLisandro Dalcin   Vec       vec_lte;
22211a84d6SLisandro Dalcin 
23211a84d6SLisandro Dalcin   PetscInt     order;
24211a84d6SLisandro Dalcin   PetscBool    adapt;
25211a84d6SLisandro Dalcin   TSStepStatus status;
26211a84d6SLisandro Dalcin } TS_BDF;
27211a84d6SLisandro Dalcin 
28211a84d6SLisandro Dalcin 
29211a84d6SLisandro Dalcin PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
30211a84d6SLisandro Dalcin {
31211a84d6SLisandro Dalcin   PetscInt k,j;
32211a84d6SLisandro Dalcin   for (k=0; k<n; k++)
33211a84d6SLisandro Dalcin     for (L[k]=1, j=0; j<n; j++)
34211a84d6SLisandro Dalcin       if (j != k)
35211a84d6SLisandro Dalcin         L[k] *= (t - T[j])/(T[k] - T[j]);
36211a84d6SLisandro Dalcin }
37211a84d6SLisandro Dalcin 
38211a84d6SLisandro Dalcin PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])
39211a84d6SLisandro Dalcin {
40211a84d6SLisandro Dalcin   PetscInt  k,j,i;
41211a84d6SLisandro Dalcin   for (k=0; k<n; k++)
42211a84d6SLisandro Dalcin     for (dL[k]=0, j=0; j<n; j++)
43211a84d6SLisandro Dalcin       if (j != k) {
44211a84d6SLisandro Dalcin         PetscReal L = 1/(T[k] - T[j]);
45211a84d6SLisandro Dalcin         for (i=0; i<n; i++)
46211a84d6SLisandro Dalcin           if (i != j && i != k)
47211a84d6SLisandro Dalcin             L *= (t - T[i])/(T[k] - T[i]);
48211a84d6SLisandro Dalcin         dL[k] += L;
49211a84d6SLisandro Dalcin       }
50211a84d6SLisandro Dalcin }
51211a84d6SLisandro Dalcin 
52211a84d6SLisandro Dalcin #undef __FUNCT__
53211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Advance"
54211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
55211a84d6SLisandro Dalcin {
56211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
57211a84d6SLisandro Dalcin   PetscInt       i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
58211a84d6SLisandro Dalcin   Vec            tail = bdf->work[n-1];
59211a84d6SLisandro Dalcin   PetscErrorCode ierr;
60211a84d6SLisandro Dalcin 
61211a84d6SLisandro Dalcin   PetscFunctionBegin;
62211a84d6SLisandro Dalcin   for (i=n-1; i>=2; i--) {
63211a84d6SLisandro Dalcin     bdf->time[i] = bdf->time[i-1];
64211a84d6SLisandro Dalcin     bdf->work[i] = bdf->work[i-1];
65211a84d6SLisandro Dalcin   }
66211a84d6SLisandro Dalcin   bdf->n       = PetscMin(bdf->n+1,n-1);
67211a84d6SLisandro Dalcin   bdf->time[1] = t;
68211a84d6SLisandro Dalcin   bdf->work[1] = tail;
69211a84d6SLisandro Dalcin   ierr = VecCopy(X,tail);CHKERRQ(ierr);
70211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
71211a84d6SLisandro Dalcin }
72211a84d6SLisandro Dalcin 
73211a84d6SLisandro Dalcin #undef __FUNCT__
74211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_VecDot"
75211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_VecDot(TS ts,PetscInt order,PetscReal t,Vec X,Vec Xdot,PetscReal *shift)
76211a84d6SLisandro Dalcin {
77211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
78211a84d6SLisandro Dalcin   PetscInt       i,n = order+1;
79211a84d6SLisandro Dalcin   PetscReal      time[7];
80211a84d6SLisandro Dalcin   Vec            vecs[7];
81211a84d6SLisandro Dalcin   PetscScalar    alpha[7];
82211a84d6SLisandro Dalcin   PetscErrorCode ierr;
83211a84d6SLisandro Dalcin 
84211a84d6SLisandro Dalcin   PetscFunctionBegin;
85211a84d6SLisandro Dalcin   n = PetscMax(n,2);
86211a84d6SLisandro Dalcin   time[0] = t; for (i=1; i<n; i++) time[i] = bdf->time[i];
87211a84d6SLisandro Dalcin   vecs[0] = X; for (i=1; i<n; i++) vecs[i] = bdf->work[i];
88211a84d6SLisandro Dalcin   LagrangeBasisDers(n,t,time,alpha);
89211a84d6SLisandro Dalcin   ierr = VecZeroEntries(Xdot);CHKERRQ(ierr);
90211a84d6SLisandro Dalcin   ierr = VecMAXPY(Xdot,n,alpha,vecs);CHKERRQ(ierr);
91211a84d6SLisandro Dalcin   if (shift) *shift = PetscRealPart(alpha[0]);
92211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
93211a84d6SLisandro Dalcin }
94211a84d6SLisandro Dalcin 
95211a84d6SLisandro Dalcin #undef __FUNCT__
96211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_VecLTE"
97211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
98211a84d6SLisandro Dalcin {
99211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
100211a84d6SLisandro Dalcin   PetscInt       i,n = order+1;
101211a84d6SLisandro Dalcin   PetscReal      *time = bdf->time;
102211a84d6SLisandro Dalcin   Vec            *vecs = bdf->work;
103211a84d6SLisandro Dalcin   PetscScalar    a[8],b[8],alpha[8];
104211a84d6SLisandro Dalcin   PetscErrorCode ierr;
105211a84d6SLisandro Dalcin 
106211a84d6SLisandro Dalcin   PetscFunctionBegin;
107211a84d6SLisandro Dalcin   LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
108211a84d6SLisandro Dalcin   LagrangeBasisDers(n+1,time[0],time,b);
109211a84d6SLisandro Dalcin   for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
110211a84d6SLisandro Dalcin   ierr = VecZeroEntries(lte);CHKERRQ(ierr);
111211a84d6SLisandro Dalcin   ierr = VecMAXPY(lte,n+1,alpha,vecs);CHKERRQ(ierr);
112211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
113211a84d6SLisandro Dalcin }
114211a84d6SLisandro Dalcin 
115211a84d6SLisandro Dalcin #undef __FUNCT__
116211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Extrapolate"
117211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
118211a84d6SLisandro Dalcin {
119211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
120211a84d6SLisandro Dalcin   PetscInt       n = order+1;
121211a84d6SLisandro Dalcin   PetscReal      *time = bdf->time+1;
122211a84d6SLisandro Dalcin   Vec            *vecs = bdf->work+1;
123211a84d6SLisandro Dalcin   PetscScalar    alpha[7];
124211a84d6SLisandro Dalcin   PetscErrorCode ierr;
125211a84d6SLisandro Dalcin 
126211a84d6SLisandro Dalcin   PetscFunctionBegin;
127211a84d6SLisandro Dalcin   n = PetscMin(n,bdf->n);
128211a84d6SLisandro Dalcin   LagrangeBasisVals(n,t,time,alpha);
129211a84d6SLisandro Dalcin   ierr = VecZeroEntries(X);CHKERRQ(ierr);
130211a84d6SLisandro Dalcin   ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr);
131211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
132211a84d6SLisandro Dalcin }
133211a84d6SLisandro Dalcin 
134211a84d6SLisandro Dalcin #undef __FUNCT__
135211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Interpolate"
136211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
137211a84d6SLisandro Dalcin {
138211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
139211a84d6SLisandro Dalcin   PetscInt       n = order+1;
140211a84d6SLisandro Dalcin   PetscReal      *time = bdf->time;
141211a84d6SLisandro Dalcin   Vec            *vecs = bdf->work;
142211a84d6SLisandro Dalcin   PetscScalar    alpha[7];
143211a84d6SLisandro Dalcin   PetscErrorCode ierr;
144211a84d6SLisandro Dalcin 
145211a84d6SLisandro Dalcin   PetscFunctionBegin;
146211a84d6SLisandro Dalcin   LagrangeBasisVals(n,t,time,alpha);
147211a84d6SLisandro Dalcin   ierr = VecZeroEntries(X);CHKERRQ(ierr);
148211a84d6SLisandro Dalcin   ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr);
149211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
150211a84d6SLisandro Dalcin }
151211a84d6SLisandro Dalcin 
152211a84d6SLisandro Dalcin #undef __FUNCT__
153211a84d6SLisandro Dalcin #define __FUNCT__ "TS_SNESSolve"
154211a84d6SLisandro Dalcin static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
155211a84d6SLisandro Dalcin {
156211a84d6SLisandro Dalcin   PetscInt       nits,lits;
157211a84d6SLisandro Dalcin   PetscErrorCode ierr;
158211a84d6SLisandro Dalcin 
159211a84d6SLisandro Dalcin   PetscFunctionBegin;
160211a84d6SLisandro Dalcin   ierr = SNESSolve(ts->snes,b,x);CHKERRQ(ierr);
161211a84d6SLisandro Dalcin   ierr = SNESGetIterationNumber(ts->snes,&nits);CHKERRQ(ierr);
162211a84d6SLisandro Dalcin   ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr);
163211a84d6SLisandro Dalcin   ts->snes_its += nits; ts->ksp_its += lits;
164211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
165211a84d6SLisandro Dalcin }
166211a84d6SLisandro Dalcin 
167211a84d6SLisandro Dalcin #undef __FUNCT__
168211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Restart"
169211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
170211a84d6SLisandro Dalcin {
171211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
172211a84d6SLisandro Dalcin   PetscErrorCode ierr;
173211a84d6SLisandro Dalcin 
174211a84d6SLisandro Dalcin   PetscFunctionBegin;
175211a84d6SLisandro Dalcin   bdf->k = 1; bdf->n = 0;
176211a84d6SLisandro Dalcin   ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr);
177211a84d6SLisandro Dalcin 
178211a84d6SLisandro Dalcin   bdf->time[0] = ts->ptime + ts->time_step/2;
179211a84d6SLisandro Dalcin   ierr = VecCopy(bdf->work[1],bdf->work[0]);CHKERRQ(ierr);
180211a84d6SLisandro Dalcin   ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr);
181211a84d6SLisandro Dalcin   ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr);
182211a84d6SLisandro Dalcin   ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr);
183211a84d6SLisandro Dalcin   ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);CHKERRQ(ierr);
184211a84d6SLisandro Dalcin   if (!*accept) PetscFunctionReturn(0);
185211a84d6SLisandro Dalcin 
186211a84d6SLisandro Dalcin   bdf->k = PetscMin(2,bdf->order); bdf->n++;
187211a84d6SLisandro Dalcin   ierr = VecCopy(bdf->work[0],bdf->work[2]);CHKERRQ(ierr);
188211a84d6SLisandro Dalcin   bdf->time[2] = bdf->time[0];
189211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
190211a84d6SLisandro Dalcin }
191211a84d6SLisandro Dalcin 
192211a84d6SLisandro Dalcin static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
193211a84d6SLisandro Dalcin 
194211a84d6SLisandro Dalcin #undef __FUNCT__
195211a84d6SLisandro Dalcin #define __FUNCT__ "TSStep_BDF"
196211a84d6SLisandro Dalcin static PetscErrorCode TSStep_BDF(TS ts)
197211a84d6SLisandro Dalcin {
198211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
199211a84d6SLisandro Dalcin   PetscInt       rejections = 0;
200211a84d6SLisandro Dalcin   PetscBool      stageok,accept = PETSC_TRUE;
201211a84d6SLisandro Dalcin   PetscReal      next_time_step = ts->time_step;
202211a84d6SLisandro Dalcin   PetscErrorCode ierr;
203211a84d6SLisandro Dalcin 
204211a84d6SLisandro Dalcin   PetscFunctionBegin;
205211a84d6SLisandro Dalcin   ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
206211a84d6SLisandro Dalcin 
207211a84d6SLisandro Dalcin   if (!ts->steprollback && !ts->steprestart) {
208211a84d6SLisandro Dalcin     bdf->k = PetscMin(bdf->k+1,bdf->order);
209211a84d6SLisandro Dalcin     ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr);
210211a84d6SLisandro Dalcin   }
211211a84d6SLisandro Dalcin 
212211a84d6SLisandro Dalcin   bdf->status = TS_STEP_INCOMPLETE;
213211a84d6SLisandro Dalcin   while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
214211a84d6SLisandro Dalcin 
215211a84d6SLisandro Dalcin     if (ts->steprestart) {
216211a84d6SLisandro Dalcin       ierr = TSBDF_Restart(ts,&stageok);CHKERRQ(ierr);
217211a84d6SLisandro Dalcin       if (!stageok) goto reject_step;
218211a84d6SLisandro Dalcin     }
219211a84d6SLisandro Dalcin 
220211a84d6SLisandro Dalcin     bdf->time[0] = ts->ptime + ts->time_step;
221211a84d6SLisandro Dalcin     ierr = TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);CHKERRQ(ierr);
222211a84d6SLisandro Dalcin     ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr);
223211a84d6SLisandro Dalcin     ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr);
224211a84d6SLisandro Dalcin     ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr);
225211a84d6SLisandro Dalcin     ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);CHKERRQ(ierr);
226211a84d6SLisandro Dalcin     if (!stageok) goto reject_step;
227211a84d6SLisandro Dalcin 
228211a84d6SLisandro Dalcin     bdf->status = TS_STEP_PENDING;
229211a84d6SLisandro Dalcin     ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr);
230211a84d6SLisandro Dalcin     ierr = TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);CHKERRQ(ierr);
231211a84d6SLisandro Dalcin     ierr = TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr);
232211a84d6SLisandro Dalcin     bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
233211a84d6SLisandro Dalcin     if (!accept) { ts->time_step = next_time_step; goto reject_step; }
234211a84d6SLisandro Dalcin 
235211a84d6SLisandro Dalcin     ierr = VecCopy(bdf->work[0],ts->vec_sol);CHKERRQ(ierr);
236211a84d6SLisandro Dalcin     ts->ptime += ts->time_step;
237211a84d6SLisandro Dalcin     ts->time_step = next_time_step;
238211a84d6SLisandro Dalcin     break;
239211a84d6SLisandro Dalcin 
240211a84d6SLisandro Dalcin   reject_step:
241211a84d6SLisandro Dalcin     ts->reject++; accept = PETSC_FALSE;
242211a84d6SLisandro Dalcin     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
243211a84d6SLisandro Dalcin       ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr);
244211a84d6SLisandro Dalcin       ts->reason = TS_DIVERGED_STEP_REJECTED;
245211a84d6SLisandro Dalcin     }
246211a84d6SLisandro Dalcin   }
247211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
248211a84d6SLisandro Dalcin }
249211a84d6SLisandro Dalcin 
250211a84d6SLisandro Dalcin #undef __FUNCT__
251211a84d6SLisandro Dalcin #define __FUNCT__ "TSInterpolate_BDF"
252211a84d6SLisandro Dalcin static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
253211a84d6SLisandro Dalcin {
254211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
255211a84d6SLisandro Dalcin   PetscErrorCode ierr;
256211a84d6SLisandro Dalcin 
257211a84d6SLisandro Dalcin   PetscFunctionBegin;
258211a84d6SLisandro Dalcin   ierr = TSBDF_Interpolate(ts,bdf->k,t,X);CHKERRQ(ierr);
259211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
260211a84d6SLisandro Dalcin }
261211a84d6SLisandro Dalcin 
262211a84d6SLisandro Dalcin #undef __FUNCT__
263211a84d6SLisandro Dalcin #define __FUNCT__ "TSEvaluateWLTE_BDF"
264211a84d6SLisandro Dalcin static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
265211a84d6SLisandro Dalcin {
266211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
267211a84d6SLisandro Dalcin   PetscInt       k = bdf->k;
268*7453f775SEmil Constantinescu   PetscReal      wltea,wlter;
269211a84d6SLisandro Dalcin   Vec            X = bdf->work[0], Y = bdf->vec_lte;
270211a84d6SLisandro Dalcin   PetscErrorCode ierr;
271211a84d6SLisandro Dalcin 
272211a84d6SLisandro Dalcin   PetscFunctionBegin;
273211a84d6SLisandro Dalcin   k = PetscMin(k,bdf->n-1);
274211a84d6SLisandro Dalcin   ierr = TSBDF_VecLTE(ts,k,Y);CHKERRQ(ierr);
275211a84d6SLisandro Dalcin   ierr = VecAXPY(Y,1,X);CHKERRQ(ierr);
276*7453f775SEmil Constantinescu   ierr = TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);CHKERRQ(ierr);
277211a84d6SLisandro Dalcin   if (order) *order = k + 1;
278211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
279211a84d6SLisandro Dalcin }
280211a84d6SLisandro Dalcin 
281211a84d6SLisandro Dalcin #undef __FUNCT__
282211a84d6SLisandro Dalcin #define __FUNCT__ "TSRollBack_BDF"
283211a84d6SLisandro Dalcin static PetscErrorCode TSRollBack_BDF(TS ts)
284211a84d6SLisandro Dalcin {
285211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
286211a84d6SLisandro Dalcin   PetscErrorCode ierr;
287211a84d6SLisandro Dalcin 
288211a84d6SLisandro Dalcin   PetscFunctionBegin;
289211a84d6SLisandro Dalcin   ierr = VecCopy(bdf->work[1],ts->vec_sol);CHKERRQ(ierr);
290211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
291211a84d6SLisandro Dalcin }
292211a84d6SLisandro Dalcin 
293211a84d6SLisandro Dalcin #undef __FUNCT__
294211a84d6SLisandro Dalcin #define __FUNCT__ "SNESTSFormFunction_BDF"
295211a84d6SLisandro Dalcin static PetscErrorCode SNESTSFormFunction_BDF(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
296211a84d6SLisandro Dalcin {
297211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
298211a84d6SLisandro Dalcin   PetscInt       k = bdf->k;
299211a84d6SLisandro Dalcin   PetscReal      t = bdf->time[0];
300211a84d6SLisandro Dalcin   Vec            V = bdf->vec_dot;
301211a84d6SLisandro Dalcin   PetscErrorCode ierr;
302211a84d6SLisandro Dalcin 
303211a84d6SLisandro Dalcin   PetscFunctionBegin;
304211a84d6SLisandro Dalcin   ierr = TSBDF_VecDot(ts,k,t,X,V,&bdf->shift);CHKERRQ(ierr);
305211a84d6SLisandro Dalcin   /* F = Function(t,X,V) */
306211a84d6SLisandro Dalcin   ierr = TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);CHKERRQ(ierr);
307211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
308211a84d6SLisandro Dalcin }
309211a84d6SLisandro Dalcin 
310211a84d6SLisandro Dalcin #undef __FUNCT__
311211a84d6SLisandro Dalcin #define __FUNCT__ "SNESTSFormJacobian_BDF"
312211a84d6SLisandro Dalcin static PetscErrorCode SNESTSFormJacobian_BDF(PETSC_UNUSED SNES snes,
313211a84d6SLisandro Dalcin                                              PETSC_UNUSED Vec X,
314211a84d6SLisandro Dalcin                                              Mat J,Mat P,
315211a84d6SLisandro Dalcin                                              TS ts)
316211a84d6SLisandro Dalcin {
317211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
318211a84d6SLisandro Dalcin   PetscReal      t = bdf->time[0];
319211a84d6SLisandro Dalcin   Vec            V = bdf->vec_dot;
320211a84d6SLisandro Dalcin   PetscReal      dVdX = bdf->shift;
321211a84d6SLisandro Dalcin   PetscErrorCode ierr;
322211a84d6SLisandro Dalcin 
323211a84d6SLisandro Dalcin   PetscFunctionBegin;
324211a84d6SLisandro Dalcin   /* J,P = Jacobian(t,X,V) */
325211a84d6SLisandro Dalcin   ierr = TSComputeIJacobian(ts,t,X,V,dVdX,J,P,PETSC_FALSE);CHKERRQ(ierr);
326211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
327211a84d6SLisandro Dalcin }
328211a84d6SLisandro Dalcin 
329211a84d6SLisandro Dalcin #undef __FUNCT__
330211a84d6SLisandro Dalcin #define __FUNCT__ "TSReset_BDF"
331211a84d6SLisandro Dalcin static PetscErrorCode TSReset_BDF(TS ts)
332211a84d6SLisandro Dalcin {
333211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
334211a84d6SLisandro Dalcin   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);
335211a84d6SLisandro Dalcin   PetscErrorCode ierr;
336211a84d6SLisandro Dalcin 
337211a84d6SLisandro Dalcin   PetscFunctionBegin;
338211a84d6SLisandro Dalcin   for (i=0; i<n; i++) {ierr = VecDestroy(&bdf->work[i]);CHKERRQ(ierr);}
339211a84d6SLisandro Dalcin   ierr = VecDestroy(&bdf->vec_dot);CHKERRQ(ierr);
340211a84d6SLisandro Dalcin   ierr = VecDestroy(&bdf->vec_lte);CHKERRQ(ierr);
341211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
342211a84d6SLisandro Dalcin }
343211a84d6SLisandro Dalcin 
344211a84d6SLisandro Dalcin #undef __FUNCT__
345211a84d6SLisandro Dalcin #define __FUNCT__ "TSDestroy_BDF"
346211a84d6SLisandro Dalcin static PetscErrorCode TSDestroy_BDF(TS ts)
347211a84d6SLisandro Dalcin {
348211a84d6SLisandro Dalcin   PetscErrorCode ierr;
349211a84d6SLisandro Dalcin 
350211a84d6SLisandro Dalcin   PetscFunctionBegin;
351211a84d6SLisandro Dalcin   ierr = TSReset_BDF(ts);CHKERRQ(ierr);
352211a84d6SLisandro Dalcin   ierr = PetscFree(ts->data);CHKERRQ(ierr);
353211a84d6SLisandro Dalcin   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);CHKERRQ(ierr);
354211a84d6SLisandro Dalcin   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);CHKERRQ(ierr);
355211a84d6SLisandro Dalcin   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",NULL);CHKERRQ(ierr);
356211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
357211a84d6SLisandro Dalcin }
358211a84d6SLisandro Dalcin 
359211a84d6SLisandro Dalcin #undef __FUNCT__
360211a84d6SLisandro Dalcin #define __FUNCT__ "TSSetUp_BDF"
361211a84d6SLisandro Dalcin static PetscErrorCode TSSetUp_BDF(TS ts)
362211a84d6SLisandro Dalcin {
363211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
364211a84d6SLisandro Dalcin   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);
365211a84d6SLisandro Dalcin   PetscErrorCode ierr;
366211a84d6SLisandro Dalcin 
367211a84d6SLisandro Dalcin   PetscFunctionBegin;
368211a84d6SLisandro Dalcin   bdf->k = bdf->n = 0;
369211a84d6SLisandro Dalcin   for (i=0; i<n; i++) {ierr = VecDuplicate(ts->vec_sol,&bdf->work[i]);CHKERRQ(ierr);}
370211a84d6SLisandro Dalcin   ierr = VecDuplicate(ts->vec_sol,&bdf->vec_dot);CHKERRQ(ierr);
371211a84d6SLisandro Dalcin   ierr = VecDuplicate(ts->vec_sol,&bdf->vec_lte);CHKERRQ(ierr);
372211a84d6SLisandro Dalcin 
373211a84d6SLisandro Dalcin   ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr);
374211a84d6SLisandro Dalcin   ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr);
375211a84d6SLisandro Dalcin   if (!bdf->adapt) {
376211a84d6SLisandro Dalcin     ierr = TSAdaptSetType(ts->adapt,TSADAPTNONE);CHKERRQ(ierr);
377211a84d6SLisandro Dalcin   } else {
378211a84d6SLisandro Dalcin     PetscReal low,high;
379211a84d6SLisandro Dalcin     ierr = TSAdaptBasicGetClip(ts->adapt,&low,&high);CHKERRQ(ierr);
380211a84d6SLisandro Dalcin     high = PetscMin(high,2.0);
381211a84d6SLisandro Dalcin     ierr = TSAdaptBasicSetClip(ts->adapt,low,high);CHKERRQ(ierr);
382211a84d6SLisandro Dalcin     if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
383211a84d6SLisandro Dalcin       ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
384211a84d6SLisandro Dalcin   }
385211a84d6SLisandro Dalcin 
386211a84d6SLisandro Dalcin   ierr = TSGetSNES(ts,&ts->snes);CHKERRQ(ierr);
387211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
388211a84d6SLisandro Dalcin }
389211a84d6SLisandro Dalcin 
390211a84d6SLisandro Dalcin #undef __FUNCT__
391211a84d6SLisandro Dalcin #define __FUNCT__ "TSSetFromOptions_BDF"
392211a84d6SLisandro Dalcin static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
393211a84d6SLisandro Dalcin {
394211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
395211a84d6SLisandro Dalcin   PetscErrorCode ierr;
396211a84d6SLisandro Dalcin 
397211a84d6SLisandro Dalcin   PetscFunctionBegin;
398211a84d6SLisandro Dalcin   ierr = PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");CHKERRQ(ierr);
399211a84d6SLisandro Dalcin   {
400211a84d6SLisandro Dalcin     PetscBool flg;
401211a84d6SLisandro Dalcin     PetscInt  order = bdf->order;
402211a84d6SLisandro Dalcin     PetscBool adapt = bdf->adapt;
403211a84d6SLisandro Dalcin     ierr = PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);CHKERRQ(ierr);
404211a84d6SLisandro Dalcin     if (flg) {ierr = TSBDFSetOrder(ts,order);CHKERRQ(ierr);}
405211a84d6SLisandro Dalcin     ierr = PetscOptionsBool("-ts_bdf_adapt","Use time-step adaptivity with the BDF method","TSBDFUseAdapt",adapt,&adapt,&flg);CHKERRQ(ierr);
406211a84d6SLisandro Dalcin     if (flg) {ierr = TSBDFUseAdapt(ts,adapt);CHKERRQ(ierr);}
407211a84d6SLisandro Dalcin   }
408211a84d6SLisandro Dalcin   ierr = PetscOptionsTail();CHKERRQ(ierr);
409211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
410211a84d6SLisandro Dalcin }
411211a84d6SLisandro Dalcin 
412211a84d6SLisandro Dalcin #undef __FUNCT__
413211a84d6SLisandro Dalcin #define __FUNCT__ "TSView_BDF"
414211a84d6SLisandro Dalcin static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
415211a84d6SLisandro Dalcin {
416211a84d6SLisandro Dalcin   TS_BDF         *bdf = (TS_BDF*)ts->data;
417211a84d6SLisandro Dalcin   PetscBool      iascii;
418211a84d6SLisandro Dalcin   PetscErrorCode ierr;
419211a84d6SLisandro Dalcin 
420211a84d6SLisandro Dalcin   PetscFunctionBegin;
421211a84d6SLisandro Dalcin   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
422211a84d6SLisandro Dalcin   if (iascii)    {
423211a84d6SLisandro Dalcin     ierr = PetscViewerASCIIPrintf(viewer,"  Order=%D\n",bdf->order);CHKERRQ(ierr);
424211a84d6SLisandro Dalcin   }
425211a84d6SLisandro Dalcin   if (ts->adapt) {ierr = TSAdaptView(ts->adapt,viewer);CHKERRQ(ierr);}
426211a84d6SLisandro Dalcin   if (ts->snes)  {ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr);}
427211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
428211a84d6SLisandro Dalcin }
429211a84d6SLisandro Dalcin 
430211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */
431211a84d6SLisandro Dalcin 
432211a84d6SLisandro Dalcin #undef __FUNCT__
433211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFSetOrder_BDF"
434211a84d6SLisandro Dalcin static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
435211a84d6SLisandro Dalcin {
436211a84d6SLisandro Dalcin   TS_BDF *bdf = (TS_BDF*)ts->data;
437211a84d6SLisandro Dalcin 
438211a84d6SLisandro Dalcin   PetscFunctionBegin;
439211a84d6SLisandro Dalcin   if (order == bdf->order) PetscFunctionReturn(0);
440211a84d6SLisandro Dalcin   if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
441211a84d6SLisandro Dalcin   bdf->order = order;
442211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
443211a84d6SLisandro Dalcin }
444211a84d6SLisandro Dalcin 
445211a84d6SLisandro Dalcin #undef __FUNCT__
446211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFGetOrder_BDF"
447211a84d6SLisandro Dalcin static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
448211a84d6SLisandro Dalcin {
449211a84d6SLisandro Dalcin   TS_BDF *bdf = (TS_BDF*)ts->data;
450211a84d6SLisandro Dalcin 
451211a84d6SLisandro Dalcin   PetscFunctionBegin;
452211a84d6SLisandro Dalcin   *order = bdf->order;
453211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
454211a84d6SLisandro Dalcin }
455211a84d6SLisandro Dalcin 
456211a84d6SLisandro Dalcin #undef __FUNCT__
457211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFUseAdapt_BDF"
458211a84d6SLisandro Dalcin static PetscErrorCode TSBDFUseAdapt_BDF(TS ts,PetscBool use)
459211a84d6SLisandro Dalcin {
460211a84d6SLisandro Dalcin   TS_BDF *bdf = (TS_BDF*)ts->data;
461211a84d6SLisandro Dalcin 
462211a84d6SLisandro Dalcin   PetscFunctionBegin;
463211a84d6SLisandro Dalcin   if (use == bdf->adapt) PetscFunctionReturn(0);
464211a84d6SLisandro Dalcin   if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()");
465211a84d6SLisandro Dalcin   bdf->adapt = use;
466211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
467211a84d6SLisandro Dalcin }
468211a84d6SLisandro Dalcin 
469211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */
470211a84d6SLisandro Dalcin 
471211a84d6SLisandro Dalcin /*MC
472211a84d6SLisandro Dalcin       TSBDF - DAE solver using BDF methods
473211a84d6SLisandro Dalcin 
474211a84d6SLisandro Dalcin   Level: beginner
475211a84d6SLisandro Dalcin 
476211a84d6SLisandro Dalcin .seealso:  TS, TSCreate(), TSSetType()
477211a84d6SLisandro Dalcin M*/
478211a84d6SLisandro Dalcin #undef __FUNCT__
479211a84d6SLisandro Dalcin #define __FUNCT__ "TSCreate_BDF"
480211a84d6SLisandro Dalcin PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
481211a84d6SLisandro Dalcin {
482211a84d6SLisandro Dalcin   TS_BDF         *bdf;
483211a84d6SLisandro Dalcin   PetscErrorCode ierr;
484211a84d6SLisandro Dalcin 
485211a84d6SLisandro Dalcin   PetscFunctionBegin;
486211a84d6SLisandro Dalcin   ts->ops->reset          = TSReset_BDF;
487211a84d6SLisandro Dalcin   ts->ops->destroy        = TSDestroy_BDF;
488211a84d6SLisandro Dalcin   ts->ops->view           = TSView_BDF;
489211a84d6SLisandro Dalcin   ts->ops->setup          = TSSetUp_BDF;
490211a84d6SLisandro Dalcin   ts->ops->setfromoptions = TSSetFromOptions_BDF;
491211a84d6SLisandro Dalcin   ts->ops->step           = TSStep_BDF;
492211a84d6SLisandro Dalcin   ts->ops->evaluatewlte   = TSEvaluateWLTE_BDF;
493211a84d6SLisandro Dalcin   ts->ops->rollback       = TSRollBack_BDF;
494211a84d6SLisandro Dalcin   ts->ops->interpolate    = TSInterpolate_BDF;
495211a84d6SLisandro Dalcin   ts->ops->snesfunction   = SNESTSFormFunction_BDF;
496211a84d6SLisandro Dalcin   ts->ops->snesjacobian   = SNESTSFormJacobian_BDF;
497211a84d6SLisandro Dalcin 
498211a84d6SLisandro Dalcin   ierr = PetscNewLog(ts,&bdf);CHKERRQ(ierr);
499211a84d6SLisandro Dalcin   ts->data = (void*)bdf;
500211a84d6SLisandro Dalcin 
501211a84d6SLisandro Dalcin   bdf->order  = 2;
502211a84d6SLisandro Dalcin   bdf->adapt  = PETSC_FALSE;
503211a84d6SLisandro Dalcin   bdf->status = TS_STEP_COMPLETE;
504211a84d6SLisandro Dalcin 
505211a84d6SLisandro Dalcin   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);CHKERRQ(ierr);
506211a84d6SLisandro Dalcin   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);CHKERRQ(ierr);
507211a84d6SLisandro Dalcin   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",TSBDFUseAdapt_BDF);CHKERRQ(ierr);
508211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
509211a84d6SLisandro Dalcin }
510211a84d6SLisandro Dalcin 
511211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */
512211a84d6SLisandro Dalcin 
513211a84d6SLisandro Dalcin #undef __FUNCT__
514211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFSetOrder"
515211a84d6SLisandro Dalcin /*@
516211a84d6SLisandro Dalcin   TSBDFSetOrder - Set the order of the BDF method
517211a84d6SLisandro Dalcin 
518211a84d6SLisandro Dalcin   Logically Collective on TS
519211a84d6SLisandro Dalcin 
520211a84d6SLisandro Dalcin   Input Parameter:
521211a84d6SLisandro Dalcin +  ts - timestepping context
522211a84d6SLisandro Dalcin -  order - order of the method
523211a84d6SLisandro Dalcin 
524211a84d6SLisandro Dalcin   Options Database:
525211a84d6SLisandro Dalcin .  -ts_bdf_order <order>
526211a84d6SLisandro Dalcin 
527211a84d6SLisandro Dalcin   Level: intermediate
528211a84d6SLisandro Dalcin 
529211a84d6SLisandro Dalcin @*/
530211a84d6SLisandro Dalcin PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
531211a84d6SLisandro Dalcin {
532211a84d6SLisandro Dalcin   PetscErrorCode ierr;
533211a84d6SLisandro Dalcin 
534211a84d6SLisandro Dalcin   PetscFunctionBegin;
535211a84d6SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
536211a84d6SLisandro Dalcin   PetscValidLogicalCollectiveInt(ts,order,2);
537211a84d6SLisandro Dalcin   ierr = PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));CHKERRQ(ierr);
538211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
539211a84d6SLisandro Dalcin }
540211a84d6SLisandro Dalcin 
541211a84d6SLisandro Dalcin #undef __FUNCT__
542211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFGetOrder"
543211a84d6SLisandro Dalcin /*@
544211a84d6SLisandro Dalcin   TSBDFGetOrder - Get the order of the BDF method
545211a84d6SLisandro Dalcin 
546211a84d6SLisandro Dalcin   Not Collective
547211a84d6SLisandro Dalcin 
548211a84d6SLisandro Dalcin   Input Parameter:
549211a84d6SLisandro Dalcin .  ts - timestepping context
550211a84d6SLisandro Dalcin 
551211a84d6SLisandro Dalcin   Output Parameter:
552211a84d6SLisandro Dalcin .  order - order of the method
553211a84d6SLisandro Dalcin 
554211a84d6SLisandro Dalcin   Level: intermediate
555211a84d6SLisandro Dalcin 
556211a84d6SLisandro Dalcin @*/
557211a84d6SLisandro Dalcin PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
558211a84d6SLisandro Dalcin {
559211a84d6SLisandro Dalcin   PetscErrorCode ierr;
560211a84d6SLisandro Dalcin 
561211a84d6SLisandro Dalcin   PetscFunctionBegin;
562211a84d6SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
563211a84d6SLisandro Dalcin   PetscValidIntPointer(order,2);
564211a84d6SLisandro Dalcin   ierr = PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));CHKERRQ(ierr);
565211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
566211a84d6SLisandro Dalcin }
567211a84d6SLisandro Dalcin 
568211a84d6SLisandro Dalcin #undef __FUNCT__
569211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFUseAdapt"
570211a84d6SLisandro Dalcin /*@
571211a84d6SLisandro Dalcin   TSBDFUseAdapt - Use time-step adaptivity with the BDF method
572211a84d6SLisandro Dalcin 
573211a84d6SLisandro Dalcin   Logically Collective on TS
574211a84d6SLisandro Dalcin 
575211a84d6SLisandro Dalcin   Input Parameter:
576211a84d6SLisandro Dalcin +  ts - timestepping context
577211a84d6SLisandro Dalcin -  use - flag to use adaptivity
578211a84d6SLisandro Dalcin 
579211a84d6SLisandro Dalcin   Options Database:
580211a84d6SLisandro Dalcin .  -ts_bdf_adapt
581211a84d6SLisandro Dalcin 
582211a84d6SLisandro Dalcin   Level: intermediate
583211a84d6SLisandro Dalcin 
584211a84d6SLisandro Dalcin .seealso: TSAdapt
585211a84d6SLisandro Dalcin @*/
586211a84d6SLisandro Dalcin PetscErrorCode TSBDFUseAdapt(TS ts,PetscBool use)
587211a84d6SLisandro Dalcin {
588211a84d6SLisandro Dalcin   PetscErrorCode ierr;
589211a84d6SLisandro Dalcin 
590211a84d6SLisandro Dalcin   PetscFunctionBegin;
591211a84d6SLisandro Dalcin   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
592211a84d6SLisandro Dalcin   PetscValidLogicalCollectiveBool(ts,use,2);
593211a84d6SLisandro Dalcin   ierr = PetscTryMethod(ts,"TSBDFUseAdapt_C",(TS,PetscBool),(ts,use));CHKERRQ(ierr);
594211a84d6SLisandro Dalcin   PetscFunctionReturn(0);
595211a84d6SLisandro Dalcin }
596