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