Homec4science
Diffusion SpecMiCP / ReactMiCP (chloride) Browse

specmicp/src/odeint/runge_kutta_step.hpp
b174f1ddf352chloride

Last Change
Raw File

runge_kutta_step.hpp
View Options

/*-------------------------------------------------------------------------------
Copyright (c) 2014,2015 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 algorithm for embedded runge kutta
#include "odeint_types.hpp"
#include "butcher_tableau.hpp"
#include "runge_kutta_step_structs.hpp"
#include "utils/options_handler.hpp"
namespace specmicp {
namespace odeint {
struct StepLength
{
static constexpr scalar_t not_available = -1.0;
scalar_t test;
scalar_t did;
scalar_t next;
scalar_t total;
StepLength(scalar_t h_test):
test(h_test), did(not_available), next(not_available), total(0.0)
{}
StepLength(scalar_t h_test, scalar_t h_total):
test(h_test), did(not_available), next(not_available), total(h_total)
{}
void get_next() {
test=next;
did=not_available;
next=not_available;
}
bool is_next_available() {return next != not_available;}
bool is_total_time_fixed() {return total > 0.0;}
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>;
EmbeddedRungeKuttaStep(
const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher):
m_rhs(rhs),
m_butcher(butcher)
{}
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;
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)
{
}
EmbeddedRungeKuttaStep45(
const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher,
EmbeddedRungeKuttaStepOptions options
):
base(rhs, butcher, options)
{
}
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