Actual source code: theta.c

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

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

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

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


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

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

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

 71:   return(0);
 72: }

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

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

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

100:   return(0);
101: }

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

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

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

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

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

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

134:   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");
135:   if (order == th->order) {
136:     if (th->endpoint) {
137:       VecCopy(th->X,U);
138:     } else {
139:       PetscReal shift = 1./(th->Theta*ts->time_step);
140:       VecAXPBYPCZ(th->Xdot,-shift,shift,0,U,th->X);
141:       VecAXPY(U,ts->time_step,th->Xdot);
142:     }
143:   } else if (order == th->order-1 && order) {
144:     VecWAXPY(U,ts->time_step,th->Xdot,th->X0);
145:   }
146:   return(0);
147: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

313:   SNESGetDM(snes,&dm);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

463:    Level: beginner

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

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



477:    This method can be applied to DAE.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

551:   Not Collective

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

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

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

562:   Level: Advanced

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

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

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

582:   Not Collective

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

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

591:   Level: Intermediate

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

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

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

610:   Not Collective

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

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

618:   Level: Advanced

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

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

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

638:   Not Collective

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

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

647:   Level: Intermediate

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

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

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

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

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

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

680:   Level: beginner

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

685: $  -ts_type theta -ts_theta_theta 1.

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

689: M*/
692: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
693: {

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

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

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

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

717:   Level: beginner

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

722: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

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

726: M*/
729: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
730: {

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