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: TSPreStage(ts, ts->ptime);
18: TSComputeRHSFunction(ts, ts->ptime, solution, update);
19: VecAYPX(update, ts->time_step, solution);
20: TSPostStage(ts, ts->ptime, 0, &solution);
21: TSAdaptCheckStage(ts->adapt, ts, ts->ptime, solution, &stageok);
22: if (!stageok) {
23: ts->reason = TS_DIVERGED_STEP_REJECTED;
24: return 0;
25: }
26: TSFunctionDomainError(ts, ts->ptime + ts->time_step, update, &stageok);
27: if (!stageok) {
28: ts->reason = TS_DIVERGED_STEP_REJECTED;
29: return 0;
30: }
32: TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
33: if (!accept) {
34: ts->reason = TS_DIVERGED_STEP_REJECTED;
35: return 0;
36: }
37: VecCopy(update, solution);
39: ts->ptime += ts->time_step;
40: ts->time_step = next_time_step;
41: return 0;
42: }
43: /*------------------------------------------------------------*/
45: static PetscErrorCode TSSetUp_Euler(TS ts)
46: {
47: TS_Euler *euler = (TS_Euler *)ts->data;
49: TSCheckImplicitTerm(ts);
50: VecDuplicate(ts->vec_sol, &euler->update);
51: TSGetAdapt(ts, &ts->adapt);
52: TSAdaptCandidatesClear(ts->adapt);
53: return 0;
54: }
56: static PetscErrorCode TSReset_Euler(TS ts)
57: {
58: TS_Euler *euler = (TS_Euler *)ts->data;
60: VecDestroy(&euler->update);
61: return 0;
62: }
64: static PetscErrorCode TSDestroy_Euler(TS ts)
65: {
66: TSReset_Euler(ts);
67: PetscFree(ts->data);
68: return 0;
69: }
70: /*------------------------------------------------------------*/
72: static PetscErrorCode TSSetFromOptions_Euler(TS ts, PetscOptionItems *PetscOptionsObject)
73: {
74: return 0;
75: }
77: static PetscErrorCode TSView_Euler(TS ts, PetscViewer viewer)
78: {
79: return 0;
80: }
82: static PetscErrorCode TSInterpolate_Euler(TS ts, PetscReal t, Vec X)
83: {
84: TS_Euler *euler = (TS_Euler *)ts->data;
85: Vec update = euler->update;
86: PetscReal alpha = (ts->ptime - t) / ts->time_step;
88: VecWAXPY(X, -ts->time_step, update, ts->vec_sol);
89: VecAXPBY(X, 1.0 - alpha, alpha, ts->vec_sol);
90: return 0;
91: }
93: static PetscErrorCode TSComputeLinearStability_Euler(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
94: {
95: *yr = 1.0 + xr;
96: *yi = xi;
97: return 0;
98: }
99: /* ------------------------------------------------------------ */
101: /*MC
102: TSEULER - ODE solver using the explicit forward Euler method
104: Level: beginner
106: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSBEULER`, `TSType`
107: M*/
108: PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts)
109: {
110: TS_Euler *euler;
112: PetscNew(&euler);
113: ts->data = (void *)euler;
115: ts->ops->setup = TSSetUp_Euler;
116: ts->ops->step = TSStep_Euler;
117: ts->ops->reset = TSReset_Euler;
118: ts->ops->destroy = TSDestroy_Euler;
119: ts->ops->setfromoptions = TSSetFromOptions_Euler;
120: ts->ops->view = TSView_Euler;
121: ts->ops->interpolate = TSInterpolate_Euler;
122: ts->ops->linearstability = TSComputeLinearStability_Euler;
123: ts->default_adapt_type = TSADAPTNONE;
124: ts->usessnes = PETSC_FALSE;
125: return 0;
126: }