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: }

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

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

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

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

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

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

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

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

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

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

106: /*MC
107:       TSEULER - ODE solver using the explicit forward Euler method

109:   Level: beginner

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

117:   PetscFunctionBegin;
118:   PetscCall(PetscNew(&euler));
119:   ts->data = (void *)euler;

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