diff --git a/src/odeint/butcher_tableau.hpp b/src/odeint/butcher_tableau.hpp index 5693758..109ebc0 100644 --- a/src/odeint/butcher_tableau.hpp +++ b/src/odeint/butcher_tableau.hpp @@ -1,107 +1,108 @@ /*------------------------------------------------------- - Module : odeint - File : butcher_tableau.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget , 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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_BUTCHERTABLEAU_HPP #define SPECMICP_ODEINT_BUTCHERTABLEAU_HPP //! \file butcher_tableau.hpp Butcher tableau for embedded runge kutta method #include +#include "common.hpp" namespace specmicp { namespace odeint { template class ButcherTableau { public: static const int RK_order = order; static const int RK_evaluations = s+1; - ButcherTableau(const std::array& aa, - const std::array, s>& ba, - const std::array& ca, - const std::array& csa + ButcherTableau(const std::array& aa, + const std::array, s>& 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, s>&& ba, - std::array&& ca, - std::array&& csa + ButcherTableau(std::array&& aa, + std::array, s>&& 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, s> m_b; - std::array m_c; - std::array m_cs; + std::array m_a; + std::array, s> m_b; + std::array m_c; + std::array m_cs; }; const ButcherTableau<5, 5> butcher_cash_karp45({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}, {1631.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} ); using ButcherTableauCashKarp_t = ButcherTableau<5, 5>; const ButcherTableau<5, 6> butcher_dormand_prince45( {1.0/5.0, 3.0/10.0, 4.0/5.0, 8.0/9.0, 1.0, 1.0}, { {{1.0/5.0, 0, 0, 0, 0, 0}, {3.0/40.0, 9.0/40.0, 0, 0, 0, 0}, {44.0/45.0, -56.0/15.0, 32.0/9.0, 0, 0, 0}, {19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0, 0, 0}, {9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/18656.0, 0}, {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0} }}, {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0, 0.0}, {5179.0/57600.0, 0.0, 7571.0/16695.0, 393.0/640.0, -92097.0/339200.0, 187.0/2100.0, 1.0/40.0} ); using ButcherTableauDormandPrince_t = ButcherTableau<5, 6>; } // 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 index 215b77a..6aa7201 100644 --- a/src/odeint/odeint_types.hpp +++ b/src/odeint/odeint_types.hpp @@ -1,100 +1,98 @@ /*------------------------------------------------------- - Module : odeint - File : odeint_types.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget , 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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_ODEINTTYPES_HPP #define SPECMICP_ODEINT_ODEINTTYPES_HPP //! \file odeint_types.hpp common types for odeint +#include "common.hpp" #include -#include +#include namespace specmicp { namespace odeint { -//! A scalar -using scalar_t = double; - namespace internal { //! \brief implementation for the vector type (which may not be a vector for dim=1) template struct vector_trait { using type = Eigen::Matrix; }; //! Specialisation for dimension 1 (scalar) template <> struct vector_trait<1> { using type = scalar_t; }; //! Specialization for dimension 2 (fixed-size vector) template <> struct vector_trait<2> { using type = Eigen::Matrix; }; //! Specialization for dimension 3 (fixed-size vector) template <> struct vector_trait<3> { using type = Eigen::Matrix; }; } // end namespace internal //! The vector type template using vector_t = typename internal::vector_trait::type; //! Return the maximum coefficient of the vector template scalar_t max_coeff(const vector_t& vec) {return vec.maxCoeff();} template <> scalar_t max_coeff<1>(const vector_t<1>& vec) {return vec;} //! 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 index d64a29b..ff49f64 100644 --- a/src/odeint/runge_kutta_step.hpp +++ b/src/odeint/runge_kutta_step.hpp @@ -1,168 +1,178 @@ /*------------------------------------------------------- - Module : odeint - File : runge_kutta_step - Author : Fabien Georget Copyright (c) 2014, Fabien Georget , 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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" namespace specmicp { namespace odeint { struct EmbeddedRungeKuttaStepOptions { - double safety; - double p_grow; - double p_shrink; - double errcon; - double max_try; - double eps; + scalar_t safety; + scalar_t p_grow; + scalar_t p_shrink; + scalar_t errcon; + scalar_t max_try; + scalar_t 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; + static constexpr scalar_t not_available = -1.0; + scalar_t test; + scalar_t did; + scalar_t next; - StepLength(double h_test): - test(h_test), did(not_avail), next(not_avail) + StepLength(scalar_t h_test): + test(h_test), did(not_available), next(not_available) {} - void get_next() {test=next; did=not_avail; next=not_avail;} - bool is_next_available() {return next != not_avail;} + void get_next() {test=next; did=not_available; next=not_available;} + bool is_next_available() {return next != not_available;} }; //! \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 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); + 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);} + void get_rhs( + scalar_t x, + const vector_t& 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 m_rhs; butcher_tableau_t 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; 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& 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 &y, const vector_t& dydx, double& x, StepLength& stepl); }; //! typedef to create a Dormant-Prince method template using DormandPrinceStep = EmbeddedRungeKuttaStep45; template DormandPrinceStep get_dormand_prince_step(const rhs_f& rhs) {return DormandPrinceStep(rhs, butcher_dormand_prince45);} //! typedef to create a Cash-Karp method template using CashKarpStep = EmbeddedRungeKuttaStep45; template CashKarpStep get_cash_karp_step(const rhs_f& rhs) {return CashKarpStep(rhs, butcher_cash_karp45);} template T get_rk_step(const rhs_f& rhs) { if (std::is_same>::value) return get_cash_karp_step(rhs); else if (std::is_same>::value) return get_dormand_prince_step(rhs); } } // end namespace odeint } // end namespace specmicp #include "runge_kutta_step.inl" #endif // SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP diff --git a/tests/specmicp/kinetics.cpp b/tests/specmicp/kinetics.cpp index 47540b5..2e31c17 100644 --- a/tests/specmicp/kinetics.cpp +++ b/tests/specmicp/kinetics.cpp @@ -1,217 +1,221 @@ /*------------------------------------------------------- - Module : tests/ - File : kinetics.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget , 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: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * 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. * Neither the name of the Princeton University 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 OWNER 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. ---------------------------------------------------------*/ - +#include "common.hpp" #include "specmicp/reduced_system_solver.hpp" #include "database/database.hpp" #include #include #include "utils/log.hpp" #include "odeint/runge_kutta_step.hpp" #include "specmicp/equilibrium_data.hpp" const double Mc3s = 228.3; const double Vmc3s = 73e-6; +using specmicp::index_t; +using specmicp::scalar_t; + class RateComputer { EIGEN_MAKE_ALIGNED_OPERATOR_NEW public: static constexpr double m_c3s0 = 500; const double R0 = 2.5e-6; RateComputer(): - m_database("data/cemdata_specmicp.js") + m_database("../data/cemdata_specmicp.js") { std::shared_ptr data = m_database.get_database(); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); m_database.swap_components(swapping); - std::vector to_keep = {0, 1, + std::vector to_keep = {0, 1, m_database.safe_label_to_id("Ca[2+]", m_database.get_database()->labels_basis), m_database.safe_label_to_id("SiO(OH)3[-]", m_database.get_database()->labels_basis) }; m_database.keep_only_components(to_keep); m_equilvar.resize(data->nb_component+data->nb_mineral); m_equilvar(0) = 1.0; m_equilvar.block(1, 0, data->nb_component-1, 1).setConstant(-7); m_equilvar.block(data->nb_component, 0, data->nb_mineral, 1).setConstant(0.); m_totconc = Eigen::VectorXd(data->nb_component); m_totconc0 = Eigen::VectorXd(data->nb_component); m_totconc << 1.0/data->molar_mass_basis(0), 1e-7, 1e-6, 2e-6; m_totconc0 << 1.0/data->molar_mass_basis(0), 1e-7, 1e-6, 2e-6; id_c3s = m_database.mineral_kinetic_label_to_id("C3S"); id_ch = m_database.mineral_label_to_id("Portlandite"); id_csht = m_database.mineral_label_to_id("CSH,tobermorite"); id_cshj = m_database.mineral_label_to_id("CSH,jennite"); id_ca = m_database.component_label_to_id("Ca[2+]"); n_particle = 3*Vmc3s*m_c3s0/(Mc3s*4*M_PI*std::pow(R0, 3.0)); } void update_totconc(double mc3s) { std::shared_ptr data = m_database.get_database(); double xi = (m_c3s0 - mc3s)/Mc3s; //std::cout << "xi : " << xi << std::endl; //std::cout <nu_mineral_kinetic.row(id_c3s).transpose() << std::endl; m_totconc.noalias() = m_totconc0 + xi*data->nu_mineral_kinetic.row(id_c3s).transpose(); } double compute_logiap() { return m_current_solution.logIAP_kinetic(id_c3s); } double compute_rate(double mc3s, double logiap) { double radius = std::pow((Vmc3s*mc3s*3)/(Mc3s*n_particle*4*M_PI) ,1.0/3.0); double surf = (3*Vmc3s/(Mc3s*radius)); //std::cout << n_particle << "\t" << radius << "\t" << surf << std::endl; return std::min((surf*mc3s)*Mc3s*1e-6*(2.1*logiap+35.7), -1e-6); } void get_rate(double _, double mc3s, double& rate) { //std::cout << m_totconc << std::endl; update_totconc(mc3s); //std::cout << m_totconc << std::endl; compute_equilibrium(); rate = compute_rate(mc3s, compute_logiap()); } double get_ca() { return m_current_solution.molality_component(id_ca); } double get_ch() { return m_current_solution.mass_mineral(id_ch); } double get_csht() { return m_current_solution.mass_mineral(id_csht); } double get_cshj() { return m_current_solution.mass_mineral(id_cshj); } double get_pH() { return m_current_solution.pH(); } void compute_equilibrium() { specmicp::ReducedSystemSolver solver(m_database.get_database(), m_totconc); + solver.get_options().solver_options.maxstep = 50; specmicp::micpsolver::MiCPPerformance perf = solver.solve(m_equilvar); if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized) { std::cout << "Failed to solve the system ! Err code " << (int) perf.return_code << std::endl; } m_current_solution = solver.get_solution(m_equilvar); } Eigen::VectorXd get_vars() {return m_equilvar;} private: Eigen::VectorXd m_equilvar; Eigen::VectorXd m_totconc0; Eigen::VectorXd m_totconc; specmicp::database::Database m_database; specmicp::EquilibriumState m_current_solution; int id_c3s; int id_ch; int id_csht; int id_cshj; int id_ca; double n_particle; }; int main() { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; RateComputer rate; double y = RateComputer::m_c3s0-0.05; double dydx; rate.get_rate(0, y, dydx); std::cout << std::setprecision(6); specmicp::odeint::rhs_f<1> rhs = std::bind(std::mem_fn(&RateComputer::get_rate), rate, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3); specmicp::odeint::CashKarpStep<1> driver = specmicp::odeint::get_cash_karp_step<1>(rhs); //specmicp::odeint::DormandPrinceStep<1> driver = specmicp::odeint::get_dormand_prince_step<1>(rhs); specmicp::odeint::StepLength stepl(1.0); double x = 0.0; int cnt = 0; while (y > 1e-2) //while (cnt < 10) { driver.rk_step(y, dydx, x, stepl); cnt +=1; stepl.get_next(); driver.get_rhs(x, y, dydx); rate.update_totconc(y); // shouldn't be needed rate.compute_equilibrium(); std::cout << x << "\t" << y << "\t" << rate.compute_logiap() << "\t" << rate.get_pH() << "\t" << rate.get_ca() << "\t" << rate.get_ch() << "\t" << rate.get_csht() << "\t" << rate.get_cshj() << std::endl; } // Eigen::VectorXd equil = rate.get_vars() ; // y+= 0.00001; // rate.get_rate(0, y, dydx); // Eigen::VectorXd equil2 = rate.get_vars(); // std::cout << equil2 - equil << std::endl; //std::cout << "logiap : " << rate.compute_logiap() << " - rate : " << rate.compute_rate(rate.compute_logiap()) << std::endl; return EXIT_SUCCESS; }