diff --git a/CMakeLists.txt b/CMakeLists.txt index 62af00b..f03c32a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,171 +1,173 @@ 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} -O2 -pg -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}/extended_system.cpp ${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}/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}/data_container.hpp + ${DATABASE_DIR}/module.hpp +) add_library(specmicp_database ${DATABASE_LIBFILE}) # ------- 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) 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) diff --git a/src/database/canonicalizer.cpp b/src/database/canonicalizer.cpp index 3550549..0aa2dcc 100644 --- a/src/database/canonicalizer.cpp +++ b/src/database/canonicalizer.cpp @@ -1,104 +1,116 @@ /*------------------------------------------------------- - Module : database - File : canonicalizer.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "canonicalizer.hpp" #include "utils/log.hpp" namespace specmicp { namespace database { +void DatabaseCanonicalizer::make_canonical() { + if (data->is_canonical == true) + { + WARNING << "Database already in canonical form, aborting"; + return; + } + make_aqueous_canonical(); + make_mineral_canonical(); + make_mineral_kinetic_canonical(); + data->is_canonical = true; +} + void DatabaseCanonicalizer::make_aqueous_canonical() { if (data->nu_aqueous.cols() == data->nb_component) { INFO << "Database already in canonical form"; return; } if (data->nu_aqueous.cols() != data->nb_component + data->nb_aqueous) { CRITICAL << "Canonicalization : invalid size for stoichiometric coefficients of aqueous species - abort"; throw invalid_database("Invalid size for stoichiometric coefficients of aqueous species"); } for (int i=1; inb_aqueous; ++i) { for (int j=0; jnu_aqueous(i, data->nb_component+j); if (stoech != 0) { data->nu_aqueous.row(i) += stoech*data->nu_aqueous.row(j); data->logk_aqueous(i) += stoech*data->logk_aqueous(j); data->nu_aqueous(i, data->nb_component+j) = 0; } } } data->nu_aqueous.conservativeResize(Eigen::NoChange_t(), data->nb_component); } void DatabaseCanonicalizer::make_mineral_canonical() { if (data->nu_mineral.cols() == data->nb_component) { INFO << "Database already in canonical form"; return; } if (data->nu_mineral.cols() != data->nb_component + data->nb_aqueous) { CRITICAL << "Canonicalization : invalid size for stoichiometric coefficients of mineral species - abort"; throw invalid_database("Invalid size for stoichiometric coefficients of mineral species"); } for (int i=1; inb_mineral; ++i) { for (int j=0; jnb_aqueous; ++j) { const double stoech = data->nu_mineral(i, data->nb_component+j); if (stoech != 0) { data->nu_mineral.block(i, 0, 1, data->nb_component) += stoech*data->nu_aqueous.row(j); data->logk_mineral(i) += stoech*data->logk_aqueous(j); data->nu_mineral(i, data->nb_component+j) = 0; } } } data->nu_mineral.conservativeResize(Eigen::NoChange_t(), data->nb_component); } void DatabaseCanonicalizer::make_mineral_kinetic_canonical() { if (data->nu_mineral_kinetic.cols() == data->nb_component) { INFO << "Database already in canonical form"; return; } if (data->nu_mineral_kinetic.cols() != data->nb_component + data->nb_aqueous) { CRITICAL << "Canonicalization : invalid size for stoichiometric coefficients of mineral species - abort"; throw invalid_database("Invalid size for stoichiometric coefficients of mineral species"); } for (int i=1; inb_mineral_kinetic; ++i) { for (int j=0; jnb_aqueous; ++j) { const double stoech = data->nu_mineral_kinetic(i, data->nb_component+j); if (stoech != 0) { data->nu_mineral_kinetic.block(i, 0, 1, data->nb_component) += stoech*data->nu_aqueous.row(j); data->logk_mineral_kinetic(i) += stoech*data->logk_aqueous(j); data->nu_mineral_kinetic(i, data->nb_component+j) = 0; } } } data->nu_mineral_kinetic.conservativeResize(Eigen::NoChange_t(), data->nb_component); } } // end namespace database } // end namespace specmicp diff --git a/src/database/canonicalizer.hpp b/src/database/canonicalizer.hpp index 4bd5402..306b7b6 100644 --- a/src/database/canonicalizer.hpp +++ b/src/database/canonicalizer.hpp @@ -1,64 +1,58 @@ /*------------------------------------------------------- - Module : database - File : canonicalizer - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_DATABASE_CANONICALIZER #define SPECMICP_DATABASE_CANONICALIZER //! \file canonicalizer Transform a database into canonical form -#include "memory" -#include "data_container.hpp" +#include "module.hpp" namespace specmicp { namespace database { //! \brief Transform a database into canonical form //! //! Canonical form means that every secondary species are written as function of the components //! //! P.S : what a great name -class DatabaseCanonicalizer +class DatabaseCanonicalizer: public DatabaseModule { public: - DatabaseCanonicalizer(std::shared_ptr raw_data): - data(raw_data) + DatabaseCanonicalizer(RawDatabasePtr thedata): + DatabaseModule(thedata) {} //! \brief Transform the database into canonical form //! //! This process is necessary before any computation - void make_canonical() { - make_aqueous_canonical(); - make_mineral_canonical(); - make_mineral_kinetic_canonical(); - } + void make_canonical(); private: //! \brief transform the aqueous system into a canonical form void make_aqueous_canonical(); //! \brief transform the mineral system into a canonical form void make_mineral_canonical(); //! \brief transform the mineral (governed by kineticS) system into a canonical form void make_mineral_kinetic_canonical(); - std::shared_ptr data; //!< The raw database }; -inline void canonicalize_database(std::shared_ptr data) +inline void canonicalize_database(DatabaseCanonicalizer::RawDatabasePtr data) { DatabaseCanonicalizer(data).make_canonical(); } } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_CANONICALIZER diff --git a/src/database/data_container.hpp b/src/database/data_container.hpp index aa68995..2a9b99a 100644 --- a/src/database/data_container.hpp +++ b/src/database/data_container.hpp @@ -1,56 +1,62 @@ /*------------------------------------------------------- - Module : database - File : data_container.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_DATABASE_DATACONTAINER_HPP #define SPECMICP_DATABASE_DATACONTAINER_HPP //! \file data_container.hpp Storage class for thermodynamics database #include "common_def.hpp" #include namespace specmicp { namespace database { //! \brief Storage class - Contains the database //! //! - Should not be accessed directly but through interfaces //! - Should be shared with a smart pointer struct DataContainer { + DataContainer(): + is_canonical(false) + {} + int nb_component; //!< Number of components == size of the basis int nb_aqueous; //!< Number of aqueous species (not taking into acount the basis) int nb_mineral; //!< Number of minerals (used in computation) int nb_mineral_kinetic; //!< Number of minerals (governed by kinetics) vector_labels_t labels_basis; //!< labels of the components std::map map_labels_basis; //!< map labels <-> id vector_labels_t labels_aqueous; //!< labels of the aqueous species vector_labels_t labels_minerals; //!< labels of the minerals vector_labels_t labels_minerals_kinetic; //!< labels of the minerals (governed by kinetics) reaction_mat_t nu_aqueous; //!< Stoechiometric coefficient for aqueous species logK_vector_t logk_aqueous; //!< LogK for aqueous species reaction_mat_t nu_mineral; //!< Stoichiometric coefficient for minerals logK_vector_t logk_mineral; //!< LogK for minerals reaction_mat_t nu_mineral_kinetic; //!< Stoichiometric coefficient for minerals (governed by kinetics) logK_vector_t logk_mineral_kinetic; //!< LogK for minerals (governed by kinetics) Eigen::Matrix param_aq; //!< Aqueous parameters + bool is_canonical; + }; } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_DATACONTAINER_HPP diff --git a/src/database/module.hpp b/src/database/module.hpp new file mode 100644 index 0000000..8b09a03 --- /dev/null +++ b/src/database/module.hpp @@ -0,0 +1,47 @@ +/*------------------------------------------------------- + + - Module : database + - File : module + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#ifndef SPECMICP_DATAABASE_MODULE_HPP +#define SPECMICP_DATAABASE_MODULE_HPP + +//! \file module base class for a database module + + +#include +#include "data_container.hpp" + +namespace specmicp { +namespace database { + + +//! \brief Base class for a database module +class DatabaseModule +{ +public: + //! the type of a smart pointer to a raw database + using RawDatabasePtr = std::shared_ptr; + + DatabaseModule() {} + DatabaseModule(RawDatabasePtr thedata): data(thedata) {} + + //! (Re)set the database + void reset_database() {data = std::make_shared();} + + //! Return true if the database is canonical + bool is_database_canonical() {return data->is_canonical;} + +protected: + RawDatabasePtr data; //! The database +}; + +} // end namespace database +} // end namespace specmicp + +#endif // SPECMICP_DATAABASE_MODULE_HPP diff --git a/src/database/reader.hpp b/src/database/reader.hpp index 52017ab..b406f52 100644 --- a/src/database/reader.hpp +++ b/src/database/reader.hpp @@ -1,130 +1,125 @@ /*------------------------------------------------------- - Module : database - File : reader.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef DATABASE_READER_H #define DATABASE_READER_H //! \file reader.hpp Read the json database #include -#include -#include "data_container.hpp" +#include "module.hpp" namespace specmicp { namespace database { //! \brief Read a json file containing a database //! //! This is a reader class, it does not transform //! (reduce to canonical form, swap the basis) the data -class DataReader +class DataReader: public DatabaseModule { public: - //! Type of a smart pointer containing the database - using DataContainerPtr = std::shared_ptr; //! \brief Default constructor //! //! Be sure to set the path to the database using set_database_path DataReader() { - data = std::make_shared(); + reset_database(); } //! \brief Constructor //! //! @param filepath string containing the path to the database DataReader(std::string filepath): m_filepath(filepath) { - data = std::make_shared(); + reset_database(); } //! \brief Set the database path void set_database_path(std::string filepath){ m_filepath = filepath; } //! \brief Parse database //! //! Throw db_invalid_syntax if a syntax error was detected. //! Throws std::invalid_argument if the path to the database is not correct void parse(); //! Return the databes - DataContainerPtr get_database() {return data;} + RawDatabasePtr get_database() {return data;} private: //! \brief Parse the metadata section //! //! we don't do much with them for now.... void parse_metadata(); //! \brief Parse the basis section //! //! Contains the list of primary species void parse_basis(); //! \brief Parse the aqueous section //! //! Contains the list of secondary species void parse_aqueous(); //! \brief Parse the mineral section //! //! Contains the list of minerals void parse_minerals(); //! \brief To be sure that the database is og void check_database(); //! \brief Check the size of the basis void parse_size_db(); std::string m_filepath; //! The path to the database Json::Value m_root; //! Root of the json structure - - DataContainerPtr data; //! The database to fill }; //! \brief Parse an equation void parse_equation(const std::string &equation, std::map& compo); //! \brief Get the charge of a species by parsing the label //! //! Examples : //! - neutral -> 0 //! - neutral[] -> 0 //! - charge[+] -> +1 //! - charge[-1] -> -1 //! - charge[+2] -> +2 //! - charge[2-] -> -2 double charge_from_label(const std::string& label); //! \brief Check if a mandatory section/attribute exists //! //! Throw db_invalid_syntax if not. inline void check_for_mandatory_value(const Json::Value& root, const std::string& attribute) { if (not root.isMember(attribute)) { throw db_invalid_syntax("Missing : mandatory parameter '" + attribute +"'."); } } //! \brief Return a mandatory section //! //! Throw db_invalid_syntax if the section does not exist inline Json::Value& get_mandatory_section(Json::Value& root, const std::string& section) { check_for_mandatory_value(root, section); return root[section]; } } // end namespace database } // end namespace specmicp #endif // DATABASE_READER_H diff --git a/src/database/selector.hpp b/src/database/selector.hpp index 352ab7d..babe1e4 100644 --- a/src/database/selector.hpp +++ b/src/database/selector.hpp @@ -1,71 +1,68 @@ /*------------------------------------------------------- - Module : database - File : selector.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_DATABASE_SELECTOR_HPP #define SPECMICP_DATABASE_SELECTOR_HPP //! \file selector.hpp select species int the database #include -#include -#include "data_container.hpp" +#include "module.hpp" namespace specmicp { namespace database { //! \brief Select species in the database //! //! This is a special class with the only function is to remove some components from the database -class DatabaseSelector +class DatabaseSelector: public DatabaseModule { public: - DatabaseSelector(std::shared_ptr thedata): - data(thedata), + DatabaseSelector(RawDatabasePtr thedata): + DatabaseModule(thedata), m_is_component_to_remove(thedata->nb_component, false) {} void remove_component(const std::vector& id_component_to_remove) { analyse_component(id_component_to_remove); select_aqueous(id_component_to_remove); select_minerals(id_component_to_remove); select_components(); } private: //! \brief Analyse wich component should be removed void analyse_component(const std::vector& id_component_to_remove); //! \brief Select the correct aqueous species void select_aqueous(const std::vector& id_component_to_remove); //! \brief Select the correct minerals void select_minerals(const std::vector& id_component_to_remove); //! \brief Select components //! //! This should be the last operation void select_components(); //! \brief Return true if the component must be removed bool is_component_to_remove(int id) {return m_is_component_to_remove[id];} int nb_component_to_keep() {return m_nb_component_to_keep;} - std::shared_ptr data; - std::vector m_is_component_to_remove; int m_nb_component_to_keep; }; } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_SELECTOR_HPP diff --git a/src/database/switch_basis.cpp b/src/database/switch_basis.cpp index bdf1742..7e19718 100644 --- a/src/database/switch_basis.cpp +++ b/src/database/switch_basis.cpp @@ -1,76 +1,82 @@ /*------------------------------------------------------- - Module : database - File : switch_basis.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "switch_basis.hpp" #include +#include "canonicalizer.hpp" namespace specmicp { namespace database { void BasisSwitcher::change_basis(std::vector& new_basis) { + if (not is_database_canonical()) + { + canonicalize_database(data); + } + Eigen::MatrixXd beta = Eigen::MatrixXd::Zero(data->nb_component, data->nb_component); Eigen::VectorXd kswap(data->nb_component); // swap coefficients, build beta matrix and kswap for (int i=0; inb_component; ++i) { if (new_basis[i] >= data->nb_component) { int j = new_basis[i] - data->nb_component; beta.row(i) = data->nu_aqueous.row(j); data->nu_aqueous.row(j) = Eigen::VectorXd::Zero(data->nb_component); data->nu_aqueous(j, i) = 1.0; kswap(i) = data->logk_aqueous(j); data->logk_aqueous(j) = 0.0; } else { beta(i, i) = 1.0; kswap(i) = 0.0; } } Eigen::MatrixXd betainv = beta.inverse(); data->nu_aqueous = data->nu_aqueous * betainv; data->logk_aqueous -= data->nu_aqueous * kswap; swap_aq_param(new_basis); swap_labels(new_basis); // now do the minerals data->nu_mineral = data->nu_mineral*betainv; data->logk_mineral -= data->nu_mineral * kswap; data->nu_mineral_kinetic = data->nu_mineral_kinetic*betainv; data->logk_mineral_kinetic -= data->nu_mineral_kinetic * kswap; } void BasisSwitcher::swap_aq_param(std::vector new_basis) { for (int i=0; i< data->nb_component; ++i) { if(new_basis[i] >= data->nb_component) { const int j = new_basis[i]; data->param_aq.row(i).swap(data->param_aq.row(j)); } } } void BasisSwitcher::swap_labels(std::vector new_basis) { for (int i=0; inb_component; ++i) { if (new_basis[i] >= data->nb_component) { std::swap(data->labels_basis[i], data->labels_aqueous[new_basis[i]-data->nb_component]); } } } } // end namespace database } // end namespace specmicp diff --git a/src/database/switch_basis.hpp b/src/database/switch_basis.hpp index 305904c..8982fbf 100644 --- a/src/database/switch_basis.hpp +++ b/src/database/switch_basis.hpp @@ -1,49 +1,46 @@ /*------------------------------------------------------- - Module : database - File : switch_basis.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_DATABASE_SWITCHBASIS_HPP #define SPECMICP_DATABASE_SWITCHBASIS_HPP //! \file switch_basis.hpp Switch basis -#include -#include "data_container.hpp" +#include "module.hpp" namespace specmicp { namespace database { -class BasisSwitcher +class BasisSwitcher: public DatabaseModule { public: - BasisSwitcher(std::shared_ptr thedata): - data(thedata) + BasisSwitcher(RawDatabasePtr thedata): + DatabaseModule(thedata) {} //! \brief Change the basis //! //! @param new_basis list of id of the new basis //! //! The new basis is a list of id, id = id_component for no swapping //! or id = id_aqueous + nb_component for swapping a secondary species void change_basis(std::vector& new_basis); private: //! \brief Swap aqueous parameters during a basis transformation void swap_aq_param(std::vector new_basis); //! \brief Swap labels - called during a basis transformation void swap_labels(std::vector new_basis); - - std::shared_ptr data; }; } // end namespace database } // end namespace specmicp #endif // SPECMICP_DATABASE_SWITCHBASIS_HPP