1: /*
  2:        Code for Timestepping with explicit Euler.
  3: */
  4: #include <petsc/private/tsimpl.h>                /*I   "petscts.h"   I*/

  6: typedef struct {
  7:   Vec update;     /* work vector where new solution is formed  */
  8: } TS_Euler;

 12: static PetscErrorCode TSStep_Euler(TS ts)
 13: {
 14:   TS_Euler       *euler = (TS_Euler*)ts->data;
 15:   Vec            solution = ts->vec_sol,update = euler->update;
 16:   PetscBool      stageok,accept = PETSC_TRUE;
 17:   PetscReal      next_time_step = ts->time_step;

 21:   TSPreStage(ts,ts->ptime);
 22:   TSComputeRHSFunction(ts,ts->ptime,solution,update);
 23:   VecAYPX(update,ts->time_step,solution);
 24:   TSPostStage(ts,ts->ptime,0,&solution);
 25:   TSAdaptCheckStage(ts->adapt,ts,ts->ptime,solution,&stageok);
 26:   if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
 27:   TSFunctionDomainError(ts,ts->ptime+ts->time_step,update,&stageok);
 28:   if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}

 30:   TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
 31:   if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
 32:   VecCopy(update,solution);

 34:   ts->ptime += ts->time_step;
 35:   ts->time_step = next_time_step;
 36:   return(0);
 37: }
 38: /*------------------------------------------------------------*/

 42: static PetscErrorCode TSSetUp_Euler(TS ts)
 43: {
 44:   TS_Euler       *euler = (TS_Euler*)ts->data;

 48:   TSCheckImplicitTerm(ts);
 49:   VecDuplicate(ts->vec_sol,&euler->update);
 50:   TSGetAdapt(ts,&ts->adapt);
 51:   TSAdaptCandidatesClear(ts->adapt);
 52:   return(0);
 53: }

 57: static PetscErrorCode TSReset_Euler(TS ts)
 58: {
 59:   TS_Euler       *euler = (TS_Euler*)ts->data;

 63:   VecDestroy(&euler->update);
 64:   return(0);
 65: }

 69: static PetscErrorCode TSDestroy_Euler(TS ts)
 70: {

 74:   TSReset_Euler(ts);
 75:   PetscFree(ts->data);
 76:   return(0);
 77: }
 78: /*------------------------------------------------------------*/

 82: static PetscErrorCode TSSetFromOptions_Euler(PetscOptionItems *PetscOptionsObject,TS ts)
 83: {
 85:   return(0);
 86: }

 90: static PetscErrorCode TSView_Euler(TS ts,PetscViewer viewer)
 91: {
 93:   return(0);
 94: }

 98: static PetscErrorCode TSInterpolate_Euler(TS ts,PetscReal t,Vec X)
 99: {
100:   TS_Euler       *euler = (TS_Euler*)ts->data;
101:   Vec            update = euler->update;
102:   PetscReal      alpha = (ts->ptime - t)/ts->time_step;

106:   VecWAXPY(X,-ts->time_step,update,ts->vec_sol);
107:   VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);
108:   return(0);
109: }

113: static PetscErrorCode TSComputeLinearStability_Euler(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
114: {
116:   *yr = 1.0 + xr;
117:   *yi = xi;
118:   return(0);
119: }
120: /* ------------------------------------------------------------ */

122: /*MC
123:       TSEULER - ODE solver using the explicit forward Euler method

125:   Level: beginner

127: .seealso:  TSCreate(), TS, TSSetType(), TSBEULER

129: M*/
132: PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts)
133: {
134:   TS_Euler       *euler;

138:   PetscNewLog(ts,&euler);
139:   ts->data = (void*)euler;

141:   ts->ops->setup           = TSSetUp_Euler;
142:   ts->ops->step            = TSStep_Euler;
143:   ts->ops->reset           = TSReset_Euler;
144:   ts->ops->destroy         = TSDestroy_Euler;
145:   ts->ops->setfromoptions  = TSSetFromOptions_Euler;
146:   ts->ops->view            = TSView_Euler;
147:   ts->ops->interpolate     = TSInterpolate_Euler;
148:   ts->ops->linearstability = TSComputeLinearStability_Euler;
149:   return(0);
150: }