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 <array>
+
+namespace specmicp {
+namespace odeint {
+
+template <int order>
+class ButcherTableau {
+public:
+ ButcherTableau(const std::array<double, order>& aa,
+ const std::array<std::array<double, order>, order>& ba,
+ const std::array<double, order+1>& ca,
+ const std::array<double, order+1>& csa
+ ):
+ m_a(aa), m_b(ba), m_c(ca), m_cs(csa)
+ {}
+ ButcherTableau(std::array<double, order>&& aa,
+ std::array<std::array<double, order>, order>&& ba,
+ std::array<double, order+1>&& ca,
+ std::array<double, order+1>&& 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<double, order> m_a;
+ std::array<std::array<double, order>, order> m_b;
+ std::array<double, order+1> m_c;
+ std::array<double, order+1> 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 <functional>
+#include <Eigen/Dense>
+
+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 <int dim>
+struct vector_trait
+{
+ using type = Eigen::Matrix<scalar_t, dim, 1>;
+};
+
+template <>
+struct vector_trait<1>
+{
+ using type = scalar_t;
+};
+
+} // end namespace internal
+
+//! The vector type
+template <int dim>
+using vector_t = typename internal::vector_trait<dim>::type;
+
+//! Type of the function f in dy/dx=f(x,y)
+template <int dim>
+using rhs_f = std::function<void (scalar_t, const vector_t<dim>&, vector_t<dim>&)>;
+
+
+} // 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 <int dim, int order>
+class EmbeddedRungeKuttaStep
+{
+public:
+ using vec_t = vector_t<dim>;
+
+ EmbeddedRungeKuttaStep(const rhs_f<dim>& rhs,
+ const ButcherTableau<order>& 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
+ ButcherTableau<order>& butcher() {return m_butcher;}
+
+private:
+ rhs_f<dim> m_rhs;
+ ButcherTableau<order> m_butcher;
+};
+
+//! Embedded Runge Kutta of fifth order
+template <int dim>
+class EmbeddedRungeKuttaStep45: EmbeddedRungeKuttaStep<dim, 5>
+{
+public:
+ using base = EmbeddedRungeKuttaStep<dim, 5>;
+ using vec_t = typename base::vec_t;
+ using base::butcher;
+ using base::get_rhs;
+
+ EmbeddedRungeKuttaStep45(const rhs_f<dim>& 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 <int dim>
+void EmbeddedRungeKuttaStep45<dim>::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);
+
+ }
+};