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 TS_Euler *euler = (TS_Euler *)ts->data; 12 Vec solution = ts->vec_sol, update = euler->update; 13 PetscBool stageok, accept = PETSC_TRUE; 14 PetscReal next_time_step = ts->time_step; 15 16 PetscFunctionBegin; 17 PetscCall(TSPreStage(ts, ts->ptime)); 18 PetscCall(TSComputeRHSFunction(ts, ts->ptime, solution, update)); 19 PetscCall(VecAYPX(update, ts->time_step, solution)); 20 PetscCall(TSPostStage(ts, ts->ptime, 0, &solution)); 21 PetscCall(TSAdaptCheckStage(ts->adapt, ts, ts->ptime, solution, &stageok)); 22 if (!stageok) { 23 ts->reason = TS_DIVERGED_STEP_REJECTED; 24 PetscFunctionReturn(0); 25 } 26 PetscCall(TSFunctionDomainError(ts, ts->ptime + ts->time_step, update, &stageok)); 27 if (!stageok) { 28 ts->reason = TS_DIVERGED_STEP_REJECTED; 29 PetscFunctionReturn(0); 30 } 31 32 PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 33 if (!accept) { 34 ts->reason = TS_DIVERGED_STEP_REJECTED; 35 PetscFunctionReturn(0); 36 } 37 PetscCall(VecCopy(update, solution)); 38 39 ts->ptime += ts->time_step; 40 ts->time_step = next_time_step; 41 PetscFunctionReturn(0); 42 } 43 /*------------------------------------------------------------*/ 44 45 static PetscErrorCode TSSetUp_Euler(TS ts) { 46 TS_Euler *euler = (TS_Euler *)ts->data; 47 48 PetscFunctionBegin; 49 PetscCall(TSCheckImplicitTerm(ts)); 50 PetscCall(VecDuplicate(ts->vec_sol, &euler->update)); 51 PetscCall(TSGetAdapt(ts, &ts->adapt)); 52 PetscCall(TSAdaptCandidatesClear(ts->adapt)); 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode TSReset_Euler(TS ts) { 57 TS_Euler *euler = (TS_Euler *)ts->data; 58 59 PetscFunctionBegin; 60 PetscCall(VecDestroy(&euler->update)); 61 PetscFunctionReturn(0); 62 } 63 64 static PetscErrorCode TSDestroy_Euler(TS ts) { 65 PetscFunctionBegin; 66 PetscCall(TSReset_Euler(ts)); 67 PetscCall(PetscFree(ts->data)); 68 PetscFunctionReturn(0); 69 } 70 /*------------------------------------------------------------*/ 71 72 static PetscErrorCode TSSetFromOptions_Euler(TS ts, PetscOptionItems *PetscOptionsObject) { 73 PetscFunctionBegin; 74 PetscFunctionReturn(0); 75 } 76 77 static PetscErrorCode TSView_Euler(TS ts, PetscViewer viewer) { 78 PetscFunctionBegin; 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode TSInterpolate_Euler(TS ts, PetscReal t, Vec X) { 83 TS_Euler *euler = (TS_Euler *)ts->data; 84 Vec update = euler->update; 85 PetscReal alpha = (ts->ptime - t) / ts->time_step; 86 87 PetscFunctionBegin; 88 PetscCall(VecWAXPY(X, -ts->time_step, update, ts->vec_sol)); 89 PetscCall(VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol)); 90 PetscFunctionReturn(0); 91 } 92 93 static PetscErrorCode TSComputeLinearStability_Euler(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi) { 94 PetscFunctionBegin; 95 *yr = 1.0 + xr; 96 *yi = xi; 97 PetscFunctionReturn(0); 98 } 99 /* ------------------------------------------------------------ */ 100 101 /*MC 102 TSEULER - ODE solver using the explicit forward Euler method 103 104 Level: beginner 105 106 .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSBEULER` 107 108 M*/ 109 PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts) { 110 TS_Euler *euler; 111 112 PetscFunctionBegin; 113 PetscCall(PetscNewLog(ts, &euler)); 114 ts->data = (void *)euler; 115 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 PetscFunctionReturn(0); 127 } 128