diff --git a/src/specmicp/equilibrium_data.hpp b/src/specmicp/equilibrium_data.hpp index 6057b6d..2b1240b 100644 --- a/src/specmicp/equilibrium_data.hpp +++ b/src/specmicp/equilibrium_data.hpp @@ -1,274 +1,274 @@ /*------------------------------------------------------- - Module : specmicp - File : equilibrium_data.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_SPECMICP_EQUILIBRIUMDATA_HPP #define SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP //! \file equilibrium_data.hpp store results from equilibrium computation #include "common.hpp" #include "database.hpp" #include "simulation_box.hpp" #define POW10(m) std::pow(10.0, m) namespace specmicp { //! \brief Store the equilibrium state of a simulation //! //! This class is used to store the output of an equilibrium solver class EquilibriumState: public SimulationBox { public: //! \brief Blank state //! //! Must be initialized after, should not be used except by internals EquilibriumState(): m_has_been_properly_initialised(false) {} //! \brief Copy constructor //! EquilibriumState(const EquilibriumState& other): m_has_been_properly_initialised(other.m_has_been_properly_initialised), m_main_variables(other.m_main_variables), m_concentration_secondary(other.m_concentration_secondary), m_ionic_strength(other.m_ionic_strength), m_loggamma(other.m_loggamma), m_data(other.m_data) {} //! \brief Create an equilibrium state //! //! This constructor is to be used for an aqueous computation template EquilibriumState(const Eigen::MatrixBase& main_vars, const Eigen::MatrixBase& secondary, const Eigen::MatrixBase& loggamma, scalar_t ionic_strength, std::shared_ptr data): m_has_been_properly_initialised(true), m_main_variables(main_vars), m_concentration_secondary(secondary), m_ionic_strength(ionic_strength), m_loggamma(loggamma), m_data(data) {} //! \brief Create an equilibrium state //! //! This constructor is to be used for an aqueous computation EquilibriumState(const Vector& main_vars, const Vector& secondary, const Vector& loggamma, scalar_t ionic_strength, std::shared_ptr data): m_has_been_properly_initialised(true), m_main_variables(main_vars), m_concentration_secondary(secondary), m_ionic_strength(ionic_strength), m_loggamma(loggamma), m_data(data) {} //! \brief Create an equilibrium state //! //! This constructor is to be used for a computation with gas phases EquilibriumState(const Vector& main_vars, const Vector& secondary, const Vector& loggamma, scalar_t ionic_strength, const Vector& pressures, scalar_t volume, std::shared_ptr data): SimulationBox(volume), m_has_been_properly_initialised(true), m_main_variables(main_vars), m_concentration_secondary(secondary), m_pressure_gas(pressures), m_ionic_strength(ionic_strength), m_loggamma(loggamma), m_data(data) {} //! \brief Return true if it contains actual data bool is_valid() {return m_has_been_properly_initialised;} // Main variables // ============== //! \brief Return the mass of free water in the system scalar_t mass_water() const {return m_main_variables(0);} //! \brief Return the molality of component i scalar_t molality_component(index_t i) const {return POW10(m_main_variables(i));} //! \brief Return log10 of the molality of component i scalar_t log_molality_component(index_t i) const {return m_main_variables(i);} //! \brief Return the number of moles of the mineral at equilibrium 'm' scalar_t moles_mineral(index_t m) const { assert(m >= 0 and m nb_mineral); return m_main_variables(m_data->nb_component+m);} //! \brief Return the mass of a mineral (kg) scalar_t mass_mineral(index_t m) const; //! \brief Return the volume of a mineral (V) scalar_t volume_mineral(index_t m) const { return m_data->molar_volume_mineral(m)*moles_mineral(m); } //! Return the volume of all minerals scalar_t volume_minerals() const; const Vector& main_variables() const {return m_main_variables;} // molality // ======== //! \brief return the pH scalar_t pH() const; //! \brief Return the molality of secondary species j scalar_t molality_secondary(int j) const {return m_concentration_secondary(j);} //! \brief Return a const reference to the vector of molality of the secondary array const Vector& molality_secondary() const {return m_concentration_secondary;} //! \brief Return total concentration of aqueous species for component i scalar_t total_aqueous_concentration_component(int i) const; //! \brief Compute total aqeuous concentration for each component void total_aqueous_concentrations(Vector& totconc) const; // Activity // ======== scalar_t ionic_strength() const {return m_ionic_strength;} //! \brief Return a const reference to the activity coefficient vector const Vector& activity_coefficient() const {return m_loggamma;} //! \brief Return activity coefficient of component i scalar_t activity_coefficient_component(int i) const { assert(i >= 1); return POW10(m_loggamma(i)); } //! \brief Return logarithm of activity coefficient of component i scalar_t loggamma_component(int i) const { assert(i >= 1 and i < m_data->nb_component); return m_loggamma(i); } //! \brief Return activity coefficient of secondary species j scalar_t activity_coefficient_secondary(int j) const { assert(j >= 0 and j < m_data->nb_aqueous); return POW10(m_loggamma(j+m_data->nb_component)); } //! \brief Return the log of the activity coefficient of secondary species j scalar_t loggamma_secondary(int j) const { assert(j >= 0 and j < m_data->nb_aqueous); return m_loggamma(j+m_data->nb_component); } //! \brief Return activity of component i (i > 1) scalar_t activity_component(int i) const { assert(i >= 1); return POW10(m_loggamma(i) + m_main_variables(i)); } //! \brief Return activity of secondary species j scalar_t activity_secondary(int j) const { assert(j >= 0 and j < m_data->nb_aqueous); return POW10(m_loggamma(j+m_data->nb_component))*m_concentration_secondary(j); } // Minerals // ======== //! \brief return log(IAP) for mineral m scalar_t logIAP(index_t m) const; //! \brief return log(IAP) for mineral "label mineral" scalar_t logIAP(const std::string& label_mineral) const; //! \brief return saturation index for mineral m scalar_t saturation_index(index_t m) const {return logIAP(m) - m_data->logk_mineral(m);} //! \brief return saturation index for mineral "label mineral" scalar_t saturation_index(const std::string& label_mineral) const; //! \brief return log(IAP) for mineral m scalar_t logIAP_kinetic(index_t m) const; //! \brief return log(IAP) for mineral "label_mineral" scalar_t logIAP_kinetic(const std::string& label_mineral) const; //! \brief return saturation index for mineral m scalar_t saturation_index_kinetic(int m) const {return logIAP(m) - m_data->logk_mineral_kinetic(m);} //! \brief return saturation index for mineral m scalar_t saturation_index_kinetic(const std::string& label_mineral) const; // Gas phase // ========= //! \brief Return the pressure of gas phase 'g' scalar_t pressure_gas(int g) const {return m_pressure_gas(g);} //! \brief Return a const reference to the vector of gas pressure const Vector& pressure_gas() const {return m_pressure_gas;} // Thermo variables // ================ //! \brief Return the temperature scalar_t temperature() const {return get_temperature();} //! \brief Return the total volume scalar_t volume() const {return get_volume();} //! \brief Return the volume of the condensed phases scalar_t condensed_phase_volume() const; //! \brief Return the volume of the gas phase scalar_t gas_phase_volume() const; // Modification // ============ void scale_condensed(scalar_t scale); // Database access // =============== - std::shared_ptr get_database() {return m_data;} + std::shared_ptr get_database() const {return m_data;} private: bool m_has_been_properly_initialised; Vector m_main_variables; Vector m_concentration_secondary; Vector m_pressure_gas; scalar_t m_ionic_strength; Vector m_loggamma; database::RawDatabasePtr m_data; }; } // end namespace specmicp #undef POW10 #endif // SPECMICP_SPECMICP_EQUILIBRIUMDATA_HPP diff --git a/src/specmicp/kinetics/kinetic_model.hpp b/src/specmicp/kinetics/kinetic_model.hpp index ea5a7cb..ab4fba7 100644 --- a/src/specmicp/kinetics/kinetic_model.hpp +++ b/src/specmicp/kinetics/kinetic_model.hpp @@ -1,47 +1,53 @@ #ifndef SPECMICP_SPECMICP_KINETICS_KINETICMODEL_HPP #define SPECMICP_SPECMICP_KINETICS_KINETICMODEL_HPP #include "common.hpp" #include "kinetic_variables.hpp" namespace specmicp { namespace kinetics { // \brief Base class for a kinetic model class KineticModel { public: - KineticModel(index_t nb_equation): - m_neq(nb_equation) + + KineticModel() + {} + + KineticModel(const std::vector& list_species): + m_list_species(list_species) {} virtual ~KineticModel() {} //! \brief Return the number of kinetic minerals included in the problem index_t get_neq() {return m_list_species.size();} //! \brief Compute the kinetic rates and store them in dydt virtual void compute_rate( scalar_t t, const Vector& y, KineticVariables& variables, Vector& dydt ) = 0; + void add_equation(index_t species) {m_list_species.push_back(species);} + //! \brief Index of species corresponding to equation 'id_equation' index_t index_species(index_t id_equation) {return m_list_species[id_equation];} //! \brief Iterator over the species vector std::vector::iterator species_begin() {return m_list_species.begin();} std::vector::iterator species_end() {return m_list_species.end();} private: - std::vector list_species; + std::vector m_list_species; }; } // end namespace kinetics } // end namespace specmicp #endif //SPECMICP_SPECMICP_KINETICS_KINETICMODEL_HPP diff --git a/src/specmicp/kinetics/kinetic_system.cpp b/src/specmicp/kinetics/kinetic_system.cpp index b04c5c5..dda24f7 100644 --- a/src/specmicp/kinetics/kinetic_system.cpp +++ b/src/specmicp/kinetics/kinetic_system.cpp @@ -1,59 +1,72 @@ #include "kinetic_system.hpp" #include "specmicp/reduced_system_solver.hpp" #include "kinetic_model.hpp" namespace specmicp { namespace kinetics { -void KineticSystem::update_total_concentrations(scalar_t dt, const Vector& y) +void KineticSystem::update_total_concentrations(const Vector& y) { - Vector update_minerals = (y - m_variables.moles_minerals())/dt; - RawDatabasePtr database = m_variables.get_database(); + //Vector y_temp = y.cwiseMax(Vector::Zero(y.rows())); + Vector update_minerals = (y - m_variables.moles_minerals()); m_variables.total_concentrations() = m_variables.total_concentrations_initial(); for (auto it=m_model->species_begin(); it!=m_model->species_end(); ++it) { - for (index_t component: database->range_component()) + for (index_t component: m_data->range_component()) { - m_variables.total_concentration(component) += - database->nu_mineral_kinetic(*it, component)*update_minerals(it); + if (m_data->nu_mineral_kinetic(*it, component) == 0.0) continue; + m_variables.total_concentration(component) -= + m_data->nu_mineral_kinetic(*it, component)*update_minerals(it-m_model->species_begin()); } } } +void KineticSystem::update_to_new_initial_condition(const Vector& y, scalar_t dt) +{ + update_total_concentrations(y); + m_variables.rate_components() += ( + m_variables.total_concentrations() + - m_variables.total_concentrations_initial() + )/dt; + m_variables.moles_minerals() = y; + m_variables.total_concentrations_initial() = m_variables.total_concentrations(); +} + void KineticSystem::compute_rates(scalar_t t, const Vector& y, Vector& dydt) { m_model->compute_rate(t, y, m_variables, dydt); } void KineticSystem::compute_equilibrium(ReducedSystemSolverOptions options) { ReducedSystemSolver solver; Vector variables; specmicp::micpsolver::MiCPPerformance perf; if (m_variables.equilibrium_state().is_valid()) { - solver = ReducedSystemSolver(m_variables.get_database(), + solver = ReducedSystemSolver(m_data, m_variables.total_concentrations(), m_variables.equilibrium_state(), options); variables = m_variables.equilibrium_state().main_variables(); perf = solver.solve(variables); } else { - solver = ReducedSystemSolver(m_variables.get_database(), + solver = ReducedSystemSolver(m_data, m_variables.total_concentrations(), options); perf = solver.solve(variables, true); } + //std::cout << variables << std::endl; if ((int) perf.return_code < 0) { std::cout << "Failed to solve the system ! Err code " << (int) perf.return_code << std::endl; } m_variables.update_equilibrium_state(solver.get_solution(variables)); } } // end namespace kinetics } // end namespace specmicp diff --git a/src/specmicp/kinetics/kinetic_system.hpp b/src/specmicp/kinetics/kinetic_system.hpp index e18c60d..a9862ba 100644 --- a/src/specmicp/kinetics/kinetic_system.hpp +++ b/src/specmicp/kinetics/kinetic_system.hpp @@ -1,48 +1,79 @@ #ifndef SPECMICP_SPECMICP_KINETICS_KINETICSYSTEM_HPP #define SPECMICP_SPECMICP_KINETICS_KINETICSYSTEM_HPP #include #include "common.hpp" #include "database.hpp" #include "kinetic_variables.hpp" #include "specmicp/reduced_system_solver_structs.hpp" namespace specmicp { namespace kinetics { class KineticModel; //! \brief Kinetics System, //! //! Wrapper for a kinetic model class KineticSystem { public: - KineticSystem(std::shared_ptr model - ): - m_model(model) + KineticSystem( + std::shared_ptr model, + const Vector& total_concentrations, + const Vector& mineral_moles, + RawDatabasePtr data + ): + m_data(data), + m_model(model), + m_variables(total_concentrations, mineral_moles) {} + + KineticSystem( + std::shared_ptr model, + const Vector& total_concentrations, + const Vector& mineral_moles, + const EquilibriumState& equilibrium_solution + ): + m_data(equilibrium_solution.get_database()), + m_model(model), + m_variables(total_concentrations, mineral_moles, equilibrium_solution) + {} + //! \brief Compute the kinetics rates to be solved //! //! Use the kinetic model provided by the user - void compute_rates(scalar_t x, const Vector& y, Vector& dydx); + void compute_rates(scalar_t x, const Vector& y, Vector& dydt); //! \brief Compute the equilibrium state of the solution void compute_equilibrium(ReducedSystemSolverOptions options=ReducedSystemSolverOptions()); - //! Update the total concentrations - void update_total_concentrations(scalar_t dt, const Vector& y); + //! \brief Update the total concentrations + //! + //! \param y vector of variables (mols of minerals) + void update_total_concentrations(const Vector& y); + + void update_to_new_initial_condition(const Vector& y, scalar_t dt); + + //! \brief Right Hand side function for the integrator + void rhs(scalar_t x, const Vector& y, Vector& dydt, ReducedSystemSolverOptions options=ReducedSystemSolverOptions()) + { + update_total_concentrations(y); + compute_equilibrium(options); + compute_rates(x, y, dydt); + } KineticVariables& variables() {return m_variables;} private: - KineticVariables m_variables; + RawDatabasePtr m_data; std::shared_ptr m_model; + KineticVariables m_variables; }; } // end namespace kinetics } // end namespace specmicp #endif //SPECMICP_SPECMICP_KINETICS_KINETICSYSTEM_HPP diff --git a/src/specmicp/kinetics/kinetic_system_solver.cpp b/src/specmicp/kinetics/kinetic_system_solver.cpp index e69de29..0199f83 100644 --- a/src/specmicp/kinetics/kinetic_system_solver.cpp +++ b/src/specmicp/kinetics/kinetic_system_solver.cpp @@ -0,0 +1,49 @@ +#include "kinetic_system_solver.hpp" + +#include "kinetic_model.hpp" +#include "odeint/runge_kutta_step.hpp" + +#include + +namespace specmicp { +namespace kinetics { + + +void KineticSystemSolver::solve(scalar_t dt, scalar_t total) +{ + double t = 0.0; + Vector y = m_system.variables().moles_minerals(); + Vector dydx(y.rows()); + dydx.setZero(); + m_system.rhs(t, y, dydx, get_options().speciation_options); + + + odeint::rhs_f rhs = std::bind(std::mem_fn(&KineticSystem::rhs), + m_system, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3, + get_options().speciation_options); + odeint::CashKarpStep driver = odeint::get_cash_karp_step(rhs); + + odeint::StepLength stepl(dt, total); + int cnt = 0; + while (t < total) + { + driver.rk_step(y, dydx, t, stepl); + m_system.update_to_new_initial_condition(y, stepl.did); + cnt +=1; + scalar_t remaining = stepl.remaining(t); + if (remaining == 0.0) break; + if (stepl.next > remaining) stepl.next = remaining; + stepl.get_next(); + m_current_dt = stepl.test; + + } + m_system.update_total_concentrations(y); + m_system.compute_equilibrium(get_options().speciation_options); + +} + +} // end namespace kinetics +} // end namespace specmicp diff --git a/src/specmicp/kinetics/kinetic_system_solver.hpp b/src/specmicp/kinetics/kinetic_system_solver.hpp index beb6af9..3476755 100644 --- a/src/specmicp/kinetics/kinetic_system_solver.hpp +++ b/src/specmicp/kinetics/kinetic_system_solver.hpp @@ -1,4 +1,51 @@ #ifndef SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVER_HPP #define SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVER_HPP +#include "kinetic_system.hpp" +#include "kinetic_system_solver_structs.hpp" + +#include "utils/options_handler.hpp" + +namespace specmicp { +namespace kinetics { + +class KineticModel; + +class KineticSystemSolver: public OptionsHandler +{ +public: + KineticSystemSolver( + std::shared_ptr model, + const Vector& total_concentrations, + const Vector& mineral_moles, + RawDatabasePtr data + ): + m_current_dt(-1), + m_system(model, total_concentrations, mineral_moles, data) + {} + + KineticSystemSolver( + std::shared_ptr model, + const Vector& total_concentrations, + const Vector& mineral_moles, + const EquilibriumState& equilibrium_solution + ): + m_current_dt(-1), + m_system(model, total_concentrations, mineral_moles, equilibrium_solution) + {} + + void solve(scalar_t dt, scalar_t total); + + scalar_t current_dt() {return m_current_dt;} + + KineticVariables& variables() {return m_system.variables();} + +private: + scalar_t m_current_dt; + KineticSystem m_system; +}; + +} // end namespace kinetics +} // end namespace specmicp + #endif //SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVER_HPP diff --git a/src/specmicp/kinetics/kinetic_variables.hpp b/src/specmicp/kinetics/kinetic_variables.hpp index fc4cf78..f036ce0 100644 --- a/src/specmicp/kinetics/kinetic_variables.hpp +++ b/src/specmicp/kinetics/kinetic_variables.hpp @@ -1,75 +1,108 @@ #ifndef SPECMICP_SPECMICP_KINETICS_KINETICVARIABLES_HPP #define SPECMICP_SPECMICP_KINETICS_KINETICVARIABLES_HPP #include "common.hpp" #include "specmicp/equilibrium_data.hpp" namespace specmicp { namespace kinetics { //! \brief Variables used in a kinetic computation class KineticVariables { public: KineticVariables( const Vector& total_concentrations, - const Vector& mineral_moles): + const Vector& mineral_moles + ): m_mineral_kinetics(mineral_moles), m_total_concentrations_initial(total_concentrations), - m_total_concentrations(total_concentrations) + m_total_concentrations(total_concentrations), + m_rate(Vector::Zero(total_concentrations.rows())) {} KineticVariables( const Vector& total_concentrations, const Vector& mineral_moles, const EquilibriumState& equilibrium_solution ): - KineticVariables(total_concentrations, mineral_moles), + m_mineral_kinetics(mineral_moles), + m_total_concentrations_initial(total_concentrations), + m_total_concentrations(total_concentrations), + m_rate(Vector::Zero(total_concentrations.rows())), m_equilibrium(equilibrium_solution) {} //! \brief Return the number of species considered in the kinetics computation index_t get_neq() {return m_mineral_kinetics.rows();} //! \brief Return the raw database RawDatabasePtr get_database() {return m_equilibrium.get_database();} + // Number of moles of minerals + // --------------------------- //! \brief Return mole number of mineral "mineral kinetic" scalar_t moles_mineral(index_t mineral_kinetic) const {return m_mineral_kinetics(mineral_kinetic);} //! \brief Return a const reference to the vector of moles number const Vector& moles_minerals() const {return m_mineral_kinetics;} //! \brief Return a reference to the vector of moles number Vector& moles_minerals() {return m_mineral_kinetics;} + // Total concentrations + // -------------------- + //! \brief Return the total concentration (at equilibrium) of 'component' scalar_t total_concentration(index_t component) const {return m_total_concentrations(component);} + //! \brief Return the total concentration (at equilibrium) of 'component' + scalar_t& total_concentration(index_t component) {return m_total_concentrations(component);} //! \brief Return a const reference to the total concentration (at equilibrium) vector const Vector& total_concentrations() const {return m_total_concentrations;} //! \brief Return a reference to the total concentration at equilibrium vector Vector& total_concentrations() {return m_total_concentrations; } + // Initial total concentrations + // ---------------------------- + //! \brief Return the total concentration (at equilibrium) of 'component' scalar_t total_concentration_initial(index_t component) const { return m_total_concentrations_initial(component);} //! \brief Return a const reference to the total concentration (at equilibrium) vector const Vector& total_concentrations_initial() const {return m_total_concentrations_initial;} //! \brief Return a reference to the total concentration at equilibrium vector Vector& total_concentrations_initial() {return m_total_concentrations_initial; } + // Rates + // ----- + + //! \brief Return the value of the flux for 'component' + scalar_t rate_component(index_t component) const { + return m_rate(component); + } + //! \brief Return a const reference to the flux vector + const Vector& rate_components() const {return m_rate;} + //! \brief Return a reference to the flux vector + Vector& rate_components() {return m_rate;} + //! \brief reset the rates + void reset_rate() {m_rate.setZero();} + + // Equilibrium + // ------------ + //! \brief Const Reference to the equilibrium state const EquilibriumState& equilibrium_state() const {return m_equilibrium;} //! \brief Reference to the equilibrium State EquilibriumState& equilibrium_state() {return m_equilibrium;} //! \brief Update the equilibrium state void update_equilibrium_state(const EquilibriumState& other) {m_equilibrium = other;} private: Vector m_mineral_kinetics; //!< Number of moles of minerals at equilibrium Vector m_total_concentrations_initial; //!< Initial total concentrations at equilibrium - Vector m_total_concentrations; //!< total concentration at equilibrium + Vector m_total_concentrations; //!< total concentration at equilibrium, + Vector m_rate; EquilibriumState m_equilibrium; //!< Current equilibrium state }; } // end namespace kinetics } // end namespace specmicp #endif // SPECMICP_SPECMICP_KINETICS_KINETICVARIABLES_HPP diff --git a/tests/specmicp/kinetics.cpp b/tests/specmicp/kinetics.cpp index 2e31c17..69990c7 100644 --- a/tests/specmicp/kinetics.cpp +++ b/tests/specmicp/kinetics.cpp @@ -1,221 +1,358 @@ /*------------------------------------------------------- - 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 "database.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") { 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, 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 << "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 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; + 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; + solver.get_options().charge_keeper = 1; 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() +void by_hand() { - - 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; +} + + +/////// using the new solver +/// +#include "specmicp/kinetics/kinetic_model.hpp" +#include "specmicp/kinetics/kinetic_variables.hpp" +#include "specmicp/kinetics/kinetic_system_solver.hpp" + +class C3SDissolution: public specmicp::kinetics::KineticModel +{ +public: + static constexpr scalar_t R0 = 2.5e-6; + + C3SDissolution(specmicp::RawDatabasePtr data, scalar_t n_c3s0) + { + specmicp::Database db(data); + add_equation(db.mineral_kinetic_label_to_id("C3S")); + m_n_c3s0 = n_c3s0; + n_particle = 3.0/(4.0*M_PI*std::pow(R0, 3.0))*Vmc3s*n_c3s0; + + } + + void compute_rate( + scalar_t _, + const specmicp::Vector& y, + specmicp::kinetics::KineticVariables& variables, + specmicp::Vector& dydt + ) + { + specmicp::Vector y_temp = y; + //specmicp::odeint::set_non_negative(y_temp); + //std::cout << "y_temp : " << y_temp << std::endl; + dydt.resizeLike(y); + double radius = R0*std::pow((y_temp(0)/m_n_c3s0), 1.0/3.0); + if (radius < 0.0 or not std::isfinite(radius)) radius = 0.0; + double surf = n_particle*4*M_PI*radius*radius; + //std::cout << n_particle << "\t" << radius << "\t" << surf << std::endl; + dydt(0) = 1e-9*std::min(surf*(2.1*variables.equilibrium_state().logIAP_kinetic(index_species(0)) +35.7), -1e-6); + } + +private: + double m_n_c3s0; + double n_particle; +}; + +void using_kinetic_solver() +{ + specmicp::Database the_database("../data/cemdata_specmicp.js"); + std::shared_ptr data = the_database.get_database(); + std::map swapping ({ + {"H[+]","HO[-]"}, + {"Si(OH)4", "SiO(OH)3[-]"} + }); + the_database.swap_components(swapping); + std::vector to_keep = {0, 1, + the_database.safe_label_to_id("Ca[2+]", the_database.get_database()->labels_basis), + the_database.safe_label_to_id("SiO(OH)3[-]", the_database.get_database()->labels_basis) + }; + the_database.keep_only_components(to_keep); + + + specmicp::Vector totconc(data->nb_component); + totconc(the_database.component_label_to_id("H2O")) = 1.0/data->molar_mass_basis_si(0); + totconc(the_database.component_label_to_id("HO[-]")) = 1e-7; + totconc(the_database.component_label_to_id("Ca[2+]")) = 1e-4; + totconc(the_database.component_label_to_id("SiO(OH)3[-]")) =2e-4; + + specmicp::ReducedSystemSolver first_solver(data, totconc); + + specmicp::Vector minerals(1); + minerals << 1e-3*(500/Mc3s); + std::cout << "Start : " << minerals(0) << std::endl; + + specmicp::Vector equilvar; + equilvar.resize(data->nb_component+data->nb_mineral); + equilvar(0) = 1.0; + equilvar.block(1, 0, data->nb_component-1, 1).setConstant(-7); + equilvar.block(data->nb_component, 0, data->nb_mineral, 1).setConstant(0.); + + first_solver.solve(equilvar); + + specmicp::Vector second_conc(data->nb_aqueous); + second_conc.setZero(); + + specmicp::Vector loggamma(data->nb_component+data->nb_aqueous); + loggamma.setZero(); + + specmicp::EquilibriumState equilsolution = first_solver.get_solution(equilvar); + + std::cout << equilvar << std::endl; + std::cout << "Initialization finished " << std::endl; + + std::shared_ptr model = std::make_shared(data, minerals(0)); + specmicp::kinetics::KineticSystemSolver solver(model, totconc, minerals, equilsolution); + solver.get_options().speciation_options.solver_options.maxstep = 50; + solver.get_options().speciation_options.solver_options.use_scaling = true; + solver.get_options().speciation_options.solver_options.fvectol = 1e-8; + solver.get_options().speciation_options.charge_keeper = the_database.component_label_to_id("HO[-]"); + solver.get_options().speciation_options.solver_options.non_monotone_linesearch = true; + solver.get_options().speciation_options.solver_options.projection_min_variable = 1e-8; + //solver.get_options().speciation_options.solver_options.threshold_cycling_linesearch = 1.0; + //solver.get_options().speciation_options.conservation_water = false; + solver.get_options().rk_options.non_negativity = true; + solver.solve(1000.0, 100000.0); + // solver.solve(1.0, 2000.0); + // solver.solve(1.0, 2000.0); + std::cout << "Solution : " << solver.variables().moles_mineral(0) << std::endl; + std::cout << "flux : " << std::endl << solver.variables().rate_components() << std::endl; + solver.variables().reset_rate(); + solver.solve(solver.current_dt(), 100000.0); + // solver.solve(1.0, 2000.0); + // solver.solve(1.0, 2000.0); + std::cout << "Solution : " << solver.variables().moles_mineral(0) << std::endl; + std::cout << "flux : " << std::endl << solver.variables().rate_components() << std::endl; + solver.variables().reset_rate(); + solver.solve(solver.current_dt(), 100000.0); + // solver.solve(1.0, 2000.0); + // solver.solve(1.0, 2000.0); + std::cout << "Solution : " << solver.variables().moles_mineral(0) << std::endl; + std::cout << "flux : " << std::endl << solver.variables().rate_components() << std::endl; + solver.variables().reset_rate(); + solver.solve(solver.current_dt(), 100000.0); + // solver.solve(1.0, 2000.0); + // solver.solve(1.0, 2000.0); + std::cout << "Solution : " << solver.variables().moles_mineral(0) << std::endl; + std::cout << "flux : " << std::endl << solver.variables().rate_components() << std::endl; +} + +int main() +{ + + specmicp::logger::ErrFile::stream() = &std::cerr; + specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; + + using_kinetic_solver(); + //by_hand(); + return EXIT_SUCCESS; }