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

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

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

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

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

 30:   ts->ptime += ts->time_step;
 31:   ts->time_step = next_time_step;
 32:   return 0;
 33: }
 34: /*------------------------------------------------------------*/

 36: static PetscErrorCode TSSetUp_Euler(TS ts)
 37: {
 38:   TS_Euler       *euler = (TS_Euler*)ts->data;

 40:   TSCheckImplicitTerm(ts);
 41:   VecDuplicate(ts->vec_sol,&euler->update);
 42:   TSGetAdapt(ts,&ts->adapt);
 43:   TSAdaptCandidatesClear(ts->adapt);
 44:   return 0;
 45: }

 47: static PetscErrorCode TSReset_Euler(TS ts)
 48: {
 49:   TS_Euler       *euler = (TS_Euler*)ts->data;

 51:   VecDestroy(&euler->update);
 52:   return 0;
 53: }

 55: static PetscErrorCode TSDestroy_Euler(TS ts)
 56: {
 57:   TSReset_Euler(ts);
 58:   PetscFree(ts->data);
 59:   return 0;
 60: }
 61: /*------------------------------------------------------------*/

 63: static PetscErrorCode TSSetFromOptions_Euler(PetscOptionItems *PetscOptionsObject,TS ts)
 64: {
 65:   return 0;
 66: }

 68: static PetscErrorCode TSView_Euler(TS ts,PetscViewer viewer)
 69: {
 70:   return 0;
 71: }

 73: static PetscErrorCode TSInterpolate_Euler(TS ts,PetscReal t,Vec X)
 74: {
 75:   TS_Euler       *euler = (TS_Euler*)ts->data;
 76:   Vec            update = euler->update;
 77:   PetscReal      alpha = (ts->ptime - t)/ts->time_step;

 79:   VecWAXPY(X,-ts->time_step,update,ts->vec_sol);
 80:   VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);
 81:   return 0;
 82: }

 84: static PetscErrorCode TSComputeLinearStability_Euler(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
 85: {
 86:   *yr = 1.0 + xr;
 87:   *yi = xi;
 88:   return 0;
 89: }
 90: /* ------------------------------------------------------------ */

 92: /*MC
 93:       TSEULER - ODE solver using the explicit forward Euler method

 95:   Level: beginner

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

 99: M*/
100: PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts)
101: {
102:   TS_Euler       *euler;

104:   PetscNewLog(ts,&euler);
105:   ts->data = (void*)euler;

107:   ts->ops->setup           = TSSetUp_Euler;
108:   ts->ops->step            = TSStep_Euler;
109:   ts->ops->reset           = TSReset_Euler;
110:   ts->ops->destroy         = TSDestroy_Euler;
111:   ts->ops->setfromoptions  = TSSetFromOptions_Euler;
112:   ts->ops->view            = TSView_Euler;
113:   ts->ops->interpolate     = TSInterpolate_Euler;
114:   ts->ops->linearstability = TSComputeLinearStability_Euler;
115:   ts->default_adapt_type   = TSADAPTNONE;
116:   ts->usessnes             = PETSC_FALSE;
117:   return 0;
118: }