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 46 static PetscErrorCode TSSetUp_Euler(TS ts) 47 { 48 TS_Euler *euler = (TS_Euler *)ts->data; 49 50 PetscFunctionBegin; 51 PetscCall(TSCheckImplicitTerm(ts)); 52 PetscCall(VecDuplicate(ts->vec_sol, &euler->update)); 53 PetscCall(TSGetAdapt(ts, &ts->adapt)); 54 PetscCall(TSAdaptCandidatesClear(ts->adapt)); 55 PetscFunctionReturn(PETSC_SUCCESS); 56 } 57 58 static PetscErrorCode TSReset_Euler(TS ts) 59 { 60 TS_Euler *euler = (TS_Euler *)ts->data; 61 62 PetscFunctionBegin; 63 PetscCall(VecDestroy(&euler->update)); 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 static PetscErrorCode TSDestroy_Euler(TS ts) 68 { 69 PetscFunctionBegin; 70 PetscCall(TSReset_Euler(ts)); 71 PetscCall(PetscFree(ts->data)); 72 PetscFunctionReturn(PETSC_SUCCESS); 73 } 74 /*------------------------------------------------------------*/ 75 76 static PetscErrorCode TSSetFromOptions_Euler(TS ts, PetscOptionItems *PetscOptionsObject) 77 { 78 PetscFunctionBegin; 79 PetscFunctionReturn(PETSC_SUCCESS); 80 } 81 82 static PetscErrorCode TSView_Euler(TS ts, PetscViewer viewer) 83 { 84 PetscFunctionBegin; 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 static PetscErrorCode TSInterpolate_Euler(TS ts, PetscReal t, Vec X) 89 { 90 TS_Euler *euler = (TS_Euler *)ts->data; 91 Vec update = euler->update; 92 PetscReal alpha = (ts->ptime - t) / ts->time_step; 93 94 PetscFunctionBegin; 95 PetscCall(VecWAXPY(X, -ts->time_step, update, ts->vec_sol)); 96 PetscCall(VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol)); 97 PetscFunctionReturn(PETSC_SUCCESS); 98 } 99 100 static PetscErrorCode TSComputeLinearStability_Euler(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi) 101 { 102 PetscFunctionBegin; 103 *yr = 1.0 + xr; 104 *yi = xi; 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107 /* ------------------------------------------------------------ */ 108 109 /*MC 110 TSEULER - ODE solver using the explicit forward Euler method 111 112 Level: beginner 113 114 .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSBEULER`, `TSType` 115 M*/ 116 PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts) 117 { 118 TS_Euler *euler; 119 120 PetscFunctionBegin; 121 PetscCall(PetscNew(&euler)); 122 ts->data = (void *)euler; 123 124 ts->ops->setup = TSSetUp_Euler; 125 ts->ops->step = TSStep_Euler; 126 ts->ops->reset = TSReset_Euler; 127 ts->ops->destroy = TSDestroy_Euler; 128 ts->ops->setfromoptions = TSSetFromOptions_Euler; 129 ts->ops->view = TSView_Euler; 130 ts->ops->interpolate = TSInterpolate_Euler; 131 ts->ops->linearstability = TSComputeLinearStability_Euler; 132 ts->default_adapt_type = TSADAPTNONE; 133 ts->usessnes = PETSC_FALSE; 134 PetscFunctionReturn(PETSC_SUCCESS); 135 } 136