1 /* 2 Code for Timestepping with explicit Euler. 3 */ 4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 5 6 typedef struct { 7 Vec update; /* work vector where new solution is formed */ 8 } TS_Euler; 9 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; 16 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 } 32 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)); 39 40 ts->ptime += ts->time_step; 41 ts->time_step = next_time_step; 42 PetscFunctionReturn(PETSC_SUCCESS); 43 } 44 45 static PetscErrorCode TSSetUp_Euler(TS ts) 46 { 47 TS_Euler *euler = (TS_Euler *)ts->data; 48 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 } 56 57 static PetscErrorCode TSReset_Euler(TS ts) 58 { 59 TS_Euler *euler = (TS_Euler *)ts->data; 60 61 PetscFunctionBegin; 62 PetscCall(VecDestroy(&euler->update)); 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 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 } 73 74 static PetscErrorCode TSSetFromOptions_Euler(TS ts, PetscOptionItems PetscOptionsObject) 75 { 76 PetscFunctionBegin; 77 PetscFunctionReturn(PETSC_SUCCESS); 78 } 79 80 static PetscErrorCode TSView_Euler(TS ts, PetscViewer viewer) 81 { 82 PetscFunctionBegin; 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 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; 91 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 } 97 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 } 105 106 /*MC 107 TSEULER - ODE solver using the explicit forward Euler method 108 109 Level: beginner 110 111 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSBEULER`, `TSType` 112 M*/ 113 PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts) 114 { 115 TS_Euler *euler; 116 117 PetscFunctionBegin; 118 PetscCall(PetscNew(&euler)); 119 ts->data = (void *)euler; 120 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 } 133