ode-toolbox

ODE integration tools
git clone https://git.0xfab.ch/ode-toolbox.git
Log | Files | Refs | README | LICENSE

BDF.h (2231B)

      1 /* File:   BDF.h */
      2 /* Date:   Mon Feb  8 11:44:23 2016 */
      3 /* Author: Fabian Wermelinger */
      4 /* Tag:    Backward differentiation formula using SUNDIALS library */
      5 /* Copyright 2016 ETH Zurich. All Rights Reserved. */
      6 #ifndef BDF_H_YF523CRX
      7 #define BDF_H_YF523CRX
      8 
      9 #ifdef _USE_SUNDIALS_
     10 #include <ODETB/TimeStepper/StepperBase.h>
     11 
     12 #include <cvode/cvode.h>
     13 #include <cvode/cvode_dense.h>
     14 #include <nvector/nvector_serial.h>
     15 
     16 #include <cassert>
     17 
     18 template <typename Tinput, typename Trhs>
     19 class BDF : public StepperBase<Tinput,Trhs>
     20 {
     21     Real m_a;
     22     void* m_cvode_mem;
     23     N_Vector m_y;
     24 
     25     void _init_sundials()
     26     {
     27         m_cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
     28         m_y = N_VMake_Serial(Tinput::SIZE, m_U.data());
     29         CVodeInit(m_cvode_mem, BDF::f, 0.0, m_y);
     30         CVodeSStolerances(m_cvode_mem, 1.0e-6, 1.0e-6);
     31         CVDense(m_cvode_mem, Trhs::SIZE);
     32         CVodeSetUserData(m_cvode_mem, static_cast<void*>(&m_a));
     33     }
     34 
     35 protected:
     36     using StepperBase<Tinput,Trhs>::m_settings;
     37     using StepperBase<Tinput,Trhs>::m_U;
     38     using StepperBase<Tinput,Trhs>::m_kernel;
     39     Trhs m_rhs;
     40 
     41 public:
     42     BDF(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
     43         StepperBase<Tinput,Trhs>(U,S,kern), m_rhs(U.size())
     44     {
     45         _init_sundials();
     46     }
     47     ~BDF() { CVodeFree(&m_cvode_mem); }
     48 
     49     virtual void step(void const* const data=nullptr)
     50     {
     51         /* TODO: (fabianw; Fri Feb 12 13:26:02 2016) this does not work yet.
     52          * Do not compile with Sundials support */
     53         /* m_a = a; */
     54         Real tret; // interface dummy
     55         CVode(m_cvode_mem, t+dt, m_y, &tret, CV_ONE_STEP);
     56         t = tret;
     57 
     58         if (m_settings.step % m_settings.reportGranularity == 0)
     59             std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
     60     }
     61 
     62     virtual inline void setIC(Real const ic, const Real t0=0)
     63     {
     64         m_U = ic;
     65         CVodeReInit(m_cvode_mem, t0, m_y);
     66     }
     67 
     68     static inline int f(Real t, N_Vector y, N_Vector ydot, void *a)
     69     {
     70         Tkernel kernel(Trhs::SIZE);
     71         kernel.compute(NV_DATA_S(y), NV_DATA_S(ydot), *static_cast<Real*>(a), t);
     72         return 0;
     73     }
     74 };
     75 #endif /* _USE_SUNDIALS_ */
     76 
     77 #endif /* BDF_H_YF523CRX */