diff --git a/CMakeLists.txt b/CMakeLists.txt index 7f0109d..b060ec9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,189 +1,212 @@ project(specmicp) cmake_minimum_required(VERSION 2.8) ####################### External Package ##################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package find_package(CxxTest) # Necessary if you want the tests ######################## Compilation flags ################################# if(UNIX) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11") - SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg -Wall -pedantic") + SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -pedantic") SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -DNDEBUG") message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}") else() message(WARNING "not tested !") endif() ################ Directories ##################### set(PROJECT_TEST_DIR ${PROJECT_SOURCE_DIR}/tests) set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # include include_directories(${EIGEN3_INCLUDE_DIR}) include_directories(${PROJECT_SOURCE_DIR}) ################ Modules ########################## # -------- MiCPSolver ---------------- set(MICPSOLVER_DIR ${PROJECT_SOURCE_DIR}/micpsolver) add_custom_target(specmisolver_inc SOURCES ${MICPSOLVER_DIR}/micpsolver_base.hpp ${MICPSOLVER_DIR}/micpsolver_base.inl ${MICPSOLVER_DIR}/micpsolver.hpp ${MICPSOLVER_DIR}/micpsolver.inl ${MICPSOLVER_DIR}/micpsolver_min.hpp ${MICPSOLVER_DIR}/micpsolver_min.inl ${MICPSOLVER_DIR}/micpsolver_structs.hpp ${MICPSOLVER_DIR}/micpprog.hpp ${MICPSOLVER_DIR}/ncp_function.hpp ${MICPSOLVER_DIR}/estimate_cond_number.hpp ) # ------- SpecMiCP --------------- set(SPECMICP_DIR ${PROJECT_SOURCE_DIR}/specmicp) add_custom_target(specmic_inc SOURCES #${SPECMICP_DIR}/.hpp ${SPECMICP_DIR}/thermodata.hpp ) set(SPECMICP_LIBFILE ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp ) add_library(specmicp ${SPECMICP_LIBFILE}) # -------- Database ---------------- set(DATABASE_DIR ${PROJECT_SOURCE_DIR}/database) set(DATABASE_LIBFILE ${DATABASE_DIR}/database.cpp ${DATABASE_DIR}/reader.cpp ${DATABASE_DIR}/selector.cpp ${DATABASE_DIR}/mineral_selector.cpp ${DATABASE_DIR}/canonicalizer.cpp ${DATABASE_DIR}/switch_basis.cpp ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/data_container.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) +# --------- ODE int ------------------ + +set(ODEINT_DIR ${PROJECT_SOURCE_DIR}/odeint) + + +add_custom_target(odeint_incl SOURCES + ${ODEINT_DIR}/odeint_types.hpp + ${ODEINT_DIR}/butcher_tableau.hpp + ${ODEINT_DIR}/runge_kutta_step.hpp + ${ODEINT_DIR}/runge_kutta_step.inl +) + # ------- Utils ------------------ set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.hpp) # ------- Data ------------------- set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) set(DATABASE_LIST ${DATA_DIR}/specmicp_database.js ${DATA_DIR}/cemdata_specmicp.js ) add_custom_target(data SOURCES ${DATABASE_LIST}) file(INSTALL ${DATABASE_LIST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data/) # ------ Documentation ----------- # add a target to generate API documentation with Doxygen find_package(Doxygen) if(DOXYGEN_FOUND) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/doc/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile @ONLY) add_custom_target(doc ${DOXYGEN_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "Generating API documentation with Doxygen" VERBATIM ) add_custom_target(citations SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib) file(INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/doc/) endif(DOXYGEN_FOUND) # --------- Test ---------------------- # CXX Test # ======== if(CXXTEST_FOUND) set(CXXTEST_USE_PYTHON TRUE) include_directories(${CXXTEST_INCLUDE_DIR}) enable_testing() # MiCP Solver # ------------ CXXTEST_ADD_TEST(test_cond_number test_cond_number.cpp ${PROJECT_TEST_DIR}/micpsolver/test_cond_number.hpp) CXXTEST_ADD_TEST(test_ncp_funtion test_ncp_function.cpp ${PROJECT_TEST_DIR}/micpsolver/test_ncp_function.hpp) CXXTEST_ADD_TEST(test_micpsolver test_micpsolver.cpp ${PROJECT_TEST_DIR}/micpsolver/test_micpsolver.hpp) # Database # -------- CXXTEST_ADD_TEST(test_data_reader test_data_reader.cpp ${PROJECT_TEST_DIR}/database/test_data_reader.hpp) target_link_libraries(test_data_reader specmicp_database jsoncpp) CXXTEST_ADD_TEST(test_data_selector test_data_selector.cpp ${PROJECT_TEST_DIR}/database/test_data_selector.hpp) target_link_libraries(test_data_selector specmicp_database jsoncpp) + # Odeint + # ------ + + CXXTEST_ADD_TEST(test_butchertableau + test_butchertableau.cpp + ${PROJECT_TEST_DIR}/odeint/test_butchertableau.hpp) + + + CXXTEST_ADD_TEST(test_embeddedrungekutta + test_embeddedrungekutta.cpp + ${PROJECT_TEST_DIR}/odeint/test_embeddedrungekutta.hpp) endif() # Test 'Perso' # ============ add_executable(thermocarbo ${PROJECT_TEST_DIR}/specmicp/thermocarbo.cpp) target_link_libraries(thermocarbo specmicp specmicp_database jsoncpp) add_executable(carboalu ${PROJECT_TEST_DIR}/specmicp/carboalu.cpp) target_link_libraries(carboalu specmicp specmicp_database jsoncpp) add_executable(alkaliactivated ${PROJECT_TEST_DIR}/specmicp/alkaliactivated.cpp) target_link_libraries(alkaliactivated specmicp specmicp_database jsoncpp) add_executable(rate_nicoleau ${PROJECT_TEST_DIR}/specmicp/rate_nicoleau.cpp) target_link_libraries(rate_nicoleau specmicp specmicp_database jsoncpp) set(DATA_TEST ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3sm.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st1.csv ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st2.csv ) add_custom_target(data_nicoleau SOURCES ${DATA_TEST}) file(INSTALL ${DATA_TEST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data_test/) diff --git a/src/odeint/butcher_tableau.hpp b/src/odeint/butcher_tableau.hpp new file mode 100644 index 0000000..abe137a --- /dev/null +++ b/src/odeint/butcher_tableau.hpp @@ -0,0 +1,77 @@ +/*------------------------------------------------------- + + - Module : odeint + - File : butcher_tableau.hpp + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#ifndef SPECMICP_ODEINT_BUTCHERTABLEAU_HPP +#define SPECMICP_ODEINT_BUTCHERTABLEAU_HPP + +//! \file butcher_tableau.hpp Butcher tableau for embedded runge kutta method + +#include + +namespace specmicp { +namespace odeint { + +template +class ButcherTableau { +public: + ButcherTableau(const std::array& aa, + const std::array, order>& ba, + const std::array& ca, + const std::array& csa + ): + m_a(aa), m_b(ba), m_c(ca), m_cs(csa) + {} + ButcherTableau(std::array&& aa, + std::array, order>&& ba, + std::array&& ca, + std::array&& csa + ): + m_a(aa), m_b(ba), m_c(ca), m_cs(csa) + {} + double a(int i) const {return m_a[i-2];} + double b(int i, int j) const {return m_b[i-2][j-1];} + double c(int i) const {return m_c[i-1];} + double cs(int i) const {return m_cs[i-1];} + +private: + std::array m_a; + std::array, order> m_b; + std::array m_c; + std::array m_cs; +}; + +const ButcherTableau<5> CashKarp45({1.0/5.0, 3.0/10.0, 3.0/5.0, 1.0, 7.0/8.0}, + { {{1.0/5.0, 0, 0, 0, 0}, + {3.0/40.0, 9.0/40.0, 0, 0, 0}, + {3.0/10.0, -9.0/10.0, 6.0/5.0, 0, 0}, + {-11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0}, + {1621.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0} + }}, + {37.0/378.0, 0.0, 250.0/621.0, 125.0/594.0, 0.0, 512.0/1771.0}, + {2825.0/27648.0, 0.0, 18575.0/48384.0, 13525.0/55296.0, 277.0/14336.0, 1.0/4.0} + ); + + +const ButcherTableau<5> CashKarp45({1.0/5.0, 3.0/10.0, 4.0/5.0, 1.0, 7.0/8.0}, + { {{1.0/5.0, 0, 0, 0, 0}, + {3.0/40.0, 9.0/40.0, 0, 0, 0}, + {3.0/10.0, -9.0/10.0, 6.0/5.0, 0, 0}, + {-11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0}, + {1621.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0} + }}, + {37.0/378.0, 0.0, 250.0/621.0, 125.0/594.0, 0.0, 512.0/1771.0}, + {2825.0/27648.0, 0.0, 18575.0/48384.0, 13525.0/55296.0, 277.0/14336.0, 1.0/4.0} + ); + +} // end namespace odeint +} // end namespace specmicp + + +#endif // SPECMICP_ODEINT_BUTCHERTABLEAU_HPP diff --git a/src/odeint/odeint_types.hpp b/src/odeint/odeint_types.hpp new file mode 100644 index 0000000..3c06086 --- /dev/null +++ b/src/odeint/odeint_types.hpp @@ -0,0 +1,54 @@ +/*------------------------------------------------------- + + - Module : odeint + - File : odeint_types.hpp + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#ifndef SPECMICP_ODEINT_ODEINTTYPES_HPP +#define SPECMICP_ODEINT_ODEINTTYPES_HPP + +//! \file odeint_types.hpp common types for odeint + +#include +#include + +namespace specmicp { +namespace odeint { + +//! A scalar +using scalar_t = double; + +namespace internal { + +//! implementation for the vector type (which may not be a vector for dim=1) +template +struct vector_trait +{ + using type = Eigen::Matrix; +}; + +template <> +struct vector_trait<1> +{ + using type = scalar_t; +}; + +} // end namespace internal + +//! The vector type +template +using vector_t = typename internal::vector_trait::type; + +//! Type of the function f in dy/dx=f(x,y) +template +using rhs_f = std::function&, vector_t&)>; + + +} // end namespace odeint +} // end namespace specmicp + +#endif // SPECMICP_ODEINT_ODEINTTYPES_HPP diff --git a/src/odeint/runge_kutta_step.hpp b/src/odeint/runge_kutta_step.hpp new file mode 100644 index 0000000..284389d --- /dev/null +++ b/src/odeint/runge_kutta_step.hpp @@ -0,0 +1,80 @@ +/*------------------------------------------------------- + + - 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 { + +//! \brief Base class for computing a step of the Range Kutta algorithm +template +class EmbeddedRungeKuttaStep +{ +public: + using vec_t = vector_t; + + EmbeddedRungeKuttaStep(const rhs_f& rhs, + const ButcherTableau& 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& y, vec_t& dydx) {return m_rhs(x, y, dydx);} + //! Return the Butcher tableau + ButcherTableau& butcher() {return m_butcher;} + +private: + rhs_f m_rhs; + ButcherTableau m_butcher; +}; + +//! Embedded Runge Kutta of fifth order +template +class EmbeddedRungeKuttaStep45: EmbeddedRungeKuttaStep +{ +public: + using base = EmbeddedRungeKuttaStep; + using vec_t = typename base::vec_t; + using base::butcher; + using base::get_rhs; + + EmbeddedRungeKuttaStep45(const rhs_f& rhs, + const ButcherTableau<5>& 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); +}; + +} // end namespace odeint +} // end namespace specmicp + +#include "runge_kutta_step.inl" + +#endif // SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP + diff --git a/src/odeint/runge_kutta_step.inl b/src/odeint/runge_kutta_step.inl new file mode 100644 index 0000000..bcf9671 --- /dev/null +++ b/src/odeint/runge_kutta_step.inl @@ -0,0 +1,50 @@ +/*------------------------------------------------------- + + - Module : odeint + - File : runge_kutta_step.cpp + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#include "runge_kutta_step.hpp" // syntaxic coloration only + +namespace specmicp { +namespace odeint { + +template +void EmbeddedRungeKuttaStep45::run_step(const vec_t& y, + const vec_t& dydx, + scalar_t x, + scalar_t timestep, + vec_t& y_out, + vec_t& y_err) +{ + vec_t ak2, ak3, ak4, ak5, ak6; + vec_t ytemp; + + ytemp = y + butcher().b(2, 1)*timestep*dydx; + get_rhs(x + butcher().a(2)*timestep, ytemp, ak2); + ytemp = y + timestep*(butcher().b(3,1)*dydx + butcher().b(3, 2)*ak2); + get_rhs(x + butcher().a(3)*timestep, ytemp, ak3); + ytemp = y + timestep*(butcher().b(4,1)*dydx + butcher().b(4, 2)*ak2 + butcher().b(4, 2)*ak3); + get_rhs(x + butcher().a(4)*timestep, ytemp, ak4); + ytemp = y + timestep*(butcher().b(5,1)*dydx + butcher().b(5, 2)*ak2 + + butcher().b(5, 3)*ak3 + butcher().b(5, 4)*ak4); + get_rhs(x + butcher().a(5)*timestep, ytemp, ak5); + ytemp = y + timestep*(butcher().b(6,1)*dydx + butcher().b(6, 2)*ak2 + + butcher().b(6, 3)*ak3 + butcher().b(6, 4)*ak4 + + butcher().b(6, 5)*ak5); + get_rhs(x + butcher().a(6)*timestep, ytemp, ak6); + + y_out = y + timestep*(butcher().c(1)*dydx + butcher().c(2)*ak2 + butcher().c(3)*ak3 + + butcher().c(4)*ak4 + butcher().c(5)*ak5 + butcher().c(6)*ak6); + y_err = timestep * ( (butcher().c(1) - butcher().cs(1))*dydx + (butcher().c(2) - butcher().cs(2))*ak2 + + (butcher().c(3) - butcher().cs(3))*ak3 + (butcher().c(4) - butcher().cs(4))*ak4 + + (butcher().c(5) - butcher().cs(5))*ak5 + (butcher().c(6) - butcher().cs(6))*ak6 + ); +} + +} // end namespace odeint +} // end namespace specmicp diff --git a/tests/odeint/test_butchertableau.hpp b/tests/odeint/test_butchertableau.hpp new file mode 100644 index 0000000..814189e --- /dev/null +++ b/tests/odeint/test_butchertableau.hpp @@ -0,0 +1,26 @@ +/*------------------------------------------------------- + + - Module : tests/odeint + - File : test_butchertableau.hpp + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#include "cxxtest/TestSuite.h" + +#include "odeint/butcher_tableau.hpp" + +class TestButcherTableau : public CxxTest::TestSuite +{ +public: + + void test_cashkarp() + { + TS_ASSERT_EQUALS(specmicp::odeint::CashKarp45.a(2), 1.0/5.0); + TS_ASSERT_EQUALS(specmicp::odeint::CashKarp45.b(6,3), 575.0/13824.0); + TS_ASSERT_EQUALS(specmicp::odeint::CashKarp45.c(4), 125.0/594.0); + TS_ASSERT_EQUALS(specmicp::odeint::CashKarp45.cs(4), 13525.0/55296.0); + } +}; diff --git a/tests/odeint/test_embeddedrungekutta.hpp b/tests/odeint/test_embeddedrungekutta.hpp new file mode 100644 index 0000000..41cc0f1 --- /dev/null +++ b/tests/odeint/test_embeddedrungekutta.hpp @@ -0,0 +1,40 @@ +/*------------------------------------------------------- + + - Module : tests/odeint + - File : test_embeddedrungekutta.hpp + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#include "cxxtest/TestSuite.h" +#include "odeint/runge_kutta_step.hpp" + + +using namespace specmicp::odeint; + +void rhs1(double x, const double& yin, double& yout) +{ + yout = -2*yin; +} + +class TestEmbeddedRungeKutta : public CxxTest::TestSuite +{ +public: + + void test_simpletest() + { + EmbeddedRungeKuttaStep45<1> algo(rhs1, CashKarp45); + double yout, yerr; + double tf = 0.01; + algo.run_step(1, -2, 0, tf, yout, yerr); + TS_ASSERT_DELTA(yout, std::exp(-2*tf), 1e-3); + + + tf = 0.05; + algo.run_step(1, -2, 0, tf, yout, yerr); + TS_ASSERT_DELTA(yout, std::exp(-2*tf), 1e-2); + + } +};