Homec4science
Diffusion SpecMiCP / ReactMiCP (gasphase) Browse

specmicp/src/odeint/runge_kutta_step.hpp
881a9cd5b582gasphase

Last Change
Raw File

runge_kutta_step.hpp
View Options

/*-------------------------------------------------------
- Module : odeint
- File : runge_kutta_step
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP
#define SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP
//! \file runge_kutta_step algorithm for embedded runge kutta
#include "odeint_types.hpp"
#include "butcher_tableau.hpp"
namespace specmicp {
namespace odeint {
struct EmbeddedRungeKuttaStepOptions
{
double safety;
double p_grow;
double p_shrink;
double errcon;
double max_try;
double eps;
EmbeddedRungeKuttaStepOptions():
safety(0.9), p_grow(-0.2), p_shrink(-0.25), errcon(1.89e-4), max_try(20), eps(1e-5)
{}
};
struct StepLength
{
const double not_avail = -1;
double test;
double did;
double next;
StepLength(double h_test):
test(h_test), did(not_avail), next(not_avail)
{}
void get_next() {test=next; did=not_avail; next=not_avail;}
bool is_next_available() {return next != not_avail;}
};
//! \brief Base class for computing a step of the Range Kutta algorithm
template <int dim, class butcher_tableau_t>
class EmbeddedRungeKuttaStep
{
public:
using vec_t = vector_t<dim>;
EmbeddedRungeKuttaStep(const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher):
m_rhs(rhs),
m_butcher(butcher)
{}
//! \brief compute a step
//!
//! \param y initial value
//! \param dydx initial value of the right hand site
//! \param x value of the integration variable
//! \param timestep the integration step
//! \param y_out value at the end of the step (y+timestep)
//! \param y_err estimation of the truncation error
void run_step(const vec_t& y, const vec_t& dydx, scalar_t x, scalar_t timestep, vec_t& y_out, vec_t& y_err);
protected:
//! Return the value of the right hand side at x for value y
void get_rhs(scalar_t x, const vector_t<dim>& y, vec_t& dydx) {return m_rhs(x, y, dydx);}
//! Return the Butcher tableau
butcher_tableau_t& butcher() {return m_butcher;}
EmbeddedRungeKuttaStepOptions& get_options() {return m_options;}
const EmbeddedRungeKuttaStepOptions& get_options() const {return m_options;}
private:
EmbeddedRungeKuttaStepOptions m_options;
rhs_f<dim> m_rhs;
butcher_tableau_t m_butcher;
};
//! Embedded Runge Kutta of fifth order
template <int dim, class butcher_tableau_t>
class EmbeddedRungeKuttaStep45: EmbeddedRungeKuttaStep<dim, butcher_tableau_t>
{
public:
using base = EmbeddedRungeKuttaStep<dim, butcher_tableau_t>;
using vec_t = typename base::vec_t;
using base::butcher;
using base::get_rhs;
using base::get_options;
//check that template argument are consistent
static_assert(butcher_tableau_t::RK_order == 5, "The butcher tableau must correspond to a Runge-Kutta method of order 5");
static_assert(butcher_tableau_t::RK_evaluations == 6 || butcher_tableau_t::RK_evaluations == 7, "The number of functions evaluations must be 6 or 7 for an embedded Runge-Kutta method");
EmbeddedRungeKuttaStep45(const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher):
base(rhs, butcher)
{
}
void run_step(const vec_t& y, const vec_t& dydx, scalar_t x, scalar_t timestep, vec_t& y_out, vec_t& y_err);
void rk_step(vector_t<dim> &y,
const vector_t<dim>& dydx,
double& x,
StepLength& stepl);
};
//! typedef to create a Dormant-Prince method
template <int dim>
using DormandPrinceStep = EmbeddedRungeKuttaStep45<dim, ButcherTableauDormandPrince_t>;
template <int dim>
DormandPrinceStep<dim> get_dormand_prince_step(const rhs_f<dim>& rhs) {return DormandPrinceStep<dim>(rhs, butcher_dormand_prince45);}
//! typedef to create a Cash-Karp method
template <int dim>
using CashKarpStep = EmbeddedRungeKuttaStep45<dim, ButcherTableauCashKarp_t>;
template <int dim>
CashKarpStep<dim> get_cash_karp_step(const rhs_f<dim>& rhs) {return CashKarpStep<dim>(rhs, butcher_cash_karp45);}
template <int dim, class T>
T get_rk_step(const rhs_f<dim>& rhs)
{
if (std::is_same<T, CashKarpStep<dim>>::value) return get_cash_karp_step(rhs);
else if (std::is_same<T, DormandPrinceStep<dim>>::value) return get_dormand_prince_step(rhs);
}
} // end namespace odeint
} // end namespace specmicp
#include "runge_kutta_step.inl"
#endif // SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP
c4science · Help