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:   PetscFunctionBegin;
 18:   PetscCall(TSPreStage(ts, ts->ptime));
 19:   PetscCall(TSComputeRHSFunction(ts, ts->ptime, solution, update));
 20:   PetscCall(VecAYPX(update, ts->time_step, solution));
 21:   PetscCall(TSPostStage(ts, ts->ptime, 0, &solution));
 22:   PetscCall(TSAdaptCheckStage(ts->adapt, ts, ts->ptime, solution, &stageok));
 23:   if (!stageok) {
 24:     ts->reason = TS_DIVERGED_STEP_REJECTED;
 25:     PetscFunctionReturn(PETSC_SUCCESS);
 26:   }
 27:   PetscCall(TSFunctionDomainError(ts, ts->ptime + ts->time_step, update, &stageok));
 28:   if (!stageok) {
 29:     ts->reason = TS_DIVERGED_STEP_REJECTED;
 30:     PetscFunctionReturn(PETSC_SUCCESS);
 31:   }

 33:   PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
 34:   if (!accept) {
 35:     ts->reason = TS_DIVERGED_STEP_REJECTED;
 36:     PetscFunctionReturn(PETSC_SUCCESS);
 37:   }
 38:   PetscCall(VecCopy(update, solution));

 40:   ts->ptime += ts->time_step;
 41:   ts->time_step = next_time_step;
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }
 44: /*------------------------------------------------------------*/

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

 50:   PetscFunctionBegin;
 51:   PetscCall(TSCheckImplicitTerm(ts));
 52:   PetscCall(VecDuplicate(ts->vec_sol, &euler->update));
 53:   PetscCall(TSGetAdapt(ts, &ts->adapt));
 54:   PetscCall(TSAdaptCandidatesClear(ts->adapt));
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

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

 62:   PetscFunctionBegin;
 63:   PetscCall(VecDestroy(&euler->update));
 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: static PetscErrorCode TSDestroy_Euler(TS ts)
 68: {
 69:   PetscFunctionBegin;
 70:   PetscCall(TSReset_Euler(ts));
 71:   PetscCall(PetscFree(ts->data));
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }
 74: /*------------------------------------------------------------*/

 76: static PetscErrorCode TSSetFromOptions_Euler(TS ts, PetscOptionItems *PetscOptionsObject)
 77: {
 78:   PetscFunctionBegin;
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: static PetscErrorCode TSView_Euler(TS ts, PetscViewer viewer)
 83: {
 84:   PetscFunctionBegin;
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

 88: static PetscErrorCode TSInterpolate_Euler(TS ts, PetscReal t, Vec X)
 89: {
 90:   TS_Euler *euler  = (TS_Euler *)ts->data;
 91:   Vec       update = euler->update;
 92:   PetscReal alpha  = (ts->ptime - t) / ts->time_step;

 94:   PetscFunctionBegin;
 95:   PetscCall(VecWAXPY(X, -ts->time_step, update, ts->vec_sol));
 96:   PetscCall(VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol));
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: static PetscErrorCode TSComputeLinearStability_Euler(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
101: {
102:   PetscFunctionBegin;
103:   *yr = 1.0 + xr;
104:   *yi = xi;
105:   PetscFunctionReturn(PETSC_SUCCESS);
106: }
107: /* ------------------------------------------------------------ */

109: /*MC
110:       TSEULER - ODE solver using the explicit forward Euler method

112:   Level: beginner

114: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSBEULER`, `TSType`
115: M*/
116: PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts)
117: {
118:   TS_Euler *euler;

120:   PetscFunctionBegin;
121:   PetscCall(PetscNew(&euler));
122:   ts->data = (void *)euler;

124:   ts->ops->setup           = TSSetUp_Euler;
125:   ts->ops->step            = TSStep_Euler;
126:   ts->ops->reset           = TSReset_Euler;
127:   ts->ops->destroy         = TSDestroy_Euler;
128:   ts->ops->setfromoptions  = TSSetFromOptions_Euler;
129:   ts->ops->view            = TSView_Euler;
130:   ts->ops->interpolate     = TSInterpolate_Euler;
131:   ts->ops->linearstability = TSComputeLinearStability_Euler;
132:   ts->default_adapt_type   = TSADAPTNONE;
133:   ts->usessnes             = PETSC_FALSE;
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }