Actual source code: theta.c

petsc-dev 2014-07-12
Report Typos and Errors
  1: /*
  2:   Code for timestepping with implicit Theta method
  3: */
  4: #define PETSC_DESIRE_COMPLEX
  5: #include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/
  6: #include <petscsnes.h>
  7: #include <petscdm.h>

  9: typedef struct {
 10:   Vec          X,Xdot;                   /* Storage for one stage */
 11:   Vec          X0;                       /* work vector to store X0 */
 12:   Vec          affine;                   /* Affine vector needed for residual at beginning of step */
 13:   PetscBool    extrapolate;
 14:   PetscBool    endpoint;
 15:   PetscReal    Theta;
 16:   PetscReal    stage_time;
 17:   TSStepStatus status;
 18:   char         *name;
 19:   PetscInt     order;
 20:   PetscReal    ccfl;               /* Placeholder for CFL coefficient relative to forward Euler */
 21:   PetscBool    adapt;  /* use time-step adaptivity ? */
 22: } TS_Theta;

 26: static PetscErrorCode TSThetaGetX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
 27: {
 28:   TS_Theta       *th = (TS_Theta*)ts->data;

 32:   if (X0) {
 33:     if (dm && dm != ts->dm) {
 34:       DMGetNamedGlobalVector(dm,"TSTheta_X0",X0);
 35:     } else *X0 = ts->vec_sol;
 36:   }
 37:   if (Xdot) {
 38:     if (dm && dm != ts->dm) {
 39:       DMGetNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
 40:     } else *Xdot = th->Xdot;
 41:   }
 42:   return(0);
 43: }


 48: static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
 49: {

 53:   if (X0) {
 54:     if (dm && dm != ts->dm) {
 55:       DMRestoreNamedGlobalVector(dm,"TSTheta_X0",X0);
 56:     }
 57:   }
 58:   if (Xdot) {
 59:     if (dm && dm != ts->dm) {
 60:       DMRestoreNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
 61:     }
 62:   }
 63:   return(0);
 64: }

 68: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine,DM coarse,void *ctx)
 69: {

 72:   return(0);
 73: }

 77: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
 78: {
 79:   TS             ts = (TS)ctx;
 81:   Vec            X0,Xdot,X0_c,Xdot_c;

 84:   TSThetaGetX0AndXdot(ts,fine,&X0,&Xdot);
 85:   TSThetaGetX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
 86:   MatRestrict(restrct,X0,X0_c);
 87:   MatRestrict(restrct,Xdot,Xdot_c);
 88:   VecPointwiseMult(X0_c,rscale,X0_c);
 89:   VecPointwiseMult(Xdot_c,rscale,Xdot_c);
 90:   TSThetaRestoreX0AndXdot(ts,fine,&X0,&Xdot);
 91:   TSThetaRestoreX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
 92:   return(0);
 93: }

 97: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm,DM subdm,void *ctx)
 98: {

101:   return(0);
102: }

106: static PetscErrorCode DMSubDomainRestrictHook_TSTheta(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
107: {
108:   TS             ts = (TS)ctx;
110:   Vec            X0,Xdot,X0_sub,Xdot_sub;

113:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
114:   TSThetaGetX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);

116:   VecScatterBegin(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
117:   VecScatterEnd(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);

119:   VecScatterBegin(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
120:   VecScatterEnd(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);

122:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
123:   TSThetaRestoreX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
124:   return(0);
125: }

129: static PetscErrorCode TSEvaluateStep_Theta(TS ts,PetscInt order,Vec U,PetscBool *done)
130: {
132:   TS_Theta       *th = (TS_Theta*)ts->data;

135:   if (order == 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"No time-step adaptivity implemented for 1st order theta method; Run with -ts_adapt_type none");
136:   if (order == th->order) {
137:     if (th->endpoint) {
138:       VecCopy(th->X,U);
139:     } else {
140:       PetscReal shift = 1./(th->Theta*ts->time_step);
141:       VecAXPBYPCZ(th->Xdot,-shift,shift,0,U,th->X);
142:       VecAXPY(U,ts->time_step,th->Xdot);
143:     }
144:   } else if (order == th->order-1 && order) {
145:     VecWAXPY(U,ts->time_step,th->Xdot,th->X0);
146:   }
147:   return(0);
148: }

152: static PetscErrorCode TSRollBack_Theta(TS ts)
153: {
154:   TS_Theta       *th = (TS_Theta*)ts->data;

158:   VecCopy(th->X0,ts->vec_sol);
159:   th->status    = TS_STEP_INCOMPLETE;
160:   return(0);
161: }

165: static PetscErrorCode TSStep_Theta(TS ts)
166: {
167:   TS_Theta            *th = (TS_Theta*)ts->data;
168:   PetscInt            its,lits,reject,next_scheme;
169:   PetscReal           next_time_step;
170:   SNESConvergedReason snesreason;
171:   PetscErrorCode      ierr;
172:   TSAdapt             adapt;
173:   PetscBool           accept;

176:   next_time_step = ts->time_step;
177:   th->status = TS_STEP_INCOMPLETE;
178:   VecCopy(ts->vec_sol,th->X0);
179:   for (reject=0; reject<ts->max_reject && !ts->reason && th->status != TS_STEP_COMPLETE; reject++,ts->reject++) {
180:     PetscReal shift = 1./(th->Theta*ts->time_step);
181:     th->stage_time = ts->ptime + (th->endpoint ? 1. : th->Theta)*ts->time_step;
182:     TSPreStep(ts);
183:     TSPreStage(ts,th->stage_time);

185:     if (th->endpoint) {           /* This formulation assumes linear time-independent mass matrix */
186:       VecZeroEntries(th->Xdot);
187:       if (!th->affine) {VecDuplicate(ts->vec_sol,&th->affine);}
188:       TSComputeIFunction(ts,ts->ptime,ts->vec_sol,th->Xdot,th->affine,PETSC_FALSE);
189:       VecScale(th->affine,(th->Theta-1.)/th->Theta);
190:     }
191:     if (th->extrapolate) {
192:       VecWAXPY(th->X,1./shift,th->Xdot,ts->vec_sol);
193:     } else {
194:       VecCopy(ts->vec_sol,th->X);
195:     }
196:     SNESSolve(ts->snes,th->affine,th->X);
197:     SNESGetIterationNumber(ts->snes,&its);
198:     SNESGetLinearSolveIterations(ts->snes,&lits);
199:     SNESGetConvergedReason(ts->snes,&snesreason);
200:     TSPostStage(ts,th->stage_time,0,&(th->X));
201:     ts->snes_its += its; ts->ksp_its += lits;
202:     TSGetAdapt(ts,&adapt);
203:     TSAdaptCheckStage(adapt,ts,&accept);
204:     if (!accept) continue;
205:     TSEvaluateStep(ts,th->order,ts->vec_sol,NULL);
206:     /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
207:     TSGetAdapt(ts,&adapt);
208:     TSAdaptCandidatesClear(adapt);
209:     TSAdaptCandidateAdd(adapt,NULL,th->order,1,th->ccfl,1.0,PETSC_TRUE);
210:     TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);

212:     if (accept) {
213:       /* ignore next_scheme for now */
214:       ts->ptime    += ts->time_step;
215:       ts->time_step = next_time_step;
216:       ts->steps++;
217:       th->status = TS_STEP_COMPLETE;
218:     } else {                    /* Roll back the current step */
219:       ts->ptime += next_time_step; /* This will be undone in rollback */
220:       th->status = TS_STEP_INCOMPLETE;
221:       TSRollBack(ts);
222:     }
223:   }
224:   return(0);
225: }

229: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
230: {
231:   TS_Theta       *th   = (TS_Theta*)ts->data;
232:   PetscReal      alpha = t - ts->ptime;

236:   VecCopy(ts->vec_sol,th->X);
237:   if (th->endpoint) alpha *= th->Theta;
238:   VecWAXPY(X,alpha,th->Xdot,th->X);
239:   return(0);
240: }

242: /*------------------------------------------------------------*/
245: static PetscErrorCode TSReset_Theta(TS ts)
246: {
247:   TS_Theta       *th = (TS_Theta*)ts->data;

251:   VecDestroy(&th->X);
252:   VecDestroy(&th->Xdot);
253:   VecDestroy(&th->X0);
254:   VecDestroy(&th->affine);
255:   return(0);
256: }

260: static PetscErrorCode TSDestroy_Theta(TS ts)
261: {

265:   TSReset_Theta(ts);
266:   PetscFree(ts->data);
267:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
268:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
269:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
270:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
271:   return(0);
272: }

274: /*
275:   This defines the nonlinear equation that is to be solved with SNES
276:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
277: */
280: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
281: {
282:   TS_Theta       *th = (TS_Theta*)ts->data;
284:   Vec            X0,Xdot;
285:   DM             dm,dmsave;
286:   PetscReal      shift = 1./(th->Theta*ts->time_step);

289:   SNESGetDM(snes,&dm);
290:   /* When using the endpoint variant, this is actually 1/Theta * Xdot */
291:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
292:   VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);

294:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
295:   dmsave = ts->dm;
296:   ts->dm = dm;
297:   TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
298:   ts->dm = dmsave;
299:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
300:   return(0);
301: }

305: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
306: {
307:   TS_Theta       *th = (TS_Theta*)ts->data;
309:   Vec            Xdot;
310:   DM             dm,dmsave;
311:   PetscReal      shift = 1./(th->Theta*ts->time_step);

314:   SNESGetDM(snes,&dm);

316:   /* th->Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
317:   TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);

319:   dmsave = ts->dm;
320:   ts->dm = dm;
321:   TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
322:   ts->dm = dmsave;
323:   TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
324:   return(0);
325: }

329: static PetscErrorCode TSSetUp_Theta(TS ts)
330: {
331:   TS_Theta       *th = (TS_Theta*)ts->data;
333:   SNES           snes;
334:   DM             dm;

337:   VecDuplicate(ts->vec_sol,&th->X);
338:   VecDuplicate(ts->vec_sol,&th->Xdot);
339:   VecDuplicate(ts->vec_sol,&th->X0);
340:   TSGetSNES(ts,&snes);
341:   TSGetDM(ts,&dm);
342:   if (dm) {
343:     DMCoarsenHookAdd(dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
344:     DMSubDomainHookAdd(dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
345:   }
346:   if (th->Theta == 0.5 && th->endpoint) th->order = 2;
347:   else th->order = 1;

349:   if (!th->adapt) {
350:     TSAdapt adapt;
351:     TSAdaptDestroy(&ts->adapt);
352:     TSGetAdapt(ts,&adapt);
353:     TSAdaptSetType(adapt,TSADAPTNONE);
354:   }
355:   return(0);
356: }
357: /*------------------------------------------------------------*/

361: static PetscErrorCode TSSetFromOptions_Theta(TS ts)
362: {
363:   TS_Theta       *th = (TS_Theta*)ts->data;

367:   PetscOptionsHead("Theta ODE solver options");
368:   {
369:     PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
370:     PetscOptionsBool("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
371:     PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
372:     PetscOptionsBool("-ts_theta_adapt","Use time-step adaptivity with the Theta method","",th->adapt,&th->adapt,NULL);
373:     SNESSetFromOptions(ts->snes);
374:   }
375:   PetscOptionsTail();
376:   return(0);
377: }

381: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
382: {
383:   TS_Theta       *th = (TS_Theta*)ts->data;
384:   PetscBool      iascii;

388:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
389:   if (iascii) {
390:     PetscViewerASCIIPrintf(viewer,"  Theta=%g\n",(double)th->Theta);
391:     PetscViewerASCIIPrintf(viewer,"  Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
392:   }
393:   SNESView(ts->snes,viewer);
394:   return(0);
395: }

399: PetscErrorCode  TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
400: {
401:   TS_Theta *th = (TS_Theta*)ts->data;

404:   *theta = th->Theta;
405:   return(0);
406: }

410: PetscErrorCode  TSThetaSetTheta_Theta(TS ts,PetscReal theta)
411: {
412:   TS_Theta *th = (TS_Theta*)ts->data;

415:   if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
416:   th->Theta = theta;
417:   return(0);
418: }

422: PetscErrorCode  TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
423: {
424:   TS_Theta *th = (TS_Theta*)ts->data;

427:   *endpoint = th->endpoint;
428:   return(0);
429: }

433: PetscErrorCode  TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
434: {
435:   TS_Theta *th = (TS_Theta*)ts->data;

438:   th->endpoint = flg;
439:   return(0);
440: }

442: #if defined(PETSC_HAVE_COMPLEX)
445: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
446: {
447:   PetscComplex z   = xr + xi*PETSC_i,f;
448:   TS_Theta     *th = (TS_Theta*)ts->data;
449:   const PetscReal one = 1.0;

452:   f   = (one + (one - th->Theta)*z)/(one - th->Theta*z);
453:   *yr = PetscRealPartComplex(f);
454:   *yi = PetscImaginaryPartComplex(f);
455:   return(0);
456: }
457: #endif


460: /* ------------------------------------------------------------ */
461: /*MC
462:       TSTHETA - DAE solver using the implicit Theta method

464:    Level: beginner

466:    Options Database:
467:       -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
468:       -ts_theta_extrapolate <flg> Extrapolate stage solution from previous solution (sometimes unstable)
469:       -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of the Theta method

471:    Notes:
472: $  -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
473: $  -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
474: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)



478:    This method can be applied to DAE.

480:    This method is cast as a 1-stage implicit Runge-Kutta method.

482: .vb
483:   Theta | Theta
484:   -------------
485:         |  1
486: .ve

488:    For the default Theta=0.5, this is also known as the implicit midpoint rule.

490:    When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:

492: .vb
493:   0 | 0         0
494:   1 | 1-Theta   Theta
495:   -------------------
496:     | 1-Theta   Theta
497: .ve

499:    For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).

501:    To apply a diagonally implicit RK method to DAE, the stage formula

503: $  Y_i = X + h sum_j a_ij Y'_j

505:    is interpreted as a formula for Y'_i in terms of Y_i and known values (Y'_j, j<i)

507: .seealso:  TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()

509: M*/
512: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
513: {
514:   TS_Theta       *th;

518:   ts->ops->reset          = TSReset_Theta;
519:   ts->ops->destroy        = TSDestroy_Theta;
520:   ts->ops->view           = TSView_Theta;
521:   ts->ops->setup          = TSSetUp_Theta;
522:   ts->ops->step           = TSStep_Theta;
523:   ts->ops->interpolate    = TSInterpolate_Theta;
524:   ts->ops->evaluatestep   = TSEvaluateStep_Theta;
525:   ts->ops->rollback       = TSRollBack_Theta;
526:   ts->ops->setfromoptions = TSSetFromOptions_Theta;
527:   ts->ops->snesfunction   = SNESTSFormFunction_Theta;
528:   ts->ops->snesjacobian   = SNESTSFormJacobian_Theta;
529: #if defined(PETSC_HAVE_COMPLEX)
530:   ts->ops->linearstability = TSComputeLinearStability_Theta;
531: #endif

533:   PetscNewLog(ts,&th);
534:   ts->data = (void*)th;

536:   th->extrapolate = PETSC_FALSE;
537:   th->Theta       = 0.5;
538:   th->ccfl        = 1.0;
539:   th->adapt       = PETSC_FALSE;
540:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
541:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
542:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
543:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
544:   return(0);
545: }

549: /*@
550:   TSThetaGetTheta - Get the abscissa of the stage in (0,1].

552:   Not Collective

554:   Input Parameter:
555: .  ts - timestepping context

557:   Output Parameter:
558: .  theta - stage abscissa

560:   Note:
561:   Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.

563:   Level: Advanced

565: .seealso: TSThetaSetTheta()
566: @*/
567: PetscErrorCode  TSThetaGetTheta(TS ts,PetscReal *theta)
568: {

574:   PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
575:   return(0);
576: }

580: /*@
581:   TSThetaSetTheta - Set the abscissa of the stage in (0,1].

583:   Not Collective

585:   Input Parameter:
586: +  ts - timestepping context
587: -  theta - stage abscissa

589:   Options Database:
590: .  -ts_theta_theta <theta>

592:   Level: Intermediate

594: .seealso: TSThetaGetTheta()
595: @*/
596: PetscErrorCode  TSThetaSetTheta(TS ts,PetscReal theta)
597: {

602:   PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
603:   return(0);
604: }

608: /*@
609:   TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).

611:   Not Collective

613:   Input Parameter:
614: .  ts - timestepping context

616:   Output Parameter:
617: .  endpoint - PETSC_TRUE when using the endpoint variant

619:   Level: Advanced

621: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
622: @*/
623: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
624: {

630:   PetscTryMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
631:   return(0);
632: }

636: /*@
637:   TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).

639:   Not Collective

641:   Input Parameter:
642: +  ts - timestepping context
643: -  flg - PETSC_TRUE to use the endpoint variant

645:   Options Database:
646: .  -ts_theta_endpoint <flg>

648:   Level: Intermediate

650: .seealso: TSTHETA, TSCN
651: @*/
652: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
653: {

658:   PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
659:   return(0);
660: }

662: /*
663:  * TSBEULER and TSCN are straightforward specializations of TSTHETA.
664:  * The creation functions for these specializations are below.
665:  */

669: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
670: {

674:   SNESView(ts->snes,viewer);
675:   return(0);
676: }

678: /*MC
679:       TSBEULER - ODE solver using the implicit backward Euler method

681:   Level: beginner

683:   Notes:
684:   TSBEULER is equivalent to TSTHETA with Theta=1.0

686: $  -ts_type theta -ts_theta_theta 1.

688: .seealso:  TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA

690: M*/
693: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
694: {

698:   TSCreate_Theta(ts);
699:   TSThetaSetTheta(ts,1.0);
700:   ts->ops->view = TSView_BEuler;
701:   return(0);
702: }

706: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
707: {

711:   SNESView(ts->snes,viewer);
712:   return(0);
713: }

715: /*MC
716:       TSCN - ODE solver using the implicit Crank-Nicolson method.

718:   Level: beginner

720:   Notes:
721:   TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.

723: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

725: .seealso:  TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA

727: M*/
730: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
731: {

735:   TSCreate_Theta(ts);
736:   TSThetaSetTheta(ts,0.5);
737:   TSThetaSetEndpoint(ts,PETSC_TRUE);
738:   ts->ops->view = TSView_CN;
739:   return(0);
740: }