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
TSStep_Euler(TS ts)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
TSSetUp_Euler(TS ts)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
TSReset_Euler(TS ts)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
TSDestroy_Euler(TS ts)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
TSSetFromOptions_Euler(TS ts,PetscOptionItems PetscOptionsObject)74 static PetscErrorCode TSSetFromOptions_Euler(TS ts, PetscOptionItems PetscOptionsObject)
75 {
76 PetscFunctionBegin;
77 PetscFunctionReturn(PETSC_SUCCESS);
78 }
79
TSView_Euler(TS ts,PetscViewer viewer)80 static PetscErrorCode TSView_Euler(TS ts, PetscViewer viewer)
81 {
82 PetscFunctionBegin;
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85
TSInterpolate_Euler(TS ts,PetscReal t,Vec X)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
TSComputeLinearStability_Euler(TS ts,PetscReal xr,PetscReal xi,PetscReal * yr,PetscReal * yi)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*/
TSCreate_Euler(TS ts)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