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: }