xref: /petsc/src/ts/impls/implicit/discgrad/tsdiscgrad.c (revision 4e278199b78715991f5c71ebbd945c1489263e6c)
1 /*
2   Code for timestepping with discrete gradient integrators
3 */
4 #include <petsc/private/tsimpl.h>                /*I   "petscts.h"   I*/
5 #include <petscdm.h>
6 
7 PetscBool DGCite = PETSC_FALSE;
8 const char DGCitation[] = "@article{Gonzalez1996,\n"
9                           "  title   = {Time integration and discrete Hamiltonian systems},\n"
10                           "  author  = {Oscar Gonzalez},\n"
11                           "  journal = {Journal of Nonlinear Science},\n"
12                           "  volume  = {6},\n"
13                           "  pages   = {449--467},\n"
14                           "  doi     = {10.1007/978-1-4612-1246-1_10},\n"
15                           "  year    = {1996}\n}\n";
16 
17 typedef struct {
18   PetscReal stage_time;
19   Vec       X0, X, Xdot;
20   PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *);
21   PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *);
22   PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *);
23 } TS_DiscGrad;
24 
25 static PetscErrorCode TSDiscGradGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
26 {
27   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;
28   PetscErrorCode ierr;
29 
30   PetscFunctionBegin;
31   if (X0) {
32     if (dm && dm != ts->dm) {ierr = DMGetNamedGlobalVector(dm, "TSDiscGrad_X0", X0);CHKERRQ(ierr);}
33     else                    {*X0  = ts->vec_sol;}
34   }
35   if (Xdot) {
36     if (dm && dm != ts->dm) {ierr = DMGetNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);CHKERRQ(ierr);}
37     else                    {*Xdot = dg->Xdot;}
38   }
39   PetscFunctionReturn(0);
40 }
41 
42 static PetscErrorCode TSDiscGradRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
43 {
44   PetscErrorCode ierr;
45 
46   PetscFunctionBegin;
47   if (X0) {
48     if (dm && dm != ts->dm) {ierr = DMRestoreNamedGlobalVector(dm, "TSDiscGrad_X0", X0);CHKERRQ(ierr);}
49   }
50   if (Xdot) {
51     if (dm && dm != ts->dm) {ierr = DMRestoreNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);CHKERRQ(ierr);}
52   }
53   PetscFunctionReturn(0);
54 }
55 
56 static PetscErrorCode DMCoarsenHook_TSDiscGrad(DM fine, DM coarse, void *ctx)
57 {
58   PetscFunctionBegin;
59   PetscFunctionReturn(0);
60 }
61 
62 static PetscErrorCode DMRestrictHook_TSDiscGrad(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
63 {
64   TS             ts = (TS) ctx;
65   Vec            X0, Xdot, X0_c, Xdot_c;
66   PetscErrorCode ierr;
67 
68   PetscFunctionBegin;
69   ierr = TSDiscGradGetX0AndXdot(ts, fine, &X0, &Xdot);CHKERRQ(ierr);
70   ierr = TSDiscGradGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c);CHKERRQ(ierr);
71   ierr = MatRestrict(restrct, X0, X0_c);CHKERRQ(ierr);
72   ierr = MatRestrict(restrct, Xdot, Xdot_c);CHKERRQ(ierr);
73   ierr = VecPointwiseMult(X0_c, rscale, X0_c);CHKERRQ(ierr);
74   ierr = VecPointwiseMult(Xdot_c, rscale, Xdot_c);CHKERRQ(ierr);
75   ierr = TSDiscGradRestoreX0AndXdot(ts, fine, &X0, &Xdot);CHKERRQ(ierr);
76   ierr = TSDiscGradRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c);CHKERRQ(ierr);
77   PetscFunctionReturn(0);
78 }
79 
80 static PetscErrorCode DMSubDomainHook_TSDiscGrad(DM dm, DM subdm, void *ctx)
81 {
82   PetscFunctionBegin;
83   PetscFunctionReturn(0);
84 }
85 
86 static PetscErrorCode DMSubDomainRestrictHook_TSDiscGrad(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
87 {
88   TS             ts = (TS) ctx;
89   Vec            X0, Xdot, X0_sub, Xdot_sub;
90   PetscErrorCode ierr;
91 
92   PetscFunctionBegin;
93   ierr = TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr);
94   ierr = TSDiscGradGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);CHKERRQ(ierr);
95 
96   ierr = VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
97   ierr = VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
98 
99   ierr = VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
100   ierr = VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
101 
102   ierr = TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr);
103   ierr = TSDiscGradRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);CHKERRQ(ierr);
104   PetscFunctionReturn(0);
105 }
106 
107 static PetscErrorCode TSSetUp_DiscGrad(TS ts)
108 {
109   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;
110   DM             dm;
111   PetscErrorCode ierr;
112 
113   PetscFunctionBegin;
114   if (!dg->X)    {ierr = VecDuplicate(ts->vec_sol, &dg->X);CHKERRQ(ierr);}
115   if (!dg->X0)   {ierr = VecDuplicate(ts->vec_sol, &dg->X0);CHKERRQ(ierr);}
116   if (!dg->Xdot) {ierr = VecDuplicate(ts->vec_sol, &dg->Xdot);CHKERRQ(ierr);}
117 
118   ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
119   ierr = DMCoarsenHookAdd(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr);
120   ierr = DMSubDomainHookAdd(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr);
121   PetscFunctionReturn(0);
122 }
123 
124 static PetscErrorCode TSSetFromOptions_DiscGrad(PetscOptionItems *PetscOptionsObject, TS ts)
125 {
126   PetscErrorCode ierr;
127 
128   PetscFunctionBegin;
129   ierr = PetscOptionsHead(PetscOptionsObject, "Discrete Gradients ODE solver options");CHKERRQ(ierr);
130   {
131     //ierr = PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSDiscGradSetTheta",th->Theta,&th->Theta,NULL);CHKERRQ(ierr);
132   }
133   ierr = PetscOptionsTail();CHKERRQ(ierr);
134   PetscFunctionReturn(0);
135 }
136 
137 static PetscErrorCode TSView_DiscGrad(TS ts,PetscViewer viewer)
138 {
139   PetscBool      iascii;
140   PetscErrorCode ierr;
141 
142   PetscFunctionBegin;
143   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
144   if (iascii) {
145     ierr = PetscViewerASCIIPrintf(viewer,"  Discrete Gradients\n");CHKERRQ(ierr);
146   }
147   PetscFunctionReturn(0);
148 }
149 
150 static PetscErrorCode TSReset_DiscGrad(TS ts)
151 {
152   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;
153   PetscErrorCode ierr;
154 
155   PetscFunctionBegin;
156   ierr = VecDestroy(&dg->X);CHKERRQ(ierr);
157   ierr = VecDestroy(&dg->X0);CHKERRQ(ierr);
158   ierr = VecDestroy(&dg->Xdot);CHKERRQ(ierr);
159   PetscFunctionReturn(0);
160 }
161 
162 static PetscErrorCode TSDestroy_DiscGrad(TS ts)
163 {
164   DM             dm;
165   PetscErrorCode ierr;
166 
167   PetscFunctionBegin;
168   ierr = TSReset_DiscGrad(ts);CHKERRQ(ierr);
169   ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
170   if (dm) {
171     ierr = DMCoarsenHookRemove(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr);
172     ierr = DMSubDomainHookRemove(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr);
173   }
174   ierr = PetscFree(ts->data);CHKERRQ(ierr);
175   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",NULL);CHKERRQ(ierr);
176   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",NULL);CHKERRQ(ierr);
177   PetscFunctionReturn(0);
178 }
179 
180 static PetscErrorCode TSInterpolate_DiscGrad(TS ts, PetscReal t, Vec X)
181 {
182   TS_DiscGrad   *dg = (TS_DiscGrad*)ts->data;
183   PetscReal      dt = t - ts->ptime;
184   PetscErrorCode ierr;
185 
186   PetscFunctionBegin;
187   ierr = VecCopy(ts->vec_sol, dg->X);CHKERRQ(ierr);
188   ierr = VecWAXPY(X, dt, dg->Xdot, dg->X);CHKERRQ(ierr);
189   PetscFunctionReturn(0);
190 }
191 
192 static PetscErrorCode TSDiscGrad_SNESSolve(TS ts, Vec b, Vec x)
193 {
194   SNES           snes;
195   PetscInt       nits, lits;
196   PetscErrorCode ierr;
197 
198   PetscFunctionBegin;
199   ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr);
200   ierr = SNESSolve(snes, b, x);CHKERRQ(ierr);
201   ierr = SNESGetIterationNumber(snes, &nits);CHKERRQ(ierr);
202   ierr = SNESGetLinearSolveIterations(snes, &lits);CHKERRQ(ierr);
203   ts->snes_its += nits;
204   ts->ksp_its  += lits;
205   PetscFunctionReturn(0);
206 }
207 
208 static PetscErrorCode TSStep_DiscGrad(TS ts)
209 {
210   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;
211   TSAdapt        adapt;
212   TSStepStatus   status          = TS_STEP_INCOMPLETE;
213   PetscInt       rejections      = 0;
214   PetscBool      stageok, accept = PETSC_TRUE;
215   PetscReal      next_time_step  = ts->time_step;
216   PetscErrorCode ierr;
217 
218   PetscFunctionBegin;
219   ierr = TSGetAdapt(ts, &adapt);CHKERRQ(ierr);
220   if (!ts->steprollback) {ierr = VecCopy(ts->vec_sol, dg->X0);CHKERRQ(ierr);}
221 
222   while (!ts->reason && status != TS_STEP_COMPLETE) {
223     PetscReal shift = 1/(0.5*ts->time_step);
224 
225     dg->stage_time = ts->ptime + 0.5*ts->time_step;
226 
227     ierr = VecCopy(dg->X0, dg->X);CHKERRQ(ierr);
228     ierr = TSPreStage(ts, dg->stage_time);CHKERRQ(ierr);
229     ierr = TSDiscGrad_SNESSolve(ts, NULL, dg->X);CHKERRQ(ierr);
230     ierr = TSPostStage(ts, dg->stage_time, 0, &dg->X);CHKERRQ(ierr);
231     ierr = TSAdaptCheckStage(adapt, ts, dg->stage_time, dg->X, &stageok);CHKERRQ(ierr);
232     if (!stageok) goto reject_step;
233 
234     status = TS_STEP_PENDING;
235     ierr = VecAXPBYPCZ(dg->Xdot, -shift, shift, 0, dg->X0, dg->X);CHKERRQ(ierr);
236     ierr = VecAXPY(ts->vec_sol, ts->time_step, dg->Xdot);CHKERRQ(ierr);
237     ierr = TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept);CHKERRQ(ierr);
238     status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
239     if (!accept) {
240       ierr = VecCopy(dg->X0, ts->vec_sol);CHKERRQ(ierr);
241       ts->time_step = next_time_step;
242       goto reject_step;
243     }
244     ts->ptime    += ts->time_step;
245     ts->time_step = next_time_step;
246     break;
247 
248   reject_step:
249     ts->reject++; accept = PETSC_FALSE;
250     if (!ts->reason && ts->max_reject >= 0 && ++rejections > ts->max_reject) {
251       ts->reason = TS_DIVERGED_STEP_REJECTED;
252       ierr = PetscInfo2(ts, "Step=%D, step rejections %D greater than current TS allowed, stopping solve\n", ts->steps, rejections);CHKERRQ(ierr);
253     }
254   }
255   PetscFunctionReturn(0);
256 }
257 
258 static PetscErrorCode TSGetStages_DiscGrad(TS ts, PetscInt *ns, Vec **Y)
259 {
260   TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
261 
262   PetscFunctionBegin;
263   if (ns) *ns = 1;
264   if (Y)  *Y  = &(dg->X);
265   PetscFunctionReturn(0);
266 }
267 
268 /*
269   This defines the nonlinear equation that is to be solved with SNES
270     G(U) = F[t0 + 0.5*dt, U, (U-U0)/dt] = 0
271 */
272 static PetscErrorCode SNESTSFormFunction_DiscGrad(SNES snes, Vec x, Vec y, TS ts)
273 {
274   TS_DiscGrad   *dg    = (TS_DiscGrad *) ts->data;
275   PetscReal      shift = 1/(0.5*ts->time_step);
276   Vec            X0, Xdot;
277   DM             dm, dmsave;
278   PetscErrorCode ierr;
279 
280   PetscFunctionBegin;
281   ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);
282   ierr = TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr);
283   ierr = VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x);CHKERRQ(ierr);
284 
285   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
286   dmsave = ts->dm;
287   ts->dm = dm;
288   ierr   = TSComputeIFunction(ts, dg->stage_time, x, Xdot, y, PETSC_FALSE);CHKERRQ(ierr);
289   ts->dm = dmsave;
290   ierr   = TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr);
291   PetscFunctionReturn(0);
292 }
293 
294 static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts)
295 {
296   TS_DiscGrad   *dg    = (TS_DiscGrad *) ts->data;
297   PetscReal      shift = 1/(0.5*ts->time_step);
298   Vec            Xdot;
299   DM             dm,dmsave;
300   PetscErrorCode ierr;
301 
302   PetscFunctionBegin;
303   ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);
304   /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */
305   ierr = TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot);CHKERRQ(ierr);
306 
307   dmsave = ts->dm;
308   ts->dm = dm;
309   ierr   = TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE);CHKERRQ(ierr);
310   ts->dm = dmsave;
311   ierr   = TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot);CHKERRQ(ierr);
312   PetscFunctionReturn(0);
313 }
314 
315 static PetscErrorCode TSDiscGradGetFormulation_DiscGrad(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *))
316 {
317   TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
318 
319   PetscFunctionBegin;
320   *Sfunc = dg->Sfunc;
321   *Ffunc = dg->Ffunc;
322   *Gfunc = dg->Gfunc;
323   PetscFunctionReturn(0);
324 }
325 
326 static PetscErrorCode TSDiscGradSetFormulation_DiscGrad(TS ts, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *))
327 {
328   TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;
329 
330   PetscFunctionBegin;
331   dg->Sfunc = Sfunc;
332   dg->Ffunc = Ffunc;
333   dg->Gfunc = Gfunc;
334   PetscFunctionReturn(0);
335 }
336 
337 /*MC
338   TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method
339 
340   Level: beginner
341 
342   Notes: This is the implicit midpoint rule, with an optional term that guarantees the discrete gradient property. This
343   timestepper applies to systems of the form
344 $ u_t = S(u) grad F(u)
345   where S(u) is a linear operator, and F is a functional of u.
346 
347 .seealso: TSCreate(), TSSetType(), TS, TSTHETA, TSDiscGradSetFormulation()
348 M*/
349 PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts)
350 {
351   TS_DiscGrad       *th;
352   PetscErrorCode ierr;
353 
354   PetscFunctionBegin;
355   ierr = PetscCitationsRegister(DGCitation, &DGCite);CHKERRQ(ierr);
356   ts->ops->reset           = TSReset_DiscGrad;
357   ts->ops->destroy         = TSDestroy_DiscGrad;
358   ts->ops->view            = TSView_DiscGrad;
359   ts->ops->setfromoptions  = TSSetFromOptions_DiscGrad;
360   ts->ops->setup           = TSSetUp_DiscGrad;
361   ts->ops->step            = TSStep_DiscGrad;
362   ts->ops->interpolate     = TSInterpolate_DiscGrad;
363   ts->ops->getstages       = TSGetStages_DiscGrad;
364   ts->ops->snesfunction    = SNESTSFormFunction_DiscGrad;
365   ts->ops->snesjacobian    = SNESTSFormJacobian_DiscGrad;
366   ts->default_adapt_type   = TSADAPTNONE;
367 
368   ts->usessnes = PETSC_TRUE;
369 
370   ierr = PetscNewLog(ts,&th);CHKERRQ(ierr);
371   ts->data = (void*)th;
372   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",TSDiscGradGetFormulation_DiscGrad);CHKERRQ(ierr);
373   ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",TSDiscGradSetFormulation_DiscGrad);CHKERRQ(ierr);
374   PetscFunctionReturn(0);
375 }
376 
377 /*@C
378   TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the formulation u_t = S grad F
379 
380   Not Collective
381 
382   Input Parameter:
383 . ts - timestepping context
384 
385   Output Parameters:
386 + Sfunc - constructor for the S matrix from the formulation
387 . Ffunc - functional F from the formulation
388 - Gfunc - constructor for the gradient of F from the formulation
389 
390   Calling sequence of Sfunc:
391 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *)
392 
393   Calling sequence of Ffunc:
394 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *)
395 
396   Calling sequence of Gfunc:
397 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *)
398 
399   Level: intermediate
400 
401 .seealso: TSDiscGradSetFormulation()
402 @*/
403 PetscErrorCode TSDiscGradGetFormulation(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *))
404 {
405   PetscErrorCode ierr;
406 
407   PetscFunctionBegin;
408   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
409   PetscValidPointer(Sfunc, 2);
410   PetscValidPointer(Ffunc, 3);
411   PetscValidPointer(Gfunc, 4);
412   ierr = PetscUseMethod(ts,"TSDiscGradGetFormulation_C",(TS,PetscErrorCode(**Sfunc)(TS,PetscReal,Vec,Mat,void*),PetscErrorCode(**Ffunc)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode(**Gfunc)(TS,PetscReal,Vec,Vec,void*)),(ts,Sfunc,Ffunc,Gfunc));CHKERRQ(ierr);
413   PetscFunctionReturn(0);
414 }
415 
416 /*@C
417   TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the formulation u_t = S(u) grad F(u)
418 
419   Not Collective
420 
421   Input Parameters:
422 + ts    - timestepping context
423 . Sfunc - constructor for the S matrix from the formulation
424 . Ffunc - functional F from the formulation
425 - Gfunc - constructor for the gradient of F from the formulation
426 
427   Calling sequence of Sfunc:
428 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *)
429 
430   Calling sequence of Ffunc:
431 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *)
432 
433   Calling sequence of Gfunc:
434 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *)
435 
436   Level: Intermediate
437 
438 .seealso: TSDiscGradGetFormulation()
439 @*/
440 PetscErrorCode TSDiscGradSetFormulation(TS ts, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec , PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx)
441 {
442   PetscErrorCode ierr;
443 
444   PetscFunctionBegin;
445   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
446   PetscValidFunction(Sfunc, 2);
447   PetscValidFunction(Ffunc, 3);
448   PetscValidFunction(Gfunc, 4);
449   ierr = PetscTryMethod(ts,"TSDiscGradSetFormulation_C",(TS,PetscErrorCode(*Sfunc)(TS,PetscReal,Vec,Mat,void*),PetscErrorCode(*Ffunc)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode(*Gfunc)(TS,PetscReal,Vec,Vec,void*)),(ts,Sfunc,Ffunc,Gfunc));CHKERRQ(ierr);
450   PetscFunctionReturn(0);
451 }
452