Actual source code: euler.c
1: #define PETSCTS_DLL
3: /*
4: Code for Timestepping with explicit Euler.
5: */
6: #include include/private/tsimpl.h
8: typedef struct {
9: Vec update; /* work vector where F(t[i],u[i]) is stored */
10: } TS_Euler;
14: static PetscErrorCode TSSetUp_Euler(TS ts)
15: {
16: TS_Euler *euler = (TS_Euler*)ts->data;
20: VecDuplicate(ts->vec_sol,&euler->update);
21: return(0);
22: }
26: static PetscErrorCode TSStep_Euler(TS ts,PetscInt *steps,PetscReal *ptime)
27: {
28: TS_Euler *euler = (TS_Euler*)ts->data;
29: Vec sol = ts->vec_sol,update = euler->update;
31: PetscInt i,max_steps = ts->max_steps;
32: PetscScalar dt = ts->time_step;
33:
35: *steps = -ts->steps;
36: TSMonitor(ts,ts->steps,ts->ptime,sol);
38: for (i=0; i<max_steps; i++) {
39: ts->ptime += ts->time_step;
40: TSComputeRHSFunction(ts,ts->ptime,sol,update);
41: VecAXPY(sol,dt,update);
42: ts->steps++;
43: TSMonitor(ts,ts->steps,ts->ptime,sol);
44: if (ts->ptime > ts->max_time) break;
45: }
47: *steps += ts->steps;
48: *ptime = ts->ptime;
49: return(0);
50: }
51: /*------------------------------------------------------------*/
54: static PetscErrorCode TSDestroy_Euler(TS ts)
55: {
56: TS_Euler *euler = (TS_Euler*)ts->data;
60: if (euler->update) {VecDestroy(euler->update);}
61: PetscFree(euler);
62: return(0);
63: }
64: /*------------------------------------------------------------*/
68: static PetscErrorCode TSSetFromOptions_Euler(TS ts)
69: {
71: return(0);
72: }
76: static PetscErrorCode TSView_Euler(TS ts,PetscViewer viewer)
77: {
79: return(0);
80: }
82: /* ------------------------------------------------------------ */
84: /*MC
85: TS_EULER - ODE solver using the explicit forward Euler method
87: Level: beginner
89: .seealso: TSCreate(), TS, TSSetType(), TS_BEULER
91: M*/
95: PetscErrorCode TSCreate_Euler(TS ts)
96: {
97: TS_Euler *euler;
101: ts->ops->setup = TSSetUp_Euler;
102: ts->ops->step = TSStep_Euler;
103: ts->ops->destroy = TSDestroy_Euler;
104: ts->ops->setfromoptions = TSSetFromOptions_Euler;
105: ts->ops->view = TSView_Euler;
107: PetscNew(TS_Euler,&euler);
108: ts->data = (void*)euler;
110: return(0);
111: }