diff --git a/src/specmicp/reaction_path.cpp b/src/specmicp/reaction_path.cpp index 00c8a03..5a70c99 100644 --- a/src/specmicp/reaction_path.cpp +++ b/src/specmicp/reaction_path.cpp @@ -1,82 +1,82 @@ /*------------------------------------------------------- - Module : specmicp - File : reaction_path.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "reaction_path.hpp" #include "thermodata.hpp" #include "micpsolver/micpsolver.hpp" #include "utils/log.hpp" #include <iostream> namespace specmicp { void ReactionPathDriver::read_database() { m_database.parse_database(m_model->database_path); m_database.make_canonical(); m_data = m_database.get_database(); } int ReactionPathDriver::label_to_id(const std::string& label, const std::vector<std::string>& list_labels) { auto it = std::find(list_labels.begin(), list_labels.end(), label); if (it == list_labels.end()) { throw std::invalid_argument("Unknow species : " + label); } return it - list_labels.begin(); } void ReactionPathDriver::dissolve_to_components() { m_tot_conc = Eigen::VectorXd(m_data->nb_component); m_tot_conc_increase = Eigen::VectorXd(m_data->nb_component); // everything is dissolved ... for (auto it=m_model->amount_components.begin(); it!=m_model->amount_components.end(); ++it) { const int id = label_to_id(it->first, m_data->labels_basis); m_tot_conc(id) = it->second.first; m_tot_conc_increase(id) = it->second.second; } std::vector<int> to_remove; to_remove.reserve(10); int new_i = 0; for (int i=0; i<m_data->nb_component; ++i) { if (m_tot_conc(i) == 0 and m_tot_conc_increase(i) == 0) { to_remove.push_back(0); continue; } m_tot_conc(new_i) = m_tot_conc(i); m_tot_conc_increase(new_i) = m_tot_conc_increase(i); ++new_i; } m_tot_conc.conservativeResize(new_i); m_tot_conc_increase.conservativeResize(new_i); m_database.remove_components(to_remove); } micpsolver::MiCPPerformance ReactionPathDriver::one_step(Eigen::VectorXd& x) { m_current_solver = ReducedSystemSolver(m_data, m_tot_conc); - m_current_solver.get_options().solver_options = m_solver_options; + m_current_solver.get_options() = m_solver_options; micpsolver::MiCPPerformance perf = m_current_solver.solve(x); if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) { ERROR << "Failed to solve the system !"; } m_tot_conc += m_tot_conc_increase; return perf; } } // end namespace specmicp diff --git a/src/specmicp/reaction_path.hpp b/src/specmicp/reaction_path.hpp index c78df8c..1022f77 100644 --- a/src/specmicp/reaction_path.hpp +++ b/src/specmicp/reaction_path.hpp @@ -1,106 +1,112 @@ /*------------------------------------------------------- - Module : specmicp - File : reaction_path.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_SPECMICP_REACTIONPATH_HPP #define SPECMICP_SPECMICP_REACTIONPATH_HPP //! \file reaction_path.hpp reaction path driver #include <memory> #include <map> #include "database/database.hpp" #include "micpsolver/micpsolver_structs.hpp" #include "reduced_system_solver.hpp" namespace specmicp { //! \brief typedef to a reaction path amount //! //! The first element is the inital amount, //! The second element is the increase to add at every step using reaction_amount_t = std::pair<double, double>; //! \brief A reaction path model //! //! Contains every data we need to run a reaction path model struct ReactionPathModel { //! The amount of components std::map<std::string, reaction_amount_t> amount_components; //! The amount of aqueous species std::map<std::string, reaction_amount_t> amount_aqueous; //! The amount of minerals std::map<std::string, reaction_amount_t> amount_minerals; //! the number of step int nb_step; //! the path to the database - optional std::string database_path; }; //! \brief Driver to solve a reaction path model class ReactionPathDriver { public: //! \brief Initialise from a model //! //! Build a database from the attribute database_path of the model - ReactionPathDriver(std::shared_ptr<ReactionPathModel> model): - m_model(model) + ReactionPathDriver(std::shared_ptr<ReactionPathModel> model, + ReducedSystemSolverOptions options=ReducedSystemSolverOptions()): + m_model(model), + m_solver_options(options) { read_database(); } //! \brief Initialise the driver from a model and a raw database ReactionPathDriver(std::shared_ptr<ReactionPathModel> model, - std::shared_ptr<database::DataContainer> data): + std::shared_ptr<database::DataContainer> data, + ReducedSystemSolverOptions options=ReducedSystemSolverOptions()): m_model(model), m_database(data), - m_data(data) + m_data(data), + m_solver_options(options) {} //! Initialize the database void read_database(); //! Dissolve everything into components void dissolve_to_components(); //! \brief Return the id of a species in 'list_labels' from its id int label_to_id(const std::string& label, const std::vector<std::string>& list_labels); //! \brief Perform one step micpsolver::MiCPPerformance one_step(Eigen::VectorXd& x); //! Set total aqueous concentrations void total_aqueous_concentration(const Eigen::VectorXd& x, Eigen::VectorXd& totaq) { m_current_solver.total_aqueous_concentration(x, totaq); } + const ReducedSystemSolverOptions& get_options() const {return m_solver_options;} + ReducedSystemSolverOptions& get_options() {return m_solver_options;} private: std::shared_ptr<ReactionPathModel> m_model; specmicp::database::Database m_database; std::shared_ptr<specmicp::database::DataContainer> m_data; Eigen::VectorXd m_tot_conc; Eigen::VectorXd m_tot_conc_increase; - micpsolver::MiCPSolverOptions m_solver_options; + ReducedSystemSolverOptions m_solver_options; ReducedSystemSolver m_current_solver; //std::shared_ptr<ThermoData> m_thermo; }; } // end namespace specmicp #endif // SPECMICP_SPECMICP_REACTIONPATH_HPP