Actual source code: euler.c

  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) {
 23:     ts->reason = TS_DIVERGED_STEP_REJECTED;
 24:     return 0;
 25:   }
 26:   TSFunctionDomainError(ts, ts->ptime + ts->time_step, update, &stageok);
 27:   if (!stageok) {
 28:     ts->reason = TS_DIVERGED_STEP_REJECTED;
 29:     return 0;
 30:   }

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

 39:   ts->ptime += ts->time_step;
 40:   ts->time_step = next_time_step;
 41:   return 0;
 42: }
 43: /*------------------------------------------------------------*/

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

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

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

 60:   VecDestroy(&euler->update);
 61:   return 0;
 62: }

 64: static PetscErrorCode TSDestroy_Euler(TS ts)
 65: {
 66:   TSReset_Euler(ts);
 67:   PetscFree(ts->data);
 68:   return 0;
 69: }
 70: /*------------------------------------------------------------*/

 72: static PetscErrorCode TSSetFromOptions_Euler(TS ts, PetscOptionItems *PetscOptionsObject)
 73: {
 74:   return 0;
 75: }

 77: static PetscErrorCode TSView_Euler(TS ts, PetscViewer viewer)
 78: {
 79:   return 0;
 80: }

 82: static PetscErrorCode TSInterpolate_Euler(TS ts, PetscReal t, Vec X)
 83: {
 84:   TS_Euler *euler  = (TS_Euler *)ts->data;
 85:   Vec       update = euler->update;
 86:   PetscReal alpha  = (ts->ptime - t) / ts->time_step;

 88:   VecWAXPY(X, -ts->time_step, update, ts->vec_sol);
 89:   VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol);
 90:   return 0;
 91: }

 93: static PetscErrorCode TSComputeLinearStability_Euler(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
 94: {
 95:   *yr = 1.0 + xr;
 96:   *yi = xi;
 97:   return 0;
 98: }
 99: /* ------------------------------------------------------------ */

101: /*MC
102:       TSEULER - ODE solver using the explicit forward Euler method

104:   Level: beginner

106: .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSBEULER`

108: M*/
109: PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts)
110: {
111:   TS_Euler *euler;

113:   PetscNew(&euler);
114:   ts->data = (void *)euler;

116:   ts->ops->setup           = TSSetUp_Euler;
117:   ts->ops->step            = TSStep_Euler;
118:   ts->ops->reset           = TSReset_Euler;
119:   ts->ops->destroy         = TSDestroy_Euler;
120:   ts->ops->setfromoptions  = TSSetFromOptions_Euler;
121:   ts->ops->view            = TSView_Euler;
122:   ts->ops->interpolate     = TSInterpolate_Euler;
123:   ts->ops->linearstability = TSComputeLinearStability_Euler;
124:   ts->default_adapt_type   = TSADAPTNONE;
125:   ts->usessnes             = PETSC_FALSE;
126:   return 0;
127: }