Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/src/specmicp_common/odeint/runge_kutta_step.hpp
d2545c482392master

Last Change
Raw File

runge_kutta_step.hpp
View Options

/* =============================================================================
Copyright (c) 2014 - 2016
F. Georget <fabieng@princeton.edu> Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
3. Neither the name of the copyright holder nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
============================================================================= */
#ifndef SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP
#define SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP
//! \file runge_kutta_step.hpp
//! \brief Stepping algorithm for the embedded runge kutta ODE integration method
#include "odeint_types.hpp"
#include "butcher_tableau.hpp"
#include "runge_kutta_step_structs.hpp"
#include "specmicp_common/options_handler.hpp"
namespace specmicp {
namespace odeint {
//! \brief Struct to store the different timesteps
struct StepLength
{
static constexpr scalar_t not_available = -1.0; //!< Value is not available
scalar_t test; //!< Proposed timestep
scalar_t did; //!< Last timestep
scalar_t next; //!< Next timestep
scalar_t total; //!< Total time
//! \brief Constructor, give the finish line
StepLength(scalar_t h_test):
test(h_test), did(not_available), next(not_available), total(0.0)
{}
//! \brief Constructor with finish line
StepLength(scalar_t h_test, scalar_t h_total):
test(h_test), did(not_available), next(not_available), total(h_total)
{}
//! \brief Reset for a new timestep
void get_next() {
test=next;
did=not_available;
next=not_available;
}
//! \brief Check if next timestep is available
bool is_next_available() {return next != not_available;}
//! \brief Return the
bool is_total_time_fixed() {return total > 0.0;}
//! \brief Return the remaining time
scalar_t remaining(scalar_t t) {return (total-t);}
};
//! \brief Base class for computing a step of the Range Kutta algorithm
template <int dim, class butcher_tableau_t>
class EmbeddedRungeKuttaStep : public OptionsHandler<EmbeddedRungeKuttaStepOptions>
{
public:
using vec_t = vector_t<dim>; //!< Type of a vector
//! \brief Constructor with default options
EmbeddedRungeKuttaStep(
const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher):
m_rhs(rhs),
m_butcher(butcher)
{}
//! \brief Constructor with custom options
EmbeddedRungeKuttaStep(
const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher,
EmbeddedRungeKuttaStepOptions& options
):
OptionsHandler(options),
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;}
private:
rhs_f<dim> m_rhs; //!< Function returning the right hand side
butcher_tableau_t m_butcher; //! the Butcher tableeau
};
//! 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>; //!< Type of the base class
using vec_t = typename base::vec_t; //!< Type of a vector
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");
//! \brief Constructor with default options
EmbeddedRungeKuttaStep45(const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher):
base(rhs, butcher)
{
}
//! \brief Constructor with custom options
EmbeddedRungeKuttaStep45(
const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher,
EmbeddedRungeKuttaStepOptions options
):
base(rhs, butcher, options)
{
}
//! \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);
//! \brief Run a runge kutta step
void rk_step(vector_t<dim> &y,
const vector_t<dim>& dydx,
double& x,
StepLength& stepl);
};
//! \typedef DormandPrinceStep
//! \brief Create a Dormant-Prince method
template <int dim>
using DormandPrinceStep = EmbeddedRungeKuttaStep45<dim, ButcherTableauDormandPrince_t>;
//! \brief Stepper for the Dormant-Prince method
template <int dim>
DormandPrinceStep<dim> get_dormand_prince_step(const rhs_f<dim>& rhs) {return DormandPrinceStep<dim>(rhs, butcher_dormand_prince45);}
//! \typedef CashKarpStep
//! \brief Create a Cash-Karp method
template <int dim>
using CashKarpStep = EmbeddedRungeKuttaStep45<dim, ButcherTableauCashKarp_t>;
//! \brief Stepper for the Cash-Karp method
template <int dim>
CashKarpStep<dim> get_cash_karp_step(const rhs_f<dim>& rhs) {return CashKarpStep<dim>(rhs, butcher_cash_karp45);}
//! \brief Stepper for Runge-Kutta algorithm
//!
//! Either Dormant-Prince or Cash-Karp
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